Full Formula for Heavy Quarkonium Energy Levels at Next-to-next-to-next-to-leading Order

We derive a full formula for the energy level of a heavy quarkonium state identified by the quantum numbers $n$, $\ell$, $s$ and $j$, up to ${\cal O}(\alpha_s^5 m)$ and ${\cal O}(\alpha_s^5 m \log \alpha_s)$ in perturbative QCD. The QCD Bethe logarithm is given in a one-parameter integral form. The rest of the formula is given as a combination of rational numbers, transcendental numbers ($\pi$, $\zeta(3)$, $\zeta(5)$) and finite sums (besides the 3-loop constant $\bar{a}_3$ of the static potential whose full analytic form is still unknown). A derivation of the formula is given.


Introduction
For a long time properties of various hadrons have been studied in order to obtain a better understanding of dynamics of the strong interaction. Among various observed hadrons, the heavy quarkonium states are unique, in the sense that these are the only known individual hadronic states whose properties can be predicted in a self-contained manner within perturbative QCD. Namely, several observables associated with individual heavy quarkonium states (such as energy levels, leptonic decay widths and transition rates) can be computed systematically in expansions in the strong coupling constant α s . Such series expansions make sense, since the large mass of the heavy quarks, m(≫ Λ QCD ), and the color-singlet nature of this bound state restricts the relevant dynamical degrees of freedom of this system to those only in an ultraviolet (UV) region and the asymptotic freedom of QCD designates expansions in a small coupling constant.
To carry out computations of higher-order terms of these expansions, development of effective field theories (EFTs) for heavy quarkonium states, such as potential-NRQCD (pNRQCD) [1] and velocity-NRQCD (vNRQCD) [2], has been crucial. It is also essential that various computational technologies have made rapid progress during the same period.
(See e.g. [3,4].) These theoretical tools enable systematic organization of computations of higher-order corrections by separating different energy scales involved in the computations in a systematic manner [5]. These developments opened up vast applications of perturbative QCD in heavy quarkonium physics. See [6] for comprehensive reviews. Some of more recent results can be found in [7,8,9,10,11,12,13].
Among the observables of heavy quarkonium states, it turns out that the energy levels can be predicted particularly accurately in perturbative QCD. The decoupling of infrared (IR) degrees of freedom from these states can be naturally incorporated in computations of the energy levels by using a short-distance mass of the heavy quarks instead of the pole mass, and this prescription improves convergence of the perturbative series of the energy levels drastically [14].
Spectroscopy of heavy quarkonium states (in particular the bottomonium states) has provided an important testing ground of perturbative QCD. The full O(α 4 s m) and O(α 4 s m log α s ) corrections to the energy levels were computed in [15]. Analyses of the bottomonium spectrum, which incorporate the IR decoupling 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 [16]. 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.
Stability of the predictions for the energy levels and their agreement with experimental data are predominantly determined by the prediction for the static energy E stat (r) = 2m pole + V QCD (r). Once the pole mass is expressed in terms of a short-distance mass, the prediction for E stat (r) in perturbative QCD agrees with lattice computations or typical phenomenological potentials in the relevant distance range [17,18]. Furthermore, by increasing the order of the perturbative expansion, the range of r, where convergence and agreement are seen, becomes wider [19].
After the development of EFTs, computations of the O(α 5 s m) and O(α 5 s m log α s ) corrections to the energy levels made progress. Within pNRQCD, the corrections consist of two parts, the contributions from the potential region and ultra-soft region. The nextto-next-to-next-to-leading order (NNNLO) Hamiltonian, which dictates the contributions from the potential region, was computed in [20] (besides the 3-loop corrections to V QCD (r), a 3 , which were computed later numerically in [21,19,22]). The contributions from the ultra-soft 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 one-parameter integral of elementary functions [23]. 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 [24], while the O(α 5 s m) and O(α 5 s m log α s ) corrections for a general Swave state (n, j) were computed in [25]. (See also [26] for earlier computations of the O(α 5 s m) corrections.) In this paper, we derive a full formula for the energy level of a heavy quarkonium state identified by the quantum numbers n, ℓ, s and j, up to O(α 5 s m) and O(α 5 s m log α s ). The QCD Bethe logarithm is given in a one-parameter integral form. The rest of the formula is given as a combination of rational numbers, transcendental numbers (π, ζ(3), ζ(5)) and finite sums (besidesā 3 whose full analytic form is still unknown). We explain details of the derivation of the formula.
We have already applied the formula to an analysis of the bottomonium spectrum [27]. The current status of the perturbative prediction for the bottomonium spectrum is practically determined by a fine-level cancellation between 2m pole and V QCD (r) in E stat (r) and depends sensitively on the (yet unknown) precise value of the order α 4 s correction to the pole-MS mass relation (d 3 ). If d 3 is tuned to improve convergence of E stat (r) maximally, we observe a reasonable agreement between the predictions and experimental data of the bottomonium spectrum within estimated perturbative uncertainties. On the other hand, the prediction becomes unstable quickly if d 3 deviates from the fine-tuned value. Hence, the status of the perturbative QCD prediction at NNNLO is not yet clear until d 3 is computed.
This paper is organized as follows. In Sec. 2 we explain a method to convert a class of infinite sums to combinations of transcendental numbers and finite sums. In Sec. 3 we explain briefly the theoretical framework (pNRQCD) used in our computation of the energy levels. The computation of NNNLO corrections from the potential region is given in Sec. 4. The computation of NNNLO corrections involving the ultra-soft region is presented in Sec. 5. The obtained general formula is presented in Sec. 6. Sec. 7 is devoted to the summary and discussion. Technical details are collected in the Appendices. We set up our conventions and notations in App. A. Generating functions used for Fourier transformations are given in App. B. Formulas useful for evaluating matrix elements in terms of Coulomb wave functions are given in App. C. Evaluation of an expectation value of a spin-tensor operator is explained in App. D. An arithmetic method for deriving a local Hamiltonian in the ultra-soft correction is given in App. E. Methods for taking an angular average of matrix elements for the ultra-soft corrections are given in App. F.

Finite Sum Formula at NNLO
The known formula [15] for the NNLO correction to the energy level includes an infinite sum: For a given n, ℓ it is easily converted to a combination of transcendental numbers and rational numbers. As far as we know, however, a corresponding expression for arbitrary n, ℓ has not been known. We explain a method to convert the above infinite sum to a combination of transcendental numbers [ζ(2) = π 2 /6 and ζ(3)], rational numbers and a finite sum.
By partial fractioning, one may write where Using the summation program Wa [28], which utilizes the summation algorithm developed in [29], one may obtain a relation where S 1 (i) = i k=1 1 k denotes the harmonic sum. Furthermore, for a non-negative integer a, one may use the relation Combining Eqs. (1), (2), (4) and (5), one obtains an expression which includes only transcendental numbers and finite sums. This type of technique to convert an infinite sum to a combination of finite sums and transcendental numbers is useful in simplifying the formula for the spectrum at NNNLO.

Theoretical Framework
In this section we explain the theoretical framework which we use to compute the heavy quarkonium energy levels.
3.1 Potential-NRQCD effective field theory pNRQCD [1] is an EFT, which describes the interactions among QQ bound states and IR gluons and light quarks. Color-singlet and octet QQ composite states are represented by the fields which are local in time and bilocal in space coordinates: where the spatial coordinates (2-component spin indices) of Q andQ are denoted by x ± 1 2 r (σ andσ). In addition, IR gluons and quarks are included as dynamical degrees of freedom. More precisely, pNRQCD is obtained after integrating out the hard (h) and soft (s) modes of the dynamical fields in full QCD [5]: The dynamical fields in pNRQCD have energy-momentum in the potential (p) and ultrasoft (us) regions: The details of the EFT derivation by the mode integration can be found in Ref. [30]. The Lagrangian of pNRQCD is given as expansions in r and 1/m: H S and H O denote the quantum mechanical Hamiltonians for the singlet and octet states, respectively, which are expressed by x, r, ∂ x , ∂ r and spin operators of Q andQ. They represent the contributions from the potential region, given as instantaneous interactions (i.e., potentials), between Q andQ, and dictate the main binding dynamics of the QQ states.
The ultra-soft gluon fields are expanded in r (multipole expansion) and expressed in terms of A a µ ( x, t) and its derivatives. The leading interaction of the singlet field S and the ultra-soft gluon is given by a dipole interaction [the last terms of Eq. (12)], where the color-electric field is given by The energy levels of the heavy quarkonium states are given as the positions of poles of the full propagator of the singlet field S. The full propagator, in multipole expansion in r, is given by the diagrams in Fig. 1. In momentum space and in the center-of-mass (c.m.) frame, 1 it is given by where U O (t) ab denotes the octet propagator defined by The first term of Eq. (13) includes only contributions from the potential region, whereas the second term includes contributions from both potential and ultra-soft regions.
Since | r| is of the order of the Bohr radius, r ∼ 1/(α s m), and the energy-momentum of an ultra-soft gluon is k 0 , k ∼ α 2 s m, the multipole expansion (such as r · ∂A a 0 , r · ∂ t A a ) generates expansion in α s . The second diagram of Fig. 1 is counted as NNNLO as there are two dipole interactions and the associated couplings g 2 ∼ α s .
We decompose the Hamiltonian for the singlet state (in its c.m. frame) as where V C (r) and V NC (ˆ p, r), respectively, denote the Coulomb and non-Coulomb potentials.
Here and hereafter, m represents the pole mass of Q orQ. Below we present each part of the potentials, which contributes to the NNNLO corrections to the quarkonium energy levels. These are determined in perturbative QCD by integrating out the hard and soft modes. For convenience, we use both coordinate-space and momentum-space representations of the potentials. In coordinate space we write each term of the potentials in the form whereˆ p = −i ∂ r is a derivative operator. The corresponding potential in momentum space is represented by with p ′ = p + q, and We employ dimensional regularization with one temporal dimension and d = D−1 = 3−2ǫ spatial dimensions. γ E = 0.5772 · · · denotes the Euler constant. The relation between both representations is given explicitly by

Coulomb potential
In purely perturbative expansion in α s , the Coulomb potential V C (r) is identical to the static QCD potential, which is defined from the expectation value of the Wilson loop and subtracting the ultra-soft contribution. At order α 4 s and beyond, V C (r) is IR divergent. Generally we separate V C (r) to the renormalized part V C,ren (r) and IR divergent part V C,div (r). In the computation of ultra-soft corrections, V C,div (r) plays the role of a counter term [1] .
We define the separation of 1/ǫ-pole and the renormalized part in momentum space to conform with the conventional MS scheme. This is also so for the ultra-soft part for consistency [see Eqs. (116) and (117)]. In momentum space, the Coulomb potential up to order α 4 s is given by The V -scheme coupling constant is given as a perturbative expansion in the MS coupling constant α s ≡ α s (µ) as where the coefficients of logarithms in P V n are determined by the renormalization group as The C 3 A term in P V 3 is a logarithm generated from a combination of the 1/ǫ pole of V C,div and (µ 2ǫ ) 3 of the three-loop integration measure in dimensional regularization. The parameters a 0 , a 1 , a 2 ,ā 3 and β 0 , β 1 , β 2 are given in App. A. The renormalized Coulomb potential in coordinate space is given by where L r = log(e 2γ E µ 2 r 2 ). For the renormalized part, Fourier transformation formula in d = 3 can be applied.

Non-Coulomb potentials
The non-Coulomb potentials also contain IR divergences, and they are separated to the renormalized and divergent parts in the same manner as the Coulomb potential. The NNNLO non-Coulomb potentials were obtained in Ref. [20]. In momentum space the renormalized non-Coulomb potentials are given as follows. where We parameterize the potential coefficients as follows: The 1/m potential coefficients at tree, one-loop and two-loop levels are given by The 1/m 2 potential coefficients at tree level are given by and at one-loop level by where L m = log(µ 2 /m 2 ) and L q = log(µ 2 /q 2 ). There are contributions from annihilation diagrams, which should be added to d δ i and d s i . The absorptive part is discarded in the computation of the energy levels.
In coordinate space the renormalized non-Coulomb potentials are given as follows.
where δ ℓ≥1 = 1 if ℓ ≥ 1 and zero otherwise. The Fourier transforms are computed using the formula in App. B, Eq. (192) and its derivatives with respect to r. "reg[1/r 3 ]" in V δ can be identified with 1/r 3 when applying to ℓ > 0 states but needs some modification when computing matrix elements for the ℓ = 0 states. We will come back to this point when explaining the relevant computation in Sec. 4.4.

IR divergent potentials
The IR divergent parts of the Coulomb and non-Coulomb potentials can be combined and are given in momentum space as This corresponds to Eq. (4.59) of Ref. [30] (see also Ref. [31]), in which δV c.t. is defined as a counter term and has the opposite sign compared to V pot div .

NNNLO Corrections: Potential Region
In this section we compute contributions to the NNNLO corrections to the energy levels which originate from the potential region. This corresponds to the pole positions of the singlet propagator as determined by the first diagram of Fig. 1 or by the first term of Eq. (13).

Potential perturbation
The leading-order Hamiltonian of the singlet bound state is taken as that of the pure Coulomb system: All the other parts of H S in Eq. (15) are treated as perturbations. We compute the perturbative expansion of the singlet Green function, which can be expressed as where we denote by V i the i-th order potential 2 . The energy is measured from the threshold, E = √ s − 2m. G (0) is the leading-order Green function given by In the following we suppress r, r ′ and the infinitesimal imaginary part i0 of the Green function denominators. Perturbative corrections to the Green function are given by potential insertions: The Green function has single poles corresponding to bound states in the complex energy plane. Suppose we are interested in the bound state, which is specified at leading-order by H 0 | n = E (0) n | n . (In this subsection n represents a set of quantum numbers specifying a bound state, for simplicity of notations.) The pole of the Green function for the corresponding bound state can be written as where represent the perturbative expansions of the energy eigenvalue and residue (wave function) of the bound state, respectively. Hence, the perturbative expansion of Eq. (59) is given by By comparison with the potential perturbation, we obtain a master formula for the energy corrections up to NNNLO. It is convenient to define the reduced Green function and its first derivative with respect to E as follows.
The energy and wave function corrections can be extracted from the double and single poles, respectively. The first-order corrections can be extracted from δ 1 G(E) and we obtain The second-order corrections to the energy and wave function can be extracted from the pole structure of δ 2 G(E) and the results of the previous order. We obtain The energy correction at third order is given by Thus, for the third order corrections to the energy levels, we need to compute a single insertion n|V 3 |n , double insertions n|V 1 G n V 2 |n , n|V 1 ∂G n V 1 |n and a triple insertion n|V 1 G n V 1 G n V 1 |n . The details of methods to compute all the insertions are explained in the following subsections.

Wave functions and Green function of Coulomb system
To apply the above master formula we need explicit representations for the wave functions and the Green function of the leading-order Coulomb system. These are given as follows.
Let us define the relevant variables as The radial part of the wave function for the color-singlet bound state is give by which is normalized as . The octet wave function is given by which is normalized as We use the following representation of the Coulomb Green function for our computation: ℓ denotes the Legendre polynomial, andr i = r i /r i represents the unit vector in the direction of r i . The Green function for the partial wave ℓ is given by an infinite sum 3 where R νℓ is the singlet bound-state wave function. The reduced Green function G nℓ for the partial wave ℓ can be computed as follows. Let Then with The first derivative of the radial wave function can be rewritten in terms of the same radial wave function and the one with the index β(= n − ℓ − 1) lowered by one:

Coulomb corrections
We compute the corrections which originate purely from the renormalized Coulomb potential Eq. (27) in the potential perturbation of Sec. 4.1. Using the reduced Green function in Sec. 4.2 and formulas in App. C, it is cumbersome but straightforward to obtain the Coulomb corrections expressed as combinations of infinite sums. Then using techniques similar to that of Sec. 2, we can express the corrections with finite sums. More explicitly, we proceed as follows. We compute E Integrals over the radial variables can be performed using the formulas Eqs.

Non-Coulomb corrections
We compute the NNNLO corrections which involve the non-Coulomb potentials. Since there are no NLO potentials in the non-Coulomb potentials, there are only double insertions and single insertion in the potential perturbation.

Double insertion: spin-dependent part
We first consider the spin-dependent part of the double insertion of potentials. This contribution is given by where denotes the NLO Coulomb potential, and represents the spin-dependent part of the NNLO non-Coulomb potentials.
In computing E NC 2,a we can replace U 2 by where The expectation value · · · is taken with respect to the (ℓ, s, j) state. The derivation of D S is given in App. D. The above replacement is justified, since V 1 and G nℓ do not change | ℓsj , hence only matrix elements proportional to ℓsj | U 2 | ℓsj = U ′ 2 appear. One can compute E NC 2,a easily similarly to the Coulomb corrections. One needs to compute the corrections for the ℓ = 0 and ℓ > 0 cases separately, due to existence of δ( r) and 1/r 3 . In particular, in the case ℓ = 0 one should set X LS = D S = 0 before evaluating the matrix elements of 1/r 3 . (See explanation for the single insertion below.)

Double insertion: spin-independent part
Next we consider the spin-independent part of the double insertion of potentials. This contribution is given by where represents the spin-independent part of the NNLO non-Coulomb potentials. E NC 2,b can be computed easily by eliminating theˆ p operator in the following way.
with V 0 = −C F α s /r and Matrix elements which are needed to evaluate the double insertions of the NLO Coulomb and U 1 are given as

Single insertion
Finally we compute the single insertion of the NNNLO non-Coulomb potentials in potential perturbation: We use the renormalized non-Coulomb potentials in coordinate space given in Eq. (44). We take only the NNNLO part as V 3,NC , namely those potentials which are not included in U 1 + U 2 . In computing the above matrix element, one needs to be careful in evaluating the contribution of V δ . The operators δ( r) and reg[1/r 3 ] are singular when their matrix elements are computed with respect to ℓ = 0 states. This is an artifact of working in coordinate space, since no singularities arise if the same matrix elements are evaluated in momentum space. Nevertheless, it is advantageous to work in coordinate space for obtaining a general formula for arbitrary (n, ℓ, s, j). Therefore, we compute in the following way. δ( r) and reg[1/r 3 ], respectively, stem from the Fourier transforms of 1 and L q = log(µ 2 /q 2 ). There are no singularities if we define them as the u → 0 limits of (µ 2 /q 2 ) u and ∂ u [(µ 2 /q 2 ) u ], respectively, in momentum space. Hence, we first take the Fourier transform of (µ 2 /q 2 ) u and ∂ u [(µ 2 /q 2 ) u ] using the formula Eq. (192), then evaluate the matrix elements in coordinate space using the formulas in App. C. We take the limit u → 0 in the end and obtain the matrix elements. In this way we confirm that the operator δ( r) can be handled in the standard way. On the other hand, one finds that the matrix elements of the operator reg[1/r 3 ] should be evaluated as n, ℓ = 0 reg 1 r 3 n, ℓ = 0 = 1 4 for the ℓ = 0 states, whereas reg[1/r 3 ] can be identified with 1/r 3 for the matrix elements with respect to ℓ > 0 states.

IR divergent contributions
The IR divergent contributions which appear first at NNNLO is given by the single insertion of the divergent potentials Eq. (52): As we will see in the next section, this entire contribution is canceled by the UV divergent part of the corrections nℓsj|H us div |nℓsj from the ultra-soft region. Therefore, we do not evaluate explicitly the matrix element for the divergent contributions.

NNNLO Corrections: Ultra-soft Region
In this section we compute contributions to the NNNLO corrections to the energy levels which originate from both ultra-soft and potential regions. This corresponds to the selfenergy of the singlet propagator in the second diagram of Fig. 1 or in the second term of Eq. (13).

Separation of local and non-local parts
The one-loop self-energy of the singlet field S is given by (13), which corresponds to replacing D t in Eq.
In Eq. (105), · · · nℓ denotes the expectation value taken with respect to the energy eigenstate of H (d) S specified by the quantum numbers (n, ℓ, s, j, j z ) (the expectation value depends only on n, ℓ). Only the leading-order terms of H S and H O are needed in our computation, whereas we need to keep ǫ = (3 − d)/2 non-zero until we extract the UV divergence of E us nℓ explicitly. The gluon emitted and reabsorbed by QQ after the time interval t is in the ultra-soft region. On the other hand, the propagator of the octet field e −iH O t contains multiple exchanges of Coulomb gluons which are in the potential region. Therefore the whole correction consists of a combination of contributions from the ultra-soft and potential regions.
We separate E us nℓ into two parts: the part given by an expectation value of a Hamiltonian (local in time) and the part given by a term known as the "QCD Bethe logarithm" (non-local in time). This is achieved in the following way. The correlation function of the color electric field (non-local in time and local in space) is evaluated in the lowest order in expansion in α s as After integrating over t and using we obtain where we have replaced E (d) n,C by the singlet Hamiltonian inside the expectation value, taking into account ordering of the operators. We rewrite X by moving all the momentum operatorsˆ p = −i ∂ r to the right of the coordinate operators r and r = | r| using commutation relations, and then we Fourier transform the result. After expanding in ǫ and using the equation of motion, we obtain (see App. E for details) After adding the contributions from the counter terms, the ultra-soft correction to the bound state energy becomes finite and is given by We have split the logarithmic term into the µ dependent part and µ independent part, where the QCD Bethe logarithm is defined as All the other parts can be evaluated analytically. In particular, the expectation value in the second term of Eq. (119) can be written as which is related to the UV divergence.

One-parameter integral form of QCD Bethe logarithm
The QCD Bethe logarithm for the S-wave states was first computed in Ref. [23]. We derive the one-parameter integral form for the Bethe logarithm for a general quantum number in the following way. We insert the completeness relation for the octet states: Decomposing the radial and angular parts of the wave function we obtain (see App. F for derivation) Using the wave functions in Sec. 4.2, the radial part can be written as The matrix element X ℓ ′ nℓ is defined as Changing the variable to x = 1/k, we find where the summation is taken over ℓ ′ = ℓ ± 1 if ℓ ≥ 1 and ℓ ′ = 1 if ℓ = 0. Here, and Each X ℓ ′ nℓ is a rational function of x. We list explicit forms for some of the lower states: Thus, for each (n, ℓ), the integrand of Eq. (127) is given by a combination of elementary functions, and the integral can be evaluated easily numerically. For instance, Numerical values of the Bethe logarithm are given in Tab. 1 for some of the lower states.

Full Formula
Combining all the corrections we present the full formula of the heavy quarkonium energy levels up to NNNLO, 5 that is, up to O(α 5 s m) and O(α 5 s m log α s ). The heavy quarkonium state is identified by the quantum numbers n, ℓ, s and j. In this section we write the pole mass of Q orQ as M Q,pole instead of m, to avoid confusions, since m is used to represent the summation index corresponding to ℓ z .  Table 1: Numerical values of the Bethe logarithm L Bethe (n, ℓ) for some of the lower states. The color factors are set as C F = 4/3 and C A = 3.
The energy level is given by where the coefficients of logarithms of µ in P i are determined by the renormalization group as The color factors, C A , C F , T F , the coefficients of the beta function, β 0 , β 1 , β 2 , as well as the constants of the static QCD potential, a 1 , a 2 ,ā 3 , to be used below are given in App. A. Hereafter, we use δ ℓ0 and δ ℓ≥1 (equal to one if ℓ ≥ 1 and zero if ℓ = 0) to separate the cases ℓ = 0 and ℓ ≥ 1. We separate the constant part as c i = c c i + c nc i , where c c i corresponds to the contributions purely from the (renormalized) Coulomb potential, while c nc i corresponds to the sum of all the rest of the contributions: + Σ τ,1 (n, ℓ) + Σ τ,2 (n, ℓ) + Σ τ,3 (n, ℓ).
Non-Coulomb corrections are classified according to different color factors: The Bethe logarithm is given as a one-parameter integral: where X ℓ ′ nℓ and Y ℓ ′ nℓ are defined as 384 with The hypergeometric function is defined as 2 F 1 (a, b; c; z) = Γ(c) Γ(a+n)Γ(b+n) Γ(c+n) z n n! . Each X ℓ ′ nℓ (x; ρ n ) is a rational function of x. Explicit expressions of X ℓ ′ nℓ (x; ρ n ) and numerical values of the Bethe logarithm L Bethe (n, ℓ) are given in Sec. 5.2, for some of the lower states.
The above formula reproduces the known results: (1) the coefficients of α 5 s m log α s for general quantum numbers (n, ℓ, s, j) [24], and (2) the NNNLO formula for the general S-wave state (n, j) [25].

Summary and Discussion
We have computed the energy levels of the heavy quarkonium states in a double expansion in α s and log α s up to order α 5 s m and α 5 s m log α s . The result is given as a general formula dependent on the quantum numbers (n, ℓ, s, j) of the bound state. The computation is performed using the theoretical framework pNRQCD.
In pNRQCD the perturbative corrections consist of contributions from the potential and ultra-soft regions. Decomposition of both contributions is realized by a multipole expansion. The corrections from the potential region (potential corrections) are identical to perturbative corrections of energy levels in quantum mechanics. The relevant Hamiltonian has been known, and it is straightforward to obtain an infinite sum formula using a known infinite-sum representation of the Green function. We used a recent technology for evaluating multiple sums to reduce infinite sums to combinations of transcendental numbers and finite sums.
Corrections involving the ultra-soft region (ultra-soft corrections) start from the order α 5 s m and α 5 s m log α s . The relevant correction in our computation is given by a one-loop self-energy of a singlet QQ state by emitting and reabsorbing an ultra-soft gluon. It can be separated to a part given by an expectation value of a Hamiltonian (local in time) and a part given by a contribution non-local in time (QCD Bethe logarithm). A UV divergence is included in the Hamiltonian, which cancels the IR divergence contained in the potential corrections. A new feature of our computation is that these are computed in an arithmetic manner, in contrast to previous computations which involve diagrammatic analysis. Hence, our method may be easier to follow for non-experts. The Bethe logarithm is given as a one-parameter integral. Its numerical values for some lower states are given. All the other parts are evaluated analytically, besidesā 3 which is as yet known only numerically.
The obtained formula is quite lengthy. One reason is that part of the formula needs to be given separately for the cases ℓ = 0 and ℓ > 0. Another reason is that there appear various different types of finite sums. (They may be reduced to more compact expressions if we understand the nature of the finite sums better.) In general, there are two ways to compute the energy levels using perturbative QCD. One way is to compute thoroughly within perturbative QCD. (Technically this is done efficiently using an EFT.) The other way is to compute by factorizing UV and IR contributions in an operator-product expansion. In the former computation, there are wellestablished prescriptions to estimate uncertainties of the prediction within perturbative QCD. Empirically estimates of perturbative uncertainties are approximated well by IR renormalons of order Λ 3 QCD a 2 X , where a X denotes the typical radius of the bound state X. In the latter computation, UV contributions are encoded in the Wilson coefficients, which are free from IR renormalons, 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 ∼ Λ 3 QCD a 2 X 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. See discussion in [27] for more details in the case of the bottomonium spectrum.

Appendices A Conventions and Notations
In this appendix we present definitions of parameters and conventions used in the main body of this paper.
The color factors for the SU(N C ) gauge group are given by The coefficients of the beta function are given by where n l denotes the number of massless quark flavors.
The strong coupling constant in the MS scheme is defined in the following way. The bare gauge coupling constant is expressed as We choose the counter terms to be only multiple poles in ǫ in dimensional regularization: where each z n (g R ) is given as a series expansion in g R . Then the strong coupling constant in the MS scheme is given by

B Formulas for Fourier Transformation
We can use the following formulas for Fourier transforms in d = 3 − 2ǫ spatial dimensions: where q = | q| and r = | r|. If we differentiate F and G by r and q, respectively, we may also obtain formulas for Fourier transforms of tensor operators.

C.1 Formulas
The generating function of the Laguerre polynomial is given by It is easy to evaluate an integral of the generating functions ∞ 0 z a+β e −z U a (z, t)U a (z, s) = Γ(1 + a + β) Expansions in s and t read

Diagonal part
If we take the diagonal part of Eq. (195), we obtain By taking the limits β → −2, −1, 0, · · · , the following relations are obtained. z a+β e −z L a n (z)L a m (z) = ∞ 0 z a+β e −z U a (z, t)U a (z, s) The off-diagonal (t n s m ) component corresponds to i + k = n, j + k = m, with constraints k ≤ n ∧ k ≤ m ∧ n = m. This gives We derive some formulas assuming n > m.
In a similar manner, by evaluating ∂ 2 r [(μr) 2(u+ǫ) /r], one can show (μr) 2(u+ǫ) r 3 i r ·ˆ p = iu r 3 (μr) 2(u+ǫ) We use this relation to eliminate r ·ˆ p from X. The first term of Eq. (225) contributes only to the finite part of H us , hence we may use the relation which holds for d = 3, to eliminate 1 r 3 r i r jpjpi . Note that the left-hand-side vanishes inside the expectation value.
We take the Fourier transform of X using the formulas in App. B. We insert the result in Eq. (111) and expand the local operators in ǫ and take the limit u → 0. Furthermore, we add and obtain eqs. (114)-(117).

F Angular Average
We explain how to take the sum over ℓ ′ z in Eq. (123). We can do this in two different ways.
We may replace the external brackets x |, | y of T x,y by ℓ, m |, | ℓ, m and obtain