{'fraction_non_alphanumeric': 0.03788636218596439, 'fraction_numerical': 0.01133825437682146, 'mean_word_length': 4.268997982515131, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improved Perturbative QCD Approach to the Bottomonium Spectrum

Recently it has been shown that the gross structure of the bottomonium spectrum is reproduced reasonably well within the non-relativistic boundstate theory based on perturbative QCD. In that calculation, however, the fine splittings and the S-P level splittings are predicted to be considerably narrower than the corresponding experimental values. We investigate the bottomonium spectrum within a specific framework based on perturbative QCD, which incorporates all the corrections up to O(alpha_S^5 m_b) and O(alpha_S^4 m_b), respectively, in the computations of the fine splittings and the S-P splittings. We find that the agreement with the experimental data for the fine splittings improves drastically due to an enhancement of the wave functions close to the origin as compared to the Coulomb wave functions. The agreement of the S-P splittings with the experimental data also becomes better. We find that natural scales of the fine splittings and the S-P splittings are larger than those of the boundstates themselves. On the other hand, the predictions of the level spacings between consecutive principal quantum numbers depend rather strongly on the scale mu of the operator \propto C_A/(m_b r^2). The agreement of the whole spectrum with the experimental data is much better than the previous predictions when mu \simeq 3-4 GeV for alpha_S(M_Z)=0.1181. There seems to be a phenomenological preference for some suppression mechanism for the above operator.


Introduction
For a long time most successful theoretical approaches to study the heavy quarkonium spectra have been those based on various phenomenological potential models. These model approaches have been capable not only of reproducing the charmonium and bottomonium spectra to a high accuracy, but also of explaining various other properties of heavy quarkonia such as their transition rates and decay rates. Through this success, these phenomenological models have established, essentially, that the heavy quarkonium states can be described well as non-relativistic boundstate systems; see e.g. [1] for one of the most recent analyses. On the other hand, the problem of the phenomenological approaches is that it is difficult to improve the theoretical predictions systematically, and that it is difficult to relate the parameters of the models to the fundamental theory.
Recently there has been new progress in explaining these heavy quarkonium spectra within the framework of non-relativistic boundstate theory based on perturbative QCD. It has been shown that, by incorporating the cancellation of O(Λ QCD ) renormalons contained in the pole mass and the static QCD potential, the gross structure of the bottomonium spectrum is reproduced reasonably well by the perturbative computation of the spectrum up to O(1/c 2 ) = O(α 4 S m b ) [2,3]. Furthermore, it has been shown that the static QCD potential, calculated in a series expansion in α S up to O(α 3 S ), agrees well with typical phenomenological potentials in the region relevant to bottomonium spectroscopy, once the O(Λ QCD ) renormalon contained in the QCD potential is cancelled against that contained in the b-quark pole mass [4,5]. Since the static QCD potential calculated by lattice simulations is consistent with the phenomenological potentials in this region [6,7], all these calculations are consistent with one another. See also [8], which made a direct comparison of lattice results and perturbative predictions of the QCD potential in a renormalon-subtracted scheme and found the same consistency.
The key concept which led to these new results from perturbative QCD, is the following. Conventionally, theoretical calculations of the energy levels of a non-relativistic quarkantiquark boundstate closely followed that of a QED boundstate such as positronium: it starts from the natural picture that, when an electron and a positron are at rest and far apart from each other, they tend to be free particles and the total energy of the system is given by the sum of the energies of the two particles (pole masses); as the electron and positron approach each other, the energy of the system decreases due to the negative potential energy, so that the total energy of the boundstate is given as the sum of the pole masses minus the binding energy. When the calculation along the same line was applied to the quark-antiquark system, however, the perturbative expansion of the boundstate energy turned out to be poorly convergent, due to the contributions from infrared (IR) gluons with wave-lengths of order Λ −1 QCD . We can regard this as reflecting the invalidity of the free quark picture when the quark and antiquark are far apart from each other. On the other hand, intuitively we expect that there should be a way to calculate the boundstate energy in which the contributions of IR gluons can be mostly eliminated. This is because when the boundstate size is sufficiently smaller than Λ −1 QCD , IR gluons cannot resolve the colour charges of the constituent particles, so that they decouple from this colour-singlet system. Indeed this idea was theoretically validated in the language of renormalons and their cancellation [9,10]. As a result, the convergence of the perturbative expansion improved dramatically, extending the predictive power of perturbative QCD beyond what could be achieved before.
Nevertheless, there remain some problems regarding the above prediction of the bottomonium spectrum from perturbative QCD. Among them, especially interesting is the problem that the fine structure in the 1P j levels as well as the splittings between the 2S and 1P j states are predicted to be considerably narrower than the corresponding experimental values. One may think that the level of agreement of the theoretical prediction with the experimental data is still consistent within errors: according to an estimate based on next-to-leading order renormalons, each energy level has a theoretical uncertainty of order Λ 3 QCD · a 2 (a is the typical size of the corresponding boundstate) which may be comparable in size to the 2S-1P j splittings and may be much larger than the fine splittings in the 1P j levels; compare the error estimate in [3]. One should note, however, that if we calculate these level splittings instead of the individual energy levels, the O(Λ 3 QCD ) renormalons should get largely cancelled when we take the differences of the energy levels. Hence, the theoretical uncertainties of the splittings can be much smaller than those of the individual energy levels, and the disagreement between the theoretical prediction and the experimental data may be quite serious. It is important to clarify whether it is possible to resolve these disagreements in the fine splittings and the S-P splittings within the context of perturbative QCD, e.g. by including higher-order corrections, or whether we need to take into account specific non-perturbative contributions for this purpose.
In this paper we investigate this problem of the fine splittings and the S-P splittings in the bottomonium spectrum within a specific framework based on perturbative QCD. We note that these splittings have been successfully reproduced by the phenomenological potential models, and that a connection between the static QCD potential and phenomenological potentials has been elucidated in [4,5]. In order to take advantage of these results, we develop a framework which enables detailed comparison of the predictions of the phenomenological models and of perturbative QCD. We also incorporate some of the higher-order corrections to the nonrelativistic Hamiltonian of the quark-antiquark system which have not been included in the analysis [3].
The basic theoretical ingredients of our analysis are as follows: (1) We take into account the cancellation of the O(Λ QCD ) leading renormalons by reexpressing the b-quark pole mass in terms of the MS mass. (2) We take a specific scheme for the perturbative expansion such that all the corrections up to O(α 4 S m b ) and O(α 5 S m b ) are incorporated in the calculation of the S-P splittings and the fine splittings, respectively. Furthermore, some of the higher-order corrections, which appear to be important for these observables, are incorporated.
The organization of the paper is as follows. In Sec. 2 we present the framework of our calculation. We examine the energy levels and the wave functions of our zeroth-order Hamiltonian in Sec. 3. The analysis of the fine splittings is given in Sec. 4 and that of the S-P level splittings in Sec. 5. Then, in Sec. 6, we compare the whole structure of the bottomonium spectrum given by our prediction, by other theoretical predictions and by the experimental data. Concluding remarks are given in Sec. 7. We derive a formula useful for our analysis in the Appendix. We first recall the non-relativistic Hamiltonian of a quark and antiquark pair given in the series expansion in 1/c up to O(1/c 2 ) which is determined from perturbative QCD. (See e.g. [11,12,13,14].) Considering an application to the bottomonium states, we assume the quark (antiquark) to be the b-quark (b-quark). The Hamiltonian is given by We choose the zeroth-order part of the Hamiltonian to be where m b is the pole mass of the b-quark, and V QCD (r) denotes the static QCD potential up to O(α 3 S ). This choice differs from the usual zeroth-order Hamiltonian of the 1/c-expansion, since terms of the QCD potential. Other operators of Eq. (1) are treated as perturbations to H 0 , all of which are O(1/c 2 ) in the usual order counting in 1/c-expansion. U + W A constitutes the O(1/c 2 ) part of the Breit Hamiltonian known from QED boundstate theory, where the spin-dependent operator is given by and the spin-independent operator is given by On the other hand, represents the operator characteristic to the non-abelian gauge theory. In this paper, unless the argument is specified explicitly, α (n l ) S denotes the strong coupling constant renormalized at the renormalization scale µ, defined in the MS scheme with n l active flavours, i.e. α (n l ) S ≡ α (n l ) S (µ); C F = 4/3 and C A = 3 are the colour factors; L and S are, respectively, the orbital-angular momentum and the total spin of the quark-antiquark pair. For the bb system, n l = 4. PSfrag replacements linear extrapolation r UV = 0.5 GeV −1 r IR = 4.5 GeV −1 Figure 1: Construction of the improved potential.

Improved Potential
In our analysis of the bottomonium spectrum and wave functions, we use an improved "potential" E imp (r) instead of 2m b + V QCD (r) in the zeroth-order Hamiltonian H 0 . This E imp (r) is constructed in the following way: We divide the range of r into three regions by introducing ultraviolet (UV) and infrared (IR) distance scales, r UV and r IR , see Fig. 1.
(i) At intermediate distances, r UV < r < r IR , where the bulk of the bottomonium wave functions are located, E imp (r) is identified with the total energy of the static bb system E bb tot (r)| µ=µ 2 (r) = [2m b + V QCD (r)]| µ=µ 2 (r) computed in [5]. E bb tot (r) is defined in Eq. (5) of that paper and depends on the parameters α In identifying E imp (r) with E bb tot (r), the scale µ is determined as a function of r according to the second prescription (µ = µ 2 (r), Eq. (14)) in [5]: since E bb tot (r) is less µ-dependent and its series expansion converges better if we choose a larger value for µ when r is smaller, and if we choose a smaller value for µ when r is larger, we consider our choice of µ to give a more accurate prediction for E bb tot (r) than choosing some fixed (r-independent) value of µ; see Figs. 2,3. The c-quark MS mass is taken as m c = 1.243 GeV [2]. We will explain how we fix m b in our analysis below. For other details, we follow the convention of Secs. IIA and IIB of [5].
(ii) At short distances, r < r UV , we use a renormalization-group improved QCD potential.
It is obtained by integrating the three-loop renormalization-group improved interquark force F (r) ≡ −dV QCD (r)/dr, following the method of Sec. 4 of [4]. There, it was shown PSfrag replacements PSfrag replacements for N = 1, 2, 3 when µ is fixed independently of r: a) µ = 1 GeV, and b) µ = 2 GeV. The scale µ that provides the best convergence of the perturbative series depends on r. that the QCD potential becomes more convergent if we improve the interquark force by means of the renormalization group and integrate it over r, rather than directly improving the QCD potential by the renormalization group. * The initial value for the renormalization group evolution of F (r) and the constant part of E imp (r) are determined such that E imp (r) becomes continuous at r = r UV up to the first derivative.
(iii) At long distances, r > r IR , we use a linear potential, E imp (r) = C 1 r + C 2 , where C 1 and C 2 are determined such that E imp (r) becomes continuous at r = r IR up to the first derivative.
The main part of the improved potential E imp (r) is that given in the intermediate-distance region (i). Indeed this part of the potential dictates the main features of the results of our analysis. According to its construction, however, E imp (r) as defined in the region (i) becomes unstable and unreliable at r < ∼ 1/m b and r > ∼ 1/Λ QCD . This is the reason why we continue it to other definitions at short and long distances. The short-and intermediate-distance parts (i) and (ii) are determined by perturbative QCD, whereas the long-distance part (iii) is not. We will show that the shape of E imp (r) in the long-distance region (iii) does not affect the results of our analysis significantly. It should be noted that, with our definition, the µ-dependence of E imp (r) has been completely eliminated. The µ-dependences mentioned later in this article are those which stem from the other terms of the Hamiltonian (see below). We urge the reader to consult [4,5] for the physics background and detailed features of E imp (r) as defined in the regions (i) and (ii). In our analysis below, we set r UV = 0.5 GeV −1 and r IR = 4.5 GeV −1 .

Perturbative Expansion
We solve the non-relativistic Schrödinger equation numerically with the improved potential, and determine the zeroth-order energy level E (0) ψ and wave function | ψ of a quarkonium state. We treat U, W A and W N A as perturbations and calculate the perturbative corrections to the energy level by We will also examine corrections induced by some of the O(1/c 3 ) operators: where γ E = 0.5771... is the Euler constant; T R = 1/2. The spin-dependent operators δU LS and δU S were derived in [15,16]. (Earlier incomplete results can be found in [17,11].) We derived δW N A from the result of [18] in the following way: We discard the logarithm originating from the IR divergence associated with the energy scale. This gives the first term of Eq. (12). We have to take into account the contributions which come from the unitary transformation, The second term of Eq. (12) is generated by this unitary transformation. We note that the non-logarithmic part of δW N A cannot be determined unambiguously, since it mixes with other O(1/c 3 ) operators through IR divergence. Our definition merely represents one possible scheme.
Only when we add all the contributions to the spectrum at O(1/c 3 ), the sum is free from IR divergence and can be defined unambiguously. Since at present we do not know the full form of the Hamiltonian up to O(1/c 3 ), ‡ this problem cannot be circumvented in any case. † The operator representation of X reads operators of the Hamiltonian have been identified in Ref. [16]. We have to specify how we treat the pole mass m b in the operators p 2 /m b , U, W A , W N A , δU and δW N A . We express the pole mass m b in terms of the MS mass m b in the series expansion in α (7) and (8) of [5]. § After that we can, in principle, reexpand these operators in α (4) S (µ), since the pole mass enters the denominators and ℓ m . There is, however, no known guiding principle how to organize these expansions, since such reexpansions cannot be carried out consistently with the expansions in 1/c. ¶ This is in contrast to the reexpansion of E bb tot (r) = 2m b + V QCD (r), which has a guiding principle by the cancellation of O(Λ QCD ) renormalons, although the reexpansion is indeed inconsistent with the 1/c-expansion. Hence, we keep the pole mass as a function of m b , α   4.190 GeV. We will examine how the uncertainty of the pole mass of order Λ QCD affects our predictions in Sec. 6.
We take the input for the strong coupling constant as α S (µ 2 (r)) for E imp (r) and α (4) S (µ) for the other operators. We evolve the coupling by solving the 3-loop renormalization-group equation numerically and match it to the 4-and 3-flavour couplings successively through the matching condition [24]. (Although the 4-loop running of the MS coupling constant is available, we consider the 3-loop running more consistent in our analysis, which incorporates corrections up to the 2-loop finite part of V QCD (r).) Our choice of the zeroth-order Hamiltonian and the way we organize the perturbative expansion is largely motivated by the success of phenomenological potential models. In fact, the above organization of the perturbative expansion follows, to a large extent, the approaches of phenomenological model analyses, if we identify our H (imp) 0 with the non-relativistic Hamiltonian, p 2 /m + V pheno (r), used in those analyses. Since E imp (r) agrees well with typical phenomenological potentials up to an additive constant, we expect that we can make close comparisons with phenomenological model analyses, and that eventually the bottomonium spectrum may be reproduced with a good accuracy. Note that our zeroth-order quarkonium wave functions, which are determined from H (imp) 0 , include some of the higher-order corrections in the usual order counting of the 1/c-expansion, since E imp (r) includes the O(1/c) and O(1/c 2 ) corrections to § These formulas were derived originally in [19,20,21].
¶ For instance, if we organize the perturbative expansions of the energy levels appropriately, renormalons contained in them are of order Λ 4 QCD , whereas if one expands the levels in 1/c, renormalons become order Λ 3 QCD and worsen the perturbative convergence [22]. This may be regarded as an explicit example of problems due to an inconsistency between 1/c power counting and renormalon cancellations (beyond those controlled by ǫ-expansion [23]).
The matching scales are taken as m b and m c , respectively.   the static QCD potential. By the same token, the zeroth-order energy levels are different from the Coulomb energy levels, and in particular they depend on the orbital angular momentum l.

Zeroth-order Energy Levels, Wave Functions and Scales
We show numerical solutions to the zeroth-order Schrödinger equation (7). Here and hereafter we take the input value for the strong coupling constant to be the present world average value α ψ vary by about 10 MeV (20 MeV) when µ is varied from 2 to 5 GeV (1 to 2 GeV). Nonetheless, if we take the difference of any of the two energy levels, the µ dependences cancel mostly. The µ-dependences of the wave functions are fairly weak.
In [2,3] the scale µ = µ ψ for each quarkonium state | ψ was fixed by minimizing the µdependence of each energy level calculated in a fixed-order perturbative expansion (minimalsensitivity prescription). The scale fixed in this way turns out to represent the physical size of the corresponding quarkonium state fairly well. This was shown in [2] by comparing the scale µ ψ and the support function defined by where R ψ (r) is the radial part of the wave function for the state | ψ . The support function represents the support in the momentum-space integral in the calculation of (the major part of) the energy level. The Coulomb wave function, evaluated with α S (µ ψ ), was used to compute the support function f ψ (q); the charm mass effects in loops were not taken into account. Here we compare the support function computed with our zeroth-order wave function and the scale fixed in [3], both of which include the charm mass effects. The resulting support functions for the S states are shown in Fig. 5 and for the P states in Fig. 6. We see that with respect PSfrag replacements 1P 2P Figure 4: r 2 |R ψ (r)| 2 for different values of µ = 1, 2, 3, 4 and 5 GeV, where R ψ (r) denotes the radial part of the zeroth-order energy eigenfunction. The area below each curve is normalized to unity. The one curve that visibly differs from the others corresponds to µ = 1 GeV. Note that the scales differ between the plots.
s (µ) is also plotted (dashed curve). Since the analysis that we advocate in this work does not attribute scales to the individual states, the scales indicated by the dotted lines are taken from [3], Table II  to the treatment in [2], the support functions are shifted slightly towards higher momentum. This is reasonable, since the wave functions calculated here are peaked closer to the origin in coordinate space than the Coulomb wave functions (see e.g. Fig. 7 below). Since the scales µ ψ are located within the IR "tails" of the corresponding support functions f ψ (q), we confirm that the above interpretaion of the minimal-sensitivity scale µ ψ is valid also when we include the charm mass effects and the O(1/c) and O(1/c 2 ) corrections of V QCD (r) to the wave functions.

Fine Splittings
In this section we examine the fine splittings in the bottomonium spectrum within our framework and compare them with the experimental data as well as with previous theoretical predictions in the literature. Within perturbative QCD, it is expected that in principle the fine splittings can be computed much more accurately than the individual energy levels. This is because the potentials which contain order Λ 3 QCD r 2 renormalons (the static QCD potential and W N A ) do not induce the fine splittings, so that these renormalons cancel in the computation of the fine splittings.
The fine splittings (or fine structure) are the level splittings among the states with the same principal quantum number n, orbital-angular momentum l > 0 and spin s but with different total angular momentum j. Experimentally, the fine splittings have been observed among the 1 3 P j [(n, l, s) = (2, 1, 1)] and 2 3 P j [(n, l, s) = (3, 1, 1)] states, we therefore examine the predictions for these splittings.
In perturbative QCD the fine splittings are induced first at O(α 4 S m b ) = O(1/c 2 ) by the operator U: where This is the fixed-order formula used in [2,3]. In this formula, the Coulomb wave function | ψ C is used to compute the expectation value. The scale dependence of ∆E 4 : three powers of α S come from the Coulomb wave function ψ C | × | ψ C , and one power comes from the operator U. In [2,3] the scale µ is fixed by the minimal-sensitivity prescription; c.f. the previous section.
In our approach, we calculate the fine splittings from ∆E U = ψ | (U +δU) | ψ . As compared to the fixed-order formula, some of the higher-order corrections are incorporated through the wave function | ψ and the operator δU.
In phenomenological approaches, one computes the fine splittings using the same operator U but using the wave functions determined from phenomenological potentials: Some of the higher-order corrections to U constrained by the Gromes relation [26] have also been incorporated. As an example we compare our results with those of a phenomenological model with a Coulomb-plus-linear potential (Cornell potential) [27]: with κ = 0.52 and a = 2.34 GeV −1 . Tab. 3 compares the expectation values of the operators U LS and U S with respect to the Coulomb wave function | ψ C , our zeroth-order wave function | ψ and the Cornell wave function | ψ Cornell . The operators U LS and U S induce the fine splittings through X LS U LS − D S U S , where (X LS , D S ) = (−2, −2), (−1, 1), (1, −1/5) for j = 0, 1, 2, respectively, and (l, s) = (1, 1). The fine splittings between adjacent levels are therefore U LS −3 U S for P 1 −P 0 and 2 U LS +6 U S /5 for P 2 − P 1 . The values of α      in [3]. Taking into account the numerical values from the table, U LS and U S give roughly the same contribution to the P 1 -P 0 -splitting, while the P 2 -P 1 -splitting is dominated by U LS . We see that the expectation values of U LS and U S with respect to our wave function | ψ and the Cornell wave function | ψ Cornell are very much larger than the expectation values with respect to the Coulomb wave function | ψ C . * The reason for this behaviour can be understood in the following way: Since the potentials E imp (r) and V Cornell (r) are steeper (i.e. the attractive forces are stronger) than the Coulomb potential in the intermediate-distance region, † the wave functions are more centered towards the origin for | ψ and | ψ Cornell than for | ψ C ; see Fig. 7 which compares the squared radial wave functions multiplied by the phase space factor for the 1P states. Therefore the wave functions r ψ and r ψ Cornell are enhanced close to the origin * This is not necessarily true for the 2P states with α S = 0.726, but this can be regarded as originating from another effect, which we explain below. Namely, the value α S = 0.726 is unrealistically large.
† The cancellation of O(Λ QCD ) renormalons suggests that this behaviour can be understood naturally in terms of the QCD force F (r) = −V ′ QCD (r) = −C F α F (1/r)/r 2 : F (r) becomes more attractive than the Coulomb force as r increases due to the running of the F -scheme coupling constant α F (1/r) [4,8]. as compared to the Coulomb wave function. The enhancement factors turn out to be large for the P -wave states. ‡ Since the expectation values of U LS and U S are determined mostly by the short-distance behaviour of the wave functions, they are enhanced by large factors. This feature can be verified in Fig. 8, where we compare the integrands when the expectation values of U LS are expressed as integrals over r; the main contributions come from distances r < ∼ 2 GeV −1 . We should stress that the enhancement of the wave functions originates mainly from the behaviour of the potentials in the intermediate-distance region and not from the shortdistance behaviour. This can be seen from the fact that the difference between ψ | U (LS,S) | ψ and ψ Cornell | U (LS,S) | ψ Cornell is much smaller than the difference between ψ | U (LS,S) | ψ and ψ C | U (LS,S) | ψ C : the essential difference between our potential E imp (r) and the Cornell potential resides in the short-distance region.
We include also the contribution of the O(1/c 3 ) operator δU into our prediction of the fine splittings. As far as we know, there is no other operator which contributes to the fine splittings at O(1/c 3 ). Considering that our wave function | ψ includes all the next-to- . Inclusion of the operator δU reduces the scale dependence of our prediction. In Tab. 4 we compare the fine splittings calculated from the matrix elements ψ | U | ψ and from ψ | U + δU | ψ . The former depend ‡ By way of example, if we squeeze the function r n such that it takes a same value at half distance, i.e. r n → (2r) n , then the enhancement factor becomes 2 n . Thus, the enhancement factor is larger for a larger n. For the P -wave states, n = 4, because 2 powers come from the wave function squared and 2 powers come from the phase space.  on µ rather strongly and are larger for smaller µ, since U is proportional to α S (µ). § We see that the scale dependence has decreased considerably by the effect of δU. The scale dependences become minimal at µ ≃ 2 GeV for both 1P j and 2P j states. We may try to reproduce the splittings calculated with ψ | U + δU | ψ at these scales by taking an appropriate choice of the scales in the lower-order predictions (the splittings calculated with ψ | U | ψ ). Then the scales become µ ≃ 2 GeV for both the 1P j and 2P j states. All these scales are larger than the scales chosen for the respective states in [3], which are fixed by minimizing the scale dependence of the individual energy levels. (µ ≃ 1 GeV for the 1P j states and µ ≃ 0.8 GeV for the 2P j states; see Fig. 6.) This feature is consistent with a naive expectation: we would expect that the latter scales represent the typical scales of the binding energies, or the inverse of the sizes of the boundstates (Fig. 6), whereas the former scales represent those probed by the operators U or U + δU, which are larger because the contributions to the matrix elements come from shorter distances (Fig. 8).
To our knowledge, so far there has been no systematic argument on the order of renormalons contained in the fine splittings. Naively we expect that the largest renormalon contained in the calculation of the fine splittings would be of order Λ 3 QCD /m 2 b ≃ 1-10 MeV. This estimate is based on the Gromes relation [26] which tells us that a part of the operator U LS is determined from the static QCD potential: 26) and the fact that dV QCD (r)/dr contains an O(Λ 3 QCD r) renormalon. We may take this as an order of magnitude estimate of uncertainties of our present predictions for the fine splittings. Compared with this error estimate, our predictions of the fine splittings calculated with ψ | U + δU | ψ in Tab. 4 are in reasonable agreement with the experimental data.
We are now able to interpret the reasons why the fine splittings of the energy levels computed in [3] turned out to be quite small for the 1P j states but not so much for the 2P j states. The first point to note is that large enhancement factors are generated by the fact that the quarkonium wave functions are more centered toward the origin if we solve the Schrödinger equation with   E imp (r) (which we believe to be more realistic) rather than with the Coulomb potential. The second point is that the natural scales to be chosen for evaluating the expectation values of the operator U are larger than the natural scales for the individual energy levels: we confirmed this by incorporating the effects of the higher-order operator δU, and the results are qualitatively consistent with a naive expectation. We find that that the first effect overwhelming the second one resulted in the quite small splittings among the 1P j states in [3]. On the other hand, the cancellation of the first and the second effect resulted in reasonable sizes of the fine splittings for the 2P j states in that paper, which, in the light of our present observations, may be regarded as rather accidental.

S-P Splittings
In this section we examine the splittings between the S-wave and P -wave states. In the Coulomb spectrum the S-wave and the P -wave states with the same principal quantum number n are degenerate. In perturbative QCD, the splittings are induced by V QCD (r) starting from O(1/c) as well as by U, W A , W N A at O(1/c 2 ). Among these operators, V QCD (r) (after cancelling the order Λ QCD renormalon) and W N A contain order Λ 3 QCD r 2 renormalons. Therefore, the order Λ 3 QCD r 2 renormalons do not cancel completely in the perturbative computation of the S-P splittings. Namely, the theoretical uncertainties of the S-P splittings are expected to be larger than those of the fine splittings.
In Tab. 5 we show the S-P splittings of our zeroth-order energy levels E   the expectation values of W A is given in the Appendix.) We expect that the contributions of these operators would be smaller than the E   S (µ). We see that the scale-dependence is reduced by including the effect of the higher-order correction δW N A , although W N A + δW N A is still unstable around µ = 1 GeV. If we choose a relatively large scale µ ≃ 3 GeV, the contributions of the higher-order corrections, δU and δW N A , become small and the scale dependences of the S-P splittings are relatively small; see Tabs. 5 and 6.
Furthermore, as shown in Fig. 9, the expectation values ψ | W N A | ψ and ψ | W N A + δW N A | ψ are dominated by short-distance contributions. Therefore, following the same line of argument as in the previous section, we expect the natural scales for these contributions to be larger than that for the binding energy.
If we compare our predictions for the S-P splittings with the predictions of the fixed-order perturbative expansion in [3], there are two competing effects, just like what we found in the case of the fine splittings: these effects are the difference of the wave functions and the difference of the scales µ in the operator. Consequently our predictions for the S-P splittings turn out to be larger for a smaller scale µ, are typically larger than the fixed-order predictions for the 2S-1P splittings, and are of similar magnitude to the fixed-order predictions for the 3S-2P splittings. Since, however, the S-P splittings are dominated by the O(1/c) correction from V QCD (r) [i.e. the contributions from the O(1/c 2 ) operators are only subleading], the differences from the fixed-order predictions are not as pronounced as in the case of the fine splittings.
As long as we calculate the S-P splittings directly, as done up to here, we see no indication of large theoretical uncertainties. However, when we examine the individual energy levels, some  for the 1P j states. The graphs peak at rather small distances, indicating that the natural scales for the expectation values are larger than those for the boundstate. The agreement of the S-P splittings between the theoretical predictions and the experimental data has improved somewhat from the fixed-order results [3]. It is, however, not as good as one would naively expect. Namely the differences of the splittings between our predictions (for the scale choice µ ≃ 2-5 GeV) and the experimental data are larger than the contributions  of the O(1/c 2 ) operators. As stated in the introduction, it is important to clarify whether the level of disagreement is still consistent within perturbative uncertainties. From the examination of the contributions to the individual energy levels, we conjecture that there would still be large theoretical uncertainties to the S-P splittings, particularly from the contributions of W N A , since at present we have no systematic argument on how much of the O(Λ 3 QCD ) renormalons cancel in the S-P splittings. That is, we consider the cancellation of the large and scale-dependent contributions from W N A in the S-P splittings to be accidental, unless we find a systematic argument in support of it. In this regard, we consider the predictions for the S-P splittings much less reliable than those for the fine splittings.
There are also other indications that the operator W N A becomes a source of instability of theoretical predictions. For instance, ref. [28] addresses W N A to be the source of the large uncertainties of the cross section for e + e − → tt close to threshold. Ref. [29] shows that IR logarithms, which are related to higher-order corrections to W N A and V QCD , generate very large corrections to the bottomonium spectrum at O(α 5 S m b ) and consequently cause a large uncertainty to it. Furthermore, refs. [30,31,32] resum the IR logs; it was found that W N A is reduced due to this resummation in the corresponding Wilson coefficients. Ref. [32] claims that the resummations lead to a stable theoretical prediction for the top threshold cross section. Ref. [33], however, reported a result which disagrees with that of [30,31].
In this paper we discard the IR logs altogether, which appear first at O(α 5 S m b log α S ) in the spectrum. This is not because we consider them unimportant, but because we believe that we are not yet in a position to treat this problem properly. The present treatment of the IR logs in [18], [30]- [33] seems to comprise following unsatisfactory aspects in the calculations of the heavy quarkonium spectrum and the top threshold cross section (besides the disagreement). Since log α S is not particularly large, it appears that the non-logarithmic terms cannot be neglected in comparison to the IR logs. These non-logarithmic terms are not yet fully known at O(α 5 S m b ) and beyond, so that we are unable to incorporate them unambiguously. Furthermore, if a resummation of the IR logs stabilizes the theoretical predictions substantially, we would like to understand the physical meaning behind it, i.e. why the resummation is important. A hint to these questions was suggested in [22]: it shows that O(Λ 3 QCD ) renormalons will be suppressed if we incorporate the offshellness of the quarks, which acts as an IR cutoff in the temporal dimension. This applies to the O(Λ 3 QCD ) renormalons contained in W N A and V QCD (and their  Table 9: Energy levels for all states including E (0) ψ , U +δU , W A , and W N A +δW N A . The parameters used are α higher-order corrections), and they are closely related to the IR logs and their resummation. We suspect that further investigations of all these problems may be a way to clarify and solve the problem of the operator W N A we face here.

The Spectrum
In this section we compare the whole bottomonium spectrum as determined experimentally with various theoretical predictions. We list the energy levels numerically in Tab. 9 and show the spectrum in Fig. 10. The levels were calculated according to the framework explained in Sec. 2, using the input parameter α  4.190 GeV to make the the prediction for the 1S state coincide with the experimental value for our favoured scale of µ = 3 GeV. ¶ . We consider this different choice of m b a more natural way to achieve coincidence with the experimental value than a simple shift of the whole spectrum; the numerical difference between these two prescriptions is actually very minor.
In Fig. 10 we compare our results to the spectrum obtained from the experimental data ("Exp"), the predictions of the fixed-order perturbative expansions ("BSV"), and the result of [1] as a typical prediction of recent phenomenological models ("Eichten-Quigg").
One can verify the conclusions of our analysis in the previous two sections: The fine splittings and the S-P splittings are larger than those of the fixed-order results for the n = 2 states, whereas they are of similar magnitude to the fixed-order results for the n = 3 states. The scaledependence of our predictions originates mostly from the scale-dependence of the operator W N A + δW N A ; c.f. the discussion in the last section. Only the gross level spacings between ¶ The change of m b by this adjustment does not alter the qualitative features of the predictions for the fine splittings and the S-P splittings discussed in previous sections.
We follow the scheme and the scale-fixing condition A of Sec. 4.2 of [3], except that we use the numerical solution to the renormalization-group equation for the strong coupling constant. The results are obtained for the input parameters α 234 GeV have been used to make the 1S state coincide with experiment for µ = 3 GeV. The columns for µ = 1, 2, 4 and 5 GeV have been shifted to achieve this coincidence. "Exp" shows the experimental values and "Eichten-Quigg" those obtained in [1]. Finally, "BSV" corresponds to the formalism of [3], the two columns represent choices of α adjacent n's are affected visibly by changes of µ between 2-5 GeV, whereas also the S-P splittings vary visibly for a smaller µ between 1-2 GeV. The size of this variation is very large considering that the scale dependence is formally an O(1/c 3 ) = O(α 5 S m b ) effect. The level spacings between consecutive n's as well as the S-P splittings increase for smaller µ. We regard the large scale-dependence generated by W N A and its higher-order corrections as the largest theoretical uncertainty of our prediction.
Let us note the effects of the operator W N A or W N A + δW N A in particular: the level spacings between consecutive n's are increased, while the S-P splittings are reduced. The reason can be understood as follows. The operator W N A (+δW N A ) generates an attractive potential proportional to 1/r 2 (with a logarithmic correction) which is particularly strong at short distances. Hence, those states which have larger wave functions close to the origin acquire larger binding energies. Since the states with lower n have larger amplitude close to the origin and therefore acquire larger binding energies, the level spacings between the adjacent n's become wider. Since the S states acquire larger binding energies than the P states, the S-P splittings becomes narrower.
Generally, our prediction of the spectrum has a better agreement with the experimental data than the fixed-order results. The agreement seems to be better for a larger scale choice, which appears reasonable, since the natural scale of the operator W N A + δW N A would be large. (The scale dependences due to other effects are much smaller.) The agreement seems to be optimal for a scale choice µ ≃ 3-4 GeV. We also examined our predictions for different values of the input α (5) S (M Z ) within the present world-average values 0.1181 ± 0.0020 [25]. We find that generally all the level spacings and splittings become larger for larger α (5) S (M Z ), since the binding energy increases. The widening of the level spacings, however, can be compensated largely by choosing a larger value for µ. The spectrum of the phenomenological model still has a better agreement with the experimental data, since it includes more parameters which can be adjusted.
In general one should carefully take into account theoretical uncertainties when comparing the whole spectrum with the experimental data. Based on the renormalon argument, each energy level has (at least) an uncertainty of order Λ 3 QCD r 2 , where r should be taken as a typical size of the quarkonium state. Theoretical uncertainties contained in E imp (r) in the region r < r IR can be represented typically by these renormalon estimates [4,5]. Ref. [3] estimated the uncertainties to be ±(5-30) MeV for the 1S state, ±(20-130) MeV for the n = 2 states, and ±(40-220) MeV for the n = 3 states. On the other hand, the level spacings (splittings) have smaller theoretical uncertainties since these theoretical uncertainties cancel, at least partly, as we discussed in the previous two sections.
Let us estimate the errors of our prediction from other sources. In the left part of In the right part of Fig. 11 we also show the effect of a variation of the long-distance part PSfrag replacements  Fig. 10. The scale dependence is decreased and the agreement with the experimental spectrum is strongly improved with respect to Fig. 10, indicating a phenomenological preference for a suppression mechanism for W N A . of the potential E imp (r) in the zeroth-order Schrödinger equation. We vary the tangent of the linear potential at r > r IR by factors of 2 and 1/2. (The first derivative of the potential then becomes discontinuous at r = r IR .) The variations of the energy levels are of order ±(2-5) MeV for the 2P states and ±15 MeV for the 3S state, while they are smaller than 0.1 MeV (and therefore not visible in Fig. 11) for the lower states. This can be easily understood because only the states with principal quantum number n = 3 have a wave function that extends to large enough distances to probe the potential in this region.
For the sake of comparison, we show in Fig. 12 the predictions for the bottomonium spectrum when we set the W N A + δW N A to 0 artificially (in this figure we use α . We see that the scale-dependence reduces (as expected) and there is a much better agreement with the experimental data. Phenomenologically this may be taken as an indication that the contributions of W N A should be suppressed by some mechanism.

Conclusions
We have examined the bottomonium spectrum within a specific framework based on perturbative QCD. The computation of the individual energy levels includes all the effects up to O(α 4 S m b ) and that of the fine splittings contains all the effects up to O(α 5 S m b ). We have also included some important higher-order corrections to the quarkonium wave functions and the energy levels through our use of E imp (r), which is the characterizing feature of our analysis. We have eliminated the scale-dependence from E imp (r) by fixing a specific scheme and we have discussed scale-depedences originating from the other terms of the Hamiltonian. The agreement of the fine splittings among the 1P j states between the theoretical prediction and the experimental data improved drastically as compared to the fixed-order prediction of [2,3]. We find that the centering of the wave functions towards the origin as compared to the Coulomb wave functions, due to the strong attractive force in the intermediate-distance region, strongly enhances the fine splittings. We also find that, in accord with a naïve expectation, the natural scales of the fine splittings are larger than those of the boundstates themselves; the latter were used in the analyses of [2,3]. The predictions for the fine splittings are stable against the variation of the scale µ and are in reasonable agreement with the experimental data both for the 1P j and 2P j states. We also examined the S-P level splittings. The agreement with the experimental data has improved as compared to the fixed-order results, but the predictions are still somewhat smaller than the experimental values. The predictions for these splittings are stable against the variation of µ between 2-5 GeV but become unstable for lower scales between 1-2 GeV. Natural scales of the S-P splittings are also found to be larger than those of the boundstates. On the other hand, the predictions of the level spacings between the adjacent n's depend rather strongly on the scale µ. This stems from the large scale dependence of the operator W N A and must be regarded as a major source of uncertainties in our predictions. We are motivated to choose a relatively large value for the scale µ in view of the dominance of short-distance contributions to ψ | W N A | ψ . (Other effects are much less scale dependent.) If we choose µ ≃ 3-4 GeV, and for α (5) S (M Z ) = 0.1181 and m b = 4.234 GeV, the agreement between our prediction and the experimental data for the whole bottomonium spectrum is fairly good, and is considerably better than the agreement between the fixed-order prediction and the experimental data. There seem to be some indications, however, that the contribution of the operator W N A reduces the stablity of the theoretical prediction and at the same time worsens the agreement between the prediction and the experimental data. At the present state we consider our predictions to be consistent with the experimental data within theoretical uncertainties.
operator identities to W A : Then one finds Hence, the expectation value can be written as All quantities on the right-hand-side can be evaluated from the radial wave function and the energy eigenvalue, which are obtained by solving the Schrödinger equation numerically.