${\overline{\rm MS}}$ renormalization of $S$-wave quarkonium wavefunctions at the origin

We compute $S$-wave quarkonium wavefunctions at the origin in the $\overline{\rm MS}$ scheme based on nonrelativistic effective field theories. We include the effects of nonperturbative long-distance behaviors of the potentials, while we determine the short-distance behaviors of the potentials in perturbative QCD. We obtain $\overline{\rm MS}$-renormalized quarkonium wavefunctions at the origin that have the correct scale dependences that are expected from perturbative QCD, so that the scale dependences cancel in physical quantities. Based on the calculation of the wavefunctions at the origin, we make model-independent predictions of decay constants and electromagnetic decay rates of $S$-wave charmonia and bottomonia, and compare them with measurements. We find that the poor convergence of perturbative QCD corrections are substantially improved when we include corrections to the wavefunctions at the origin in the calculation of decay constants and decay rates.


Introduction
Heavy quarkonium production and decay processes are multiscale problems that are sensitive to both short-distance and long-distance natures of QCD. As many of these processes have been measured experimentally, a great amount of effort has been made towards understanding them theoretically [1,2]. Much of the heavy quarkonium phenomenology is based on nonrelativistic effective field theories, which provide factorization formalisms that separate the perturbative short-distance physics from nonperturbative long-distance physics.
In the nonrelativistic QCD (NRQCD) factorization formalism [3,4], production or decay rates of a heavy quarkonium are given by sums of products of the perturbatively calculable short-distance coefficients (SDCs) and long-distance matrix elements (LDMEs). The LDMEs are nonperturbative quantities that correspond to the probability to find a heavy quark Q and a heavy antiquarkQ inside a quarkonium. The LDMEs have known scalings in v, the typical heavy-quark velocity inside the quarkonium, and the sum is organized in powers of v. While the SDCs can be computed in perturbative QCD, accurate determinations of the LDMEs, especially the ones that appear at the lowest orders in v, are also important in making predictions of production or decay rates of heavy quarkonia. So far, many phenomenological studies on heavy quarkonium production and decay have relied on model calculations of the LDMEs [5][6][7][8][9]. One major disadvantage of model calculations is that in general, they do not reproduce the correct ultraviolet (UV) behaviors of the LDMEs that are predicted in perturbative QCD. Perturbative QCD calculations show that the LDMEs contain ultraviolet (UV) divergences, which require renormalization [4]. That is, the LDMEs are renormalization scheme dependent. The scale associated with the renormalization of the LDMEs is often called the NRQCD factorization scale. The SDCs also depend on the scheme in which the LDMEs are renormalized, in the way that the scheme dependence cancels between the SDCs and the LDMEs in the factorization formula. It has been found from perturbative QCD calculations that strong dependencies on the factorization scale start to appear in the SDCs from two loops [10][11][12][13]. Therefore, in order to make accurate predictions based on NRQCD, it is critically important to determine the LDMEs that exhibit the correct scale dependence. Since perturbative QCD calculations are most conveniently done in dimensional regularization (DR), the SDCs are usually computed in the MS scheme. In order to be consistent with the MS calculations of the SDCs, the LDMEs must also be determined in the MS scheme.
While lattice QCD determinations of certain LDMEs exist [14][15][16][17][18][19], these calculations are usually done in quenched lattice QCD, and their results have large uncertainties. Moreover, the relations between the LDMEs in lattice and continuum are known only at one-loop level. Hence, existing lattice QCD determinations are not accurate enough to reproduce the factorization-scale dependence that is expected in perturbative QCD.
It has been known that NRQCD LDMEs can be computed from quarkonium wavefunctions at the origin [4]. Rigorous formulations for quarkonium wavefunctions have been developed in the potential NRQCD (pNRQCD) effective field theory approach [20][21][22]. This formalism provides a Schrödinger formulation, from which the quarkonium wavefunctions can be computed. The potentials that appear in the Schrödinger equation have fieldtheoretical definitions in terms of Wilson loops, and they can be computed nonperturbatively in lattice QCD [23,24]. While this makes possible the nonperturbative determination of quarkonium wavefunctions, there are still challenges in computing the wavefunctions at the origin from first principles. One major challenge is that the wavefunctions at the origin involve divergences that require renormalization. These divergences are related closely to the UV divergences that appear in the LDMEs. In order to obtain the MS-renormalized LDMEs, the wavefunctions at the origin must also be renormalized in the same scheme. The problem is that this requires dimensionally regulated calculations, which are difficult to be done outside of perturbation theory. This is because calculations in DR are most conveniently done in momentum space, while nonperturbative determinations of the potentials are done in position space. For this reason, computations of quarkonium wavefunctions at the origin to two-loop accuracy have only been done within perturbative QCD [25][26][27][28][29][30][31][32][33][34][35], where the nonperturbative, long-distance behavior of the potentials are ignored. However, many charmonium and bottomonium states are non Coulombic, so that their wavefunctions are sensitive to the nonperturbative long-distance behavior of the potentials. In such cases, the nonperturbative behavior of the potentials cannot be neglected.
While a direct momentum-space calculation of the wavefunctions at the origin in DR with nonperturbative potentials may be difficult, position-space calculations are possible if we regulate the divergences in position space. Then, we only need to convert the positionspace regularization to the MS scheme in order to obtain the MS-renormalized wavefunctions at the origin. The conversion from position-space regularization to the MS scheme may be computed in perturbative QCD, because this depends only on the divergent shortdistance behavior of the potentials, which are determined completely by perturbative QCD. This makes possible the nonperturbative calculations of MS-renormalized wavefunctions at the origin based on first principles.
In this paper, we compute the MS-renormalized quarkonium wavefunctions at the origin for S-wave charmonium and bottomonium states. We compute the wavefunctions at the origin in two steps. First, we compute the quarkonium wavefunctions in position space, by using potentials that have nonperturbative long-distance behaviors that are determined in lattice QCD, while the potentials are given by perturbative QCD at short distances. We use quantum-mechanical perturbation theory to first order in the expansion in powers of 1/m to compute the wavefunctions at the origin. Then, we convert the position-space regularization to the MS scheme. Using the MS-renormalized quarkonium wavefunctions at the origin that we obtain, we determine the NRQCD LDMEs in the strongly coupled pNRQCD formalism, which is valid for non Coulombic quarkonia. Based on the determinations of the LDMEs, we make model-independent predictions of decay constants and electromagnetic decay rates of S-wave charmonium and bottomonium states. In the NRQCD factorization formulas, we include loop corrections at leading order in v to two-loop accuracy, as well as corrections of order α 0 s v 2 , and, when available, we also include corrections of order α s v 2 . We restrict the calculation of wavefunctions to S-wave states, because the two-loop corrections to the SDCs are generally not available for the production or decay rates of quarkonia with higher orbital angular momentum. This paper is organized as follows. In section 2, we review the definitions of NRQCD LDMEs and the relations with wavefunctions at the origin. We outline the calculation of quarkonium wavefunctions in position space in sec. 3, which allow nonperturbative calculations of the wavefunctions at the origin with a position-space regulator. We compute the scheme conversion from position-space regularization to the MS scheme in sec. 4. In sec. 5, we compute the MS-renormalized wavefunctions at the origin, as well as electromagnetic decay rates and decay constants of S-wave charmonium and bottomonium states, which we compare with experimental measurements and lattice QCD determinations. We conclude in sec. 6.

NRQCD long-distance matrix elements
In this section, we review the definitions of NRQCD LDMEs involving S-wave quarkonia and their relations to quarkonium wavefunctions that appear in pNRQCD.
In NRQCD factorization formulas for electromagnetic decay rates and exclusive electromagnetic production rates of vector quarkonium V = J/ψ or Υ, the following LDME appears at leading order in v: where ψ and χ are Pauli spinor fields that annihilate and create a heavy quark and a heavy antiquark, respectively, and is the polarization vector of the quarkonium. We take nonrelativistic normalization for the state |V . At relative order v 2 , the following LDME appears: where D = ∇ − ig s A, and χ † ← → D ψ = χ † Dψ − (Dχ) † ψ. The leading-order LDME depends on the factorization scale Λ from its renormalization. This factorization scale dependence is given by the following evolution equation [4,10,11,13] 3) where m is the heavy quark pole mass, C F = (N 2 c − 1)/(2N c ), C A = N c , and N c = 3 is the number of colors. We reproduce the anomalous dimensions on the right-hand side of eq. (2.3) from NRQCD loop calculations in appendix A.
In strongly coupled pNRQCD, which is valid in the regime mv Λ QCD mv 2 , the leading-order LDME is given by [36,37] 0|χ † · σψ|V 4) while the order-v 2 LDME is given by Here, Ψ V (r) is the quarkonium wavefunction for the V state, E V is the binding energy, and c F is the short-distance coefficient associated with the spin-dependent operator in the NRQCD Lagrangian [38]. The Ψ V (r) and E V are eigenfunctions and eigenenergies of the Schrödinger equation, which we introduce in the next section. The E 1 , E (2,em) 3 , and B 1 are nonperturbative gluonic correlators of mass dimension two, which scale like Λ 2 QCD . The quantity E 3 is a dimensionless gluonic correlator. The field-theoretical definitions of the gluonic correlators can be found in refs. [36,39]. Equation (2.4) implies that the factorization scale dependence must come from E 3 and |Ψ V (0)|, because the gluonic correlators of mass dimension two are scaleless power divergent in perturbative QCD, and E V is finite. The order-α s v 2 contribution to the anomalous dimension in eq. (2.3) is consistent with the known scale dependence of E 3 [37,39]. As a result, the two-loop anomalous dimension in the first term on the right-hand side of eq. (2.3) must come from the scale dependence of the wavefunction at the origin |Ψ V (0)|.
Analogously, in electromagnetic decay rates and exclusive electromagnetic production rates of pseudoscalar quarkonium P = η c or η b , the following LDME appears in factorization formulas at leading order in v: 0|χ † ψ|P , (2.6) and at relative order v 2 , the LDME 0|χ † (− i 2 ← → D ) 2 ψ|P appears. The factorization scale dependence of the leading-order LDME is given by [4,12,13] We reproduce the anomalous dimensions on the right-hand side of eq. (2.7) from NRQCD loop calculations in appendix A. The leading-order LDME is given in strongly coupled pNRQCD by [36,37] and the order-v 2 LDME is given by where Ψ P (r) is the quarkonium wavefunction for the P state, and E P is the corresponding binding energy. Similarly to the case of vector quarkonia, the known scale dependence of E 3 coincides with the second term on the right-hand side of eq. (2.7), and so, the two-loop anomalous dimension in the first term on the right-hand side of eq. (2.7) must come from the scale dependence of the wavefunction at the origin |Ψ P (0)|. In order to obtain the wavefunctions at the origin |Ψ V (0)| and |Ψ P (0)| in the MS scheme with the correct dependence on the scale, it is necessary to compute the quarkonium wavefunctions from the Schrödinger equation with the potential that has the correct short-distance behavior that is expected from perturbative QCD. For many charmonium and bottomonium states, it is also necessary to include the nonperturbative long-distance behavior of the potential that is not captured in perturbative QCD calculations. As we have discussed previously, in order to include the long-distance behavior of the potential, it is most convenient to work in position space, where the divergences are regulated by a position-space regulator. In the following sections, we discuss our strategy to compute the wavefunctions at the origin with a position-space regulator, and compute the conversion from the position-space regularization to the MS scheme.

S-wave quarkonium wavefunctions in position space
In pNRQCD, the quarkonium wavefunctions Ψ n (r) are determined from a Schrödinger equation, which can be written in the form where r is the distance between the Q and theQ, E n is the binding energy for the state n, and V (r, ∇) is the potential, which is determined order by order in the expansion in powers of 1/m. Since we are only interested in S-wave states, the wavefunctions depend only on r = |r|. At leading order in 1/m, the potential V (r, ∇) is given by the static potential V (0) (r), which has a nonperturbative definition in terms of a Wilson loop [21,[40][41][42]. For r 1/Λ QCD , the static potential is determined by perturbative QCD, which gives V (0) (r) = −α s C F /r at leading order in α s . The potential including the correction terms of order 1/m and 1/m 2 can be written generically as 2) where we include only the contributions relevant for S-wave states up to order 1/m 2 , and the ellipsis represent terms of higher orders in 1/m. Here, S is the QQ spin, S 2 = 2 for the spin-triplet state, and S 2 = 0 for the spin-singlet state. Because unitary transformations can reshuffle the 1/m terms with the 1/m 2 terms in the potential, if we want to include corrections of order 1/m, we must also include the correction terms of order 1/m 2 in the potential. Even though the last term in eq. (3.2), which originates from the relativistic correction to the kinetic energy, is suppressed by 1/m 3 , this term must be regarded as a 1/m contribution, because a power of −∇ 2 /m can be traded with a power of E n − V (0) (r) by using the Schrödinger equation. In the same way, higher order corrections to the kinetic energy of the form ∇ 2n /m 2n−1 for n ≥ 3 are suppressed by at least 1/m 2 . If we assume that E n and V (0) (r) are of order mv 2 , these higher order corrections are suppressed by higher powers of v compared to the order 1/m and 1/m 2 terms included in eq. (3.2). Hence, we neglect the higher order corrections to the kinetic energy of the form ∇ 2n /m 2n−1 for n ≥ 3.
The form of the potential beyond the static one depends on the scheme in which the matching between NRQCD and pNRQCD is done. Nonperturbative definitions of the potentials that appear in eq. (3.2) are found from Wilson-loop matching [23,24]. On the other hand, it is necessary to employ the results from on-shell matching [43,44] in order to obtain dimensionally regulated wavefunctions at the origin that is consistent with the SDCs, because dimensionally regulated calculations of the SDCs are also done by matching on-shell amplitudes in QCD and NRQCD. We note that the potentials are gauge invariant in both cases. The different forms of the potentials can be related by unitary transformations [23,36,45].
The order-1/m and 1/m 2 potentials in eq. (3.2) are to be considered perturbations in the quantum-mechanical perturbation theory (QMPT), where the wavefunctions and binding energies are first computed at leading order from the Schrödinger equation without including the corrections to the potential that are suppressed by powers of 1/m. Then, the corrections of higher orders in 1/m are included by using the Rayleigh-Schrödinger perturbation theory.
If we ignore the terms suppressed by powers of 1/m and keep only the static potential in eq. (3.2), then S-wave wavefunctions are finite at the origin r = 0; this is because the behavior of the wavefunctions near r = 0 is determined by the short-distance behavior of the static potential, which diverges like 1/r at r = 0. Therefore, the corrections to the wavefunctions at the origin are finite if the corrections come from potentials that diverge at most like 1/r. On the other hand, the 1/m and 1/m 2 terms in eq. (3.2) can produce divergences in the wavefunctions at the origin if they diverge faster than the static potential at r = 0. We list the short-distance behavior of the potentials from on-shell matching and Wilson-loop matching in appendix B. The divergent behavior of the wavefunctions at r = 0 can be inferred from nonrelativistic quantum mechanics. Since the 1/m potential V (1) (r) diverges like 1/r 2 at r = 0, the first order correction to the wavefunctions from the 1/m potential produces a logarithmic divergence that is proportional to α 2 s log r at r = 0. As we will see later, the velocity-dependent potential V (2) p 2 (r) and the relativistic correction to the kinetic energy −∇ 4 /(4m 3 ) also produce logarithmic divergences that are similar to the correction from the 1/m potential. The 1/m 2 potential includes the delta function δ (3) (r), and this produces at first order in the QMPT a power divergence proportional to 1/(mr) to the wavefunctions at r = 0.
We compute the wavefunctions and the corrections of higher orders in 1/m in the following way. We define the leading-order (LO) potential V LO (r) from the static potential V (0) (r) by subtracting the perturbative corrections of order α 2 s and beyond, but keeping the long-distance nonperturbative behavior. The specific form of the LO potential that we use will be given in sec. 5. This makes V LO (r) behave like −α s C F /r at short distances, while it coincides with the static potential at long distances. The perturbative corrections of higher orders in α s will be included as perturbations in the QMPT. The LO wavefunction Ψ LO n (r) and the binding energy E LO n satisfy the LO Schrödinger equation To first order in the QMPT, the wavefunction Ψ(r) for an S-wave state n is given in terms of the LO wavefunction Ψ LO n (r) by where δV (r, ∇) = V (r, ∇) − V LO (r).Ĝ n (r , r) is the reduced Green's function for the eigenstate n, which is defined bŷ where the sum runs over all eigenstates of the LO Schrödinger equation except for the state n. Although the sum includes states with nonzero orbital angular momentum, only S-wave states contribute to the integral in eq. (3.5) due to the rotational symmetry of δV (r, ∇). The reduced Green's function is related to the Green's function G(r , r; E) bŷ while G(r , r; E) is defined for arbitrary complex E by The Green's function satisfies the equation The Green's function in position space can be computed by using the formal definition in eq. (3.8), or by solving the differential equation in eq. (3.9). Note that the reduced Green's function can be computed from G(r , r; E) by usinĝ We note that the reduced Green's function satisfies d 3 rĜ n (r , r)Ψ LO n (r) = 0, (3.12) which follows from the orthogonality of wavefunctions. The vanishing of eq. (3.12) also follows from the fact that adding or subtracting a constant to the potential V (r, ∇) have no effect on the wavefunctions Ψ n (r). The corrections to the wavefunction δΨ n (r ) can be computed from eq. (3.5). The corrections from the velocity-dependent potential and the relativistic correction to the kinetic energy contain ∇ 2 , which can be reduced by using the Schrödinger equation and eq. (3.10). The correction from the velocity-dependent potential reads It is clear that the first term in the last line is finite at r = 0, since Ψ LO n (r) is regular at r = 0, and V (2) p 2 (r) diverges like 1/r at r = 0. While in perturbative QCD the quantity V (2) p 2 (0) is a scaleless power divergence and therefore vanishes, the nonperturbative definition of the velocity-dependent potential in terms of Wilson loops implies that V (2) p 2 (0) is proportional to a gluonic correlator iE 2 , which is defined in refs. [36,37]. Similarly, the correction from the −∇ 4 /(4m 3 ) term is given by Again, it is clear that the first term in the last line is finite at r = 0. The LO potential at zero distance V LO (0) vanishes in perturbative QCD, because it is a scaleless power divergence. This quantity also vanishes nonperturbatively, which follows from the exact vanishing of the static potential at zero distance [36]. By using the expressions in eqs. (3.13) and (3.14), we can write the first order correction to the wavefunction as where δV C (r) = V (0) (r) − V LO (r), and We note that for an S-wave state n, each term in eq. (3.15) is a function of r = |r |, and does not depend on the angles of r . Because of the explicit rotational symmetry of δV (r, ∇) and δV(r), eq. (3.15) is unchanged if we only include S-wave states in the definition of the reduced Green's function in eq. (3.6).
In eq. (3.15), the divergences at r = 0 are contained in the first integral. The second integral in eq. (3.15) is finite at r = 0, because the terms in the curly brackets diverge at most like 1/r at r = 0. The UV divergence in the first integral can be cut off by setting r = r 0 with r 0 > 0, instead of setting r = 0. This defines the finite-r regularization [46], which is the position-space regularization that we use in this paper. We note that a similar version of position-space regularization has been used in refs. [26,28,29,31]. We define the correction to the S-wave wavefunctions at the origin in the finite-r regularization by δΨ n (0) r 0 = − d 3 rĜ n (r , r)δV(r)Ψ LO n (r) where the subscript r 0 implies that the divergences are regulated by a finite distance r 0 between the Q andQ. Here, Ψ n (0)| r 0 = Ψ LO n (0) + δΨ n (0)| r 0 , and we used V LO (0) = 0 following the exact vanishing of the static potential at r = 0 [36].
In order to obtain the wavefunctions at the origin in the MS scheme, we need to compute the scheme conversion from finite-r regularization to DR. This is given by the difference between the two different schemes in the divergent integral d 3 rĜ n (0, r)δV(r)Ψ LO n (r). We define the scheme conversion coefficient δZ through the relation where the divergent integral in the MS scheme is first computed in momentum space in d = 4 − 2 spacetime dimensions, and then the UV poles are subtracted according to the prescription for renormalization in the MS scheme. The prescription that we use for the MS scheme is defined by associating a factor of (Λ 2 e γ E 4π ) with each loop integration, and subtracting the 1/ poles after evaluating the loop integrals in DR and expanding in powers of . Here, γ E is the Euler-Mascheroni constant. Then, Λ is the renormalization scale in the MS scheme. We note that in the calculation of the scheme conversion, we employ the potentials determined from on-shell matching in order to ensure the consistency with the SDCs computed in DR. Since Ψ LO n (r) is regular at r = 0, we can replace Ψ LO n (r) in the integrand by Ψ LO n (0) without affecting the right-hand side of eq. (3.19), because this affects only the finite parts of the divergent integrals, which cancel in δZ. Therefore, we can write δZ as δZ = d 3 rĜ n (0, r)δV(r) We compute δZ in the next section.
As we have argued based on the divergent behavior of the wavefunctions at the origin from nonrelativistic quantum mechanics, the corrections to the wavefunctions at the origin that are divergent at r 0 = 0 in the finite-r regularization can contain contributions that are not suppressed by powers of 1/m, even though the corrections come from 1/m and 1/m 2 potentials. For example, the logarithmically divergent correction from the 1/m potential is proportional to α 2 s log r 0 , which is not suppressed by any power of 1/m. This behavior will be confirmed in the calculation of δZ in the next section. The appearance of such large corrections is in accordance with the fact that, unless r m −1 , the 1/m and 1/m 2 potentials can overpower the static potential at short distances, so that the expansion in powers of 1/m is no longer valid. These corrections can potentially jeopardize the nonrelativistic power counting, unless they are subtracted through renormalization. Since δZ reproduces the divergent small r 0 dependence of the wavefunctions at the origin in the finite-r regularization, the divergences in Ψ n (0)| r 0 at small r 0 are subtracted completely by the scheme conversion δZ × Ψ LO n (0), and hence, the nonrelativistic power counting is restored in the MS-renormalized wavefunctions at the origin.

S-wave quarkonium wavefunctions at the origin in the MS scheme
In this section, we compute the scheme conversion coefficient δZ defined in eq. (3.20), which converts finite-r regularization to the MS scheme. We note that, since the reduced Green's function can be written as a linear combination of the Green's function G(r , r; E) for different E by using eq. (3.11), it is sufficient to compute (4.1) We will later show that δZ E is independent of E, and therefore, coincides with δZ for all S-wave states n. We first compute eq. (4.1) in perturbative QCD, and then show that the nonperturbative long-distance behaviors of the potentials do not affect the result.

Green's function in dimensional regularization
In order to compute the divergent integral in eq. where we use the shorthand Then, the divergent integral d 3 r G(0, r; E)δV(r) in DR can be written as where δṼ(p) is the momentum-space counterpart of δV(r) in d dimensions. Explicit expressions for δṼ(p) in DR will be given in the next section. The finite-r regularized integral can also be expressed in terms of the momentum-space Green's function as wheren is an arbitrary unit vector. Invariance of δṼ(p) under rotations ensures that the right-hand side of eq. (4.5) is independent ofn. The integrals in the finite-r regularization can be computed at d = 4 because the UV divergence is regulated by r 0 > 0. We note that the UV divergences in eqs. (4.4) and (4.5) come from the behavior of G(p , p; E) and δṼ(p) at large p and p . We compute the Green's function in momentum space in order to determine its large-momentum behavior. The momentum-space Green's function satisfies the Lippmann-Schwinger equation whereṼ LO (k) is the LO potential in momentum space. A formal solution of eq. (4.6) can be found iteratively, which reads where the first term comes from the free propagation of the QQ, and the second term corresponds to a single exchange of the LO potential between the Q and theQ. The quantity T encodes two or more exchanges of the LO potential to all orders: where k n+1 = p − p for each n. The formal solution in eq. (4.7) is organized so that the large p and p behavior in each term becomes less divergent as the number of exchanges of the LO potential increases [47]. This greatly simplifies the calculation of δZ, since the divergent contributions in eqs. (4.4) and (4.5) come only from the first few terms in eq. (4.7), and the non-divergent contributions coming from higher numbers of exchanges of the LO potential cancel in eq. (4.1). Hence, for the purpose of computing δZ, it suffices to consider only the first few terms in eq. (4.7).

Potentials in dimensional regularization
A necessary ingredient in computing δZ is the potentials in momentum space in d spacetime dimensions. In order to obtain the correct d-dimensional expressions, it is necessary to compute the potentials in the on-shell matching scheme in momentum space, which is done in perturbative QCD. The momentum-space potentialṼ (p , p) appears in the d-dimensional momentum-space Schrödinger equation in the form whereΨ n (p) is the momentum-space wavefunction. In d = 4 dimensions, the momentumspace potential is related to the position-space counterpart V (r, ∇) bỹ (4.10) The d-dimensional expression forṼ (p , p) to two-loop accuracy has been obtained in refs. [30,47], which we display here: where q = p − p, and the scale Λ comes from associating a factor of Λ 2 e γ E 4π with each loop integration. The constant c is defined by where in the last equality, we expanded in powers of up to order . The constant s depends on S; for spin triplet (S 2 = 2), and for spin singlet (S 2 = 0), The constant s for the spin singlet case can be obtained by using the results in ref. [47] to compute the d-dimensional spin projection according to the treatment of Pauli matrices in DR in ref. [48]. Equation (4.14) agrees through order with refs. [12,46]. To orderaccuracy, s + 2 can be written in terms of S 2 as Equation (4.11) implies that the LO potential in momentum space is given bỹ which is valid in d = 4 − 2 dimensions. The term δṼ C (q 2 ) corresponds to the loop corrections to the static potential, for which the explicit expressions can be found in refs. [30,47,49,50]. Since the corrections from δṼ C (q 2 ) to the wavefunctions at the origin are finite, we do not need to consider this term in the calculation of δZ. We note that eq. (4.11) reproduces the position-space expressions in eq. (B.3). Now we obtain the d-dimensional expression for δṼ(p) from eq. (4.11) by repeating the calculations in eqs. (3.13), (3.14), and (3.15) in momentum space in d spacetime dimensions. After a straightforward calculation, we obtain The first term corresponds to the 1/m potential, while the second term comes from the spin-dependent potential. The third and the fourth terms correspond to the corrections from the velocity-dependent potential and the relativistic correction to the kinetic energy, respectively. We see that at d = 4, the Fourier transform of δṼ(q) is exactly the positionspace counterpart δV(r) in eq. (3.16) at short distances.

Scheme conversion
Now we compute δZ using the d-dimensional momentum-space expressions of the Green's function in eq. (4.7) and δṼ(q) in eq. (4.17). We first compute the contribution to δZ E from the 1/m potential, which comes from the first term in eq. (4.17). The free propagation term in the formal solution for the Green's function gives the following contribution in DR: with the integral over p. Since the integral in eq. (4.18) is logarithmically divergent, the contributions from one or more exchanges of the LO potential in the Green's function are finite and do not contribute to δZ. The same quantity in finite-r regularization can be computed using the momentum-space expression in eq. (4.5), which gives The contribution to δZ E from the 1/m potential can be found by subtracting eq. (4.19) from eq. (4.18), and then subtracting the 1/ pole. Since the dependence on E cancels between the dimensionally regulated integral and the finite-r regularized integral, we obtain Here, Λ is now the MS scale associated with the renormalization of the wavefunction at the origin.
We now compute the contribution to δZ coming from the velocity-dependent potential, which comes from the third term in eq. (4.17). SinceṼ LO (k) = −4πα s C F /k 2 in perturbative QCD, the first term in the last line in eq. (4.17) can be evaluated as where again we associated a factor of Λ 2 e γ E 4π with the integral over k. Apart from an -dependent factor, eq. (4.21) is the same as the 1/m potential in perturbative QCD. Hence, Again, Λ is now the MS scale associated with the renormalization of the wavefunction at the origin. The computation of the contribution from the relativistic correction to the kinetic energy is similar to the case of the velocity-dependent potential. The last term in eq. (4.17) can be computed in perturbative QCD as This is just 1/4 times the result of the velocity-dependent potential in eq. (4.21). Hence, we obtain Finally, we consider the contribution from the spin-dependent potential in eq. (4.17). The free propagation term in the Green's function gives, in DR, This integral is power UV divergent; while this is not apparent in the last line of eq. (4.25) because power divergences are subtracted automatically in DR, the divergence can still be identified in the integrand. This implies that the contribution from the second term in eq. (4.7) is also divergent. This contribution reads, in DR, where we associated a factor of Λ 2 e γ E 4π with each loop integration. This integral is logarithmically divergent, and hence, the contributions from two or more exchanges of the LO potential are finite. Now we compute the divergent integrals in finite-r regularization, where the free propagating contribution gives The contribution from one exchange of the LO potential is The spin-dependent contribution to δZ can then be obtained from eqs. (4.25), (4.26), (4.27), and (4.28). We see again that the dependences on E cancel between the dimensionally regulated integrals and the finite-r regularized integrals. After subtracting the 1/ pole, we obtain We note that calculation of the spin-dependent contribution has been done in ref. [46]. Equation (45) of ref. [46] can be obtained by subtracting eq. (4.28) from eq. (4.26), dividing by a factor of πα s C F (s + 2)/m 2 , subtracting the 1/ pole, and setting Λ = e −γ E /r 0 . The complete result for δZ, which can be obtained by combining eqs.
where we use the shorthand L Λ = log(Λr 0 e γ E ). The scale Λ is the MS scale associated with the renormalization of the wavefunction at the origin. While the contribution from the spin-dependent potential has been obtained in ref. [46], the contributions from the 1/m potential, the velocity-dependent potential, and the relativistic correction to the kinetic energy are new. This result is accurate up to order α 2 s , and is also sufficiently accurate to reproduce the divergent small-r 0 behavior of the finite-r regularized wavefunctions at the origin Ψ n (0)| r 0 computed to first order in the QMPT using eq. (3.17). We note that in the calculation of δZ E , the E dependences cancel between the MSrenormalized integrals and the finite-r regularized integrals. This cancellation also occurs in the individual contributions in eqs. (4.20), (4.22), (4.24), and (4.29). We argue that this cancellation is not accidental; since a shift in E modifies only the finite piece of the divergent integral d 3 r G(0, r; E)δV(r), changes in E do not affect the difference between DR and finite-r regularization. Therefore, the scheme conversion coefficient δZ is given by eq. (4.30) for all S-wave states n.
The cancellation of the E dependence in δZ E follows from the fact that δZ E depends only on the behavior of the integrandG(p , p; E)δṼ(p) at large p and p . This implies that δZ is also unaffected by any modifications to the integrand that preserve the largemomentum behavior. We note that in position space, inclusion of the nonperturbative long-distance contribution to the potential can be done by adding functions of r that are regular at r = 0. In momentum space, this is equivalent to modifyingṼ LO (q) and δṼ(q) by adding functions of q that decrease faster than 1/q 2 at large q. Since such modifications do not affect the large-momentum behavior of the integrandG(p , p; E)δṼ(p), they do not affect δZ. As a result, eq. (4.30) is still valid, through order-α 2 s accuracy, even when the nonperturbative long-distance contributions are included in the potential. 1 The same argument can be made in position space: the short-distance divergence of the integral d 3 r G(0, r; E)δV(r) is determined completely in perturbative QCD, because the shortdistance behavior of δV(r) is unaffected by nonperturbative effects, and the divergent shortdistance behavior of the position-space Green's function G(r , r; E) is determined only by the short-distance behavior of V LO (r).
From the relation Ψ n (0)| MS = Ψ n (0)| r 0 −δZ ×Ψ LO n (0) we see that the scale dependence of Ψ n (0)| MS is determined by δZ, because Ψ n (0)| r 0 and Ψ LO n (0) do not depend on Λ. This allows us to compute the anomalous dimension of S-wave quarkonium wavefunctions at the origin. For spin triplet, we obtain and for spin singlet, The anomalous dimensions of the S-wave wavefunctions at the origin in eqs.

Unitary transformation
As we have discussed previously, different forms of the potential can be obtained by using unitary transformations. Since a different form of the potential can lead to a different behavior of the finite-r regularized wavefunctions at the origin Ψ n (0)| r 0 , the expression for δZ in eq. (4.30) is valid only when potentials from on-shell matching are used. On the other hand, if we want to include corrections from the nonperturbative long-distance behaviors of the potentials beyond leading order in 1/m, it is necessary to employ the nonperturbative definitions of the 1/m and 1/m 2 potentials from Wilson-loop matching. The wavefunctions computed in Wilson-loop matching must then be converted to wavefunctions in on-shell matching in order to compute the MS-renormalized wavefunctions at the origin using the relation in eq. (3.18). 2 In perturbative QCD, the explicit form of the unitary transformation that is necessary to obtain the potentials in Wilson-loop matching [eq. (B.4)] from the potentials in on-shell matching [eq. (B. 3)] has been derived in ref. [23,45]. If we define where the superscripts OS and WL imply that the potentials are obtained from on-shell matching and Wilson-loop matching, respectively. Since the differences in the 1/m and 1/m 2 potentials between on-shell matching and Wilson loop matching are known only at short distances, the precise form of W (r) can be obtained only near r = 0. Since the wavefunctions at the origin in finite-r regularization depend only on the wavefunctions at short distances, it suffices to determine U (r) for small r. We obtain where we neglect the correction from δV C (r). Since we consider δV C (r) as perturbations in the QMPT, the correction to W (r) from δV C (r) corresponds to a piece of the second order correction in QMPT from insertions of δV C (r) and V (1) (r)/m. Since we work at first order in the QMPT, we neglect this correction. Hence, to relative order 1/m, If Ψ OS n (r) is a solution of the Schrödinger equation with the potentials from on-shell matching, the wavefunction Ψ WL n (r) that satisfies the Schrödinger equation from Wilson loop matching is given by We note that near r = 0 this relation also holds for the finite-r regularized wavefunctions at the origin, because the finite-r regularized wavefunction at the origin reproduces the divergent small r 0 behavior of the wavefunction Ψ n (r) at |r| = r 0 . By using the relation in eq. (3.18), we can compute the unitary transformation of Ψ OS n (0)| r 0 as where the last equality follows from the fact that the difference between Ψ n (0)| MS and Ψ LO n (0) is suppressed by at least 1/m, because the divergent corrections in Ψ OS n (0)| r 0 that are not suppressed by powers of 1/m are subtracted by the scheme conversion. Equation (4.39) implies that We note that the difference between Ψ OS n (0)| r 0 and Ψ WL n (0)| r 0 does not depend on the longdistance behavior of the 1/m potential, up to corrections that are suppressed by 1/m 2 . This gives rise to the following approximate relation which is obtained by dividing eq. (4.40) by Ψ LO n (0) and taking the 1/m and 1/m 2 potentials given in perturbative QCD. We see from eqs. (4.40) and (4.41) that the difference between Ψ OS n (0)| r 0 and Ψ WL n (0)| r 0 comes only from the difference in the potential at short distances that is determined in perturbative QCD. This implies that, for the purpose of computing Ψ OS n (0)| r 0 including the nonperturbative long-distance behavior of the 1/m potential, it suffices to use the following prescription so that while the potential at short distances is given by the expressions from on-shell matching [eq. (B. 3)], the long-distance behavior is given by Wilson loop matching. We use eq. (4.42) to include the nonperturbative long-distance behavior of the 1/m potential in computing Ψ OS n (0)| r 0 . If we also wanted to include the long-distance nonperturbative behavior of the potentials of order 1/m 2 , we would have needed to include order 1/m 2 contributions in the unitary transformation that we neglected in eq. (4.37). This contribution, on the other hand, includes divergences at r = 0 that come from second order corrections to the wavefunctions at the origin. In turn, the second order corrections to the wavefunctions at the origin produce divergences at relative order α 3 s , which is beyond the accuracy of this paper. Hence, we neglect the long-distance nonperturbative behavior of the potentials of order 1/m 2 .

Numerical results
We now compute wavefunctions at the origin in the MS scheme for S-wave charmonia and bottomonia. We first compute the wavefunctions at the origin in the finite-r regularization from eq. (3.17), including the nonperturbative long-distance contribution to the static potential. We also consider the nonperturbative long-distance contribution to the 1/m potential by using the prescription given in eq. (4.42). Then, the MS-renormalized wavefunctions at the origin can be obtained from the relation in eq. (3.18), where δZ is given by eq. (4.30). The validity of this numerical procedure is verified by comparing with the known analytical results from perturbative QCD calculations in appendix D.
We compute decay constants and electromagnetic decay rates of S-wave charmonia and bottomonia based on the values of the wavefunctions at the origin that we obtain. The decay constants that we consider are defined in QCD by for a vector quarkonium V , where the Dirac spinor Q is a heavy quark field in QCD, m V is the mass of the quarkonium V , and is the polarization vector for the state |V . In the QCD definition of f V , the state |V is normalized relativistically. We also consider the decay constants of pseudoscalar quarkonium P , which are defined by where p µ is the 4-momentum of the quarkonium P , and the state |P is normalized relativistically. The decay constants f V and f P are renormalization scheme and scale independent. We list the NRQCD factorization formulas and the SDCs for the decay constants in appendix C. We note that the decay constants are given at leading order in α s and v by where m P is the mass of the quarkonium P . The decay constant f V is related to the leptonic decay rate of V by where α is the QED coupling constant, and e Q is the fractional charge of the heavy quark Q. For pseudoscalar quarkonium, f P cannot be compared directly with an experimentally measurable quantity, although it appears in hard exclusive production rates of pseudoscalar quarkonium P [52][53][54][55], and also have been studied in lattice QCD [56,57]. We also compute the two-photon decay rates of pseudoscalar quarkonia, which can be compared with measurements. We list the NRQCD factorization formula and the SDCs for the two-photon decay rate in appendix C.

Numerical inputs
We list the numerical inputs that we use in the numerical calculations in this section.

Heavy quark mass and the strong coupling
The heavy quark mass m that appears in the Schrödinger equation, as well as the scheme conversion coefficient in eq. (4.30) is the heavy quark pole mass, which suffers from renormalon ambiguity [58]. In order to avoid this issue, we use the modified renormalon subtracted (RS ) mass m RS , which is related to the pole mass m by [59] where ν f is the scale associated with the renormalon subtraction, and δm RS (ν f ) is given as a series in α s . At leading order in α s (order α 2 s ), δm RS (ν f ) is given by Here, the constant N m is known numerically as N m = 0.5626(260) [60], and the constants b and c k are determined by the QCD β function. Explicit formulas for b and c k are given in ref. [59]. We truncate the series in eq. (5.6) by including terms up to k = 2, which is equivalent to considering the running of α s at 4-loop accuracy. In principle, the QCD renormalization scale at which α s is computed in eq. (5.6) can be different from ν f , and the scale dependence is compensated by corrections of higher orders in α s . Because the leading renormalon ambiguity of order Λ QCD in the pole mass is not present in the RS mass, accurate values of the RS masses have been obtained for charm and bottom. In order to use the RS masses in calculation of the wavefunctions, we replace the heavy quark pole mass m by m RS + δm RS , and expand perturbatively in powers of δm RS . Since δm RS begins at order α 2 s , it suffices to consider only the correction coming from the kinetic energy in the Schrödinger equation. 3 The correction to the wavefunctions at the origin from the RS subtraction term is given by which is finite. The corrections associated with the 1/m and 1/m 2 potentials can be neglected, because they are suppressed by higher powers of α s . Hence, in computing the 3 The advantage of the use of the RS mass, compared to other renormalon subtraction schemes [59,61,62], is that the renormalon subtraction term δm RS (ν f ) begins at order α 2 s , rather than order αs. This makes it easier to organize the corrections from δm RS (ν f ) in powers of αs. For example, in the RS scheme, the corrections to the binding energies from δm RS (ν f ) contribute to decay constants and decay rates from order α 2 s v 2 , which we ignore at the current level of accuracy.
wavefunctions at the origin, it suffices to replace m by m RS (ν f ) everywhere, and then compensate for the difference by adding the correction term in eq. (5.7) to the wavefunctions at the origin, where δm RS is truncated at order α 2 s . We adopt the values for the RS charm and bottom masses determined in ref. [60] for ν f = 2 GeV, which are given by We compute α s in the MS scheme at a fixed QCD renormalization scale µ R , except when we consider resummation of logarithms in the loop corrections to the potentials. This is in order to facilitate exact order by order cancellation of the dependence on the factorization scale Λ in NRQCD factorization formulas, which requires the same α s to be used in the SDCs and in the calculations of the wavefunctions at the origin. Reference [60] gives ranges of µ R in which the theoretical determinations of η c and η b masses have mild dependences on the scale; these are given by µ R = 2.5 +1.5 −1.0 GeV for charm, and µ R = 5 +3 −3 GeV for bottom. We compute the numerical values of α s in the MS scheme at these scales using RunDec at 4-loop accuracy [63,64].

Static potential from lattice QCD
Precise nonperturbative determinations of the static potential V (0) (r) have been done in unquenched lattice QCD. We take the following parameterization of the static potential in ref. [65] given by where V 0 , σ, e, and g are determined in fits to lattice measurements. We take the results in the physical limit given in table V of ref. [65], where the central values read V 0 = 0.760/a, σ = (0.171/a) 2 , e = 0.368, with a −1 = 2.68 GeV. We ignore the last term in eq. (5.9), which comes from short distance lattice artifacts, because we only need the lattice QCD result at long distances. We match eq. (5.9) with the perturbative QCD expression in eq. (B.1) at r = r match , where the r dependences of the two expressions agree well. While the renormalization scale dependence in V (0) (r) pert cancels order by order in perturbative QCD, it is known that the convergence of the r dependence of V (0) (r) pert is poor when the renormalization scale is fixed [46]. The convergence can be improved if we resum the logarithms that are associated with the running of α s , which can be done by choosing the renormalization scale to be proportional to 1/r at short distances. In order to avoid the renormalization scale being too small, we set the r-dependent renormalization scale to be µ r = (r −2 + µ 2 R ) 1/2 , so that µ r > µ R , while µ r ≈ 1/r at short distances.
We define the nonperturbative long-distance contribution to the static potential as where V (0) (r) pert, resum is computed at the renormalization scale µ r , and ∆V (0) is chosen so that the right-hand side vanishes at r = r match . We choose r −1 match = 1.5 GeV, which is where the slopes of V (0) (r) lattice and V (0) (r) pert, resum are approximately same. Since V (0) (r) long vanishes for r < r match , we obtain the following expression for the static potential that is valid for both short and long distances: so that V (0) (r) coincides with the perturbative QCD expression for r < r match , while it reproduces the lattice QCD determination for r > r match . In fig. 1 we compare the unquenched lattice QCD results in ref. [65] with the expression for V (0) (r) in eq. (5.11).
The perturbative QCD expressions of the static potential depends on the number of light quark flavors n f , which we take to be n f = 3 for charm, and n f = 4 for bottom. We note that the matching of perturbative QCD and lattice QCD for the pNRQCD potentials have been done in a similar way in refs. [66,67] for heavy quarkonium spectroscopy. Eq. (5.11) implies that the leading-order potential is given by where in the first term on the right-hand side, α s is evaluated at a fixed renormalization scale µ R . The Coulombic correction term δV C (r) = V (0) (r) − V LO (r) that appears in the second line of eq. (3.17) is given by where in the last term, α s is evaluated at a fixed renormalization scale µ R . We note that δV C (r) contain contributions whose net effect is to shift the Coulomb strength of the LO potential at relative order α s , and so, the correction to the wavefunctions at the origin from δV C (r) begins at relative order α s . Although we work through first order in the QMPT, second order corrections from δV C (r) might be important, because this is of order α 2 s . Computing the second order Coulombic correction can also be useful in testing the convergence of the Coulombic corrections. The second order Coulombic correction to the wavefunctions at the origin is given in terms of the reduced Green's functions by In computing the second order Coulombic corrections, we neglect the a 2 term in the static potential [eq. (B.1)], because at second order in the QMPT, this term contributes at relative order α 3 s .

1/m potential from lattice QCD
We use a similar strategy as the previous section to determine the nonperturbative longdistance contribution to the 1/m potential. Unlike the static potential, nonperturbative determinations of the 1/m potential is available only from quenched lattice QCD. We use the parametrization in ref. [68] given by where A = 0.297 and σ (1) = 0.142 GeV 2 . This parametrization, which is based on the longdistance behavior expected from effective string theory in ref. [69], is obtained in ref. [68] from quenched lattice QCD results at lattice coupling β = 6.0. We match eq. (5.15) with the perturbative QCD expression at r = r match . Since we do not consider loop corrections to the 1/m potential, the expression for the 1/m potential V (1) (r) WL pert at leading order in α s depends on the choice of the renormalization scale. In order to reduce the scale dependence, we resum the logarithms that are associated with the running of α s by setting with µ r = (r −2 + µ 2 R ) 1/2 . We define the nonperturbative long-distance contribution to the 1/m potential by where ∆V (1) is chosen so that the right-hand side vanishes at r = r match . We choose r −1 match = 1.5 GeV. Since V (1) (r) WL long vanishes for r < r match , we obtain an expression for the 1/m potential that is valid for both short and long distances given by We compare the lattice QCD determination in eq. (5.15) with the expression for V (1) (r) WL in eq. (5.18) in fig. 2. Based on the argument given in section 4.4, we obtain the expression for the 1/m potential in on-shell matching that is valid for computation of wavefunctions at the origin, given by where in the first term on the right-hand side, α s is computed at a fixed renormalization scale µ R . We use this form of the 1/m potential in the calculation of the wavefunctions at the origin.

Reduced Green's function
We compute the reduced Green's functionĜ n (r , r) numerically by using two different methods, which are valid in different regimes of r and r . In the first method, which is valid for small r and r , we compute the Green's function in position space numerically by using the method given in ref. [70]. We only need to compute the S-wave contribution, which is defined by including only the S-wave states in the sum in eq. (3.8). This contribution can be written as where r < = min(|r|, |r |), r > = max(|r|, |r |), and the superscript S denotes the S-wave contribution. The functions u < and u > are two independent solutions of the differential equation with the following boundary condition so that u < (r)/r is regular at r = 0, while u > (r) is square integrable. We determine the functions u < and u > by numerically solving the differential equation for a given E. The reduced Green's function can then be obtained by using the relation in eq. (3.11), where we take the limit numerically. We note that, if E coincides with an eigenenergy of the LO Schrödinger equation E LO n , then the corresponding wavefunction Ψ LO n (r) is proportional to u < (r)/r. This means that u < (r) is square integrable if E = E LO n , and in such case, the square-integrable solution u > (r) does not exist. Hence, the limit in eq. (3.11) must be taken with care, because the numerical solution for u > (r) becomes unstable if E is too close to E LO n . When we compute the reduced Green's functions numerically using eq. (3.11), we set η = 10 −3 GeV.
Since the first method involves computing u < (r < ) by solving a differential equation with initial conditions at r < = 0, the method becomes unreliable when r and r are both large. For large r and r , we compute the reduced Green's function by using the formal solution in eq. (3.6), where we truncate the series by including only a limited number of the lowest eigensolutions of the LO Schrödinger equation. This method in turn becomes unreliable at small r and r . For example, if the LO potential V LO (r) is linear in r at long distances, the eigenenergies of highly excited S-wave states increase linearly with increasing principal quantum number, and the LO wavefunctions at the origin are constant in the principal quantum number. Hence, the series in eq. (3.6) diverges like ∞ n 1/n at r = r = 0. This implies that the truncated series becomes unreliable at small r and r .
We combine the reduced Green's function at long and short distances bŷ whereĜ n (r , r)| short is computed by using eqs. (5.20) and (3.11),Ĝ n (r , r)| long is computed by truncating the series in eq. (3.6), and b(r) is a smooth function that satisfies b(0) = 1 and b(∞) = 0, so that eq. (5.23) is reliable for all r and r . We define b(r) by with r b = 1 GeV −1 . The validity of the reduced Green's function obtained in eq. (5.23) can be tested by numerically checking the relations for k = n. We note that, due to the boundary condition u > (0) = 1, it is evident that G S (0, r; E) develops a power divergence given by m/(4πr) near r = 0. It has been shown in ref. [46] that if the LO potential is given by V LO (r) = −α s C F /r at short distances, u > (r)/r also contains a logarithmic divergence given by −α s C F m log r. Therefore, near r = 0, the Green's function behaves like where the ellipsis represent contributions that are finite at r = 0. This shows that the divergent small r behavior of G S (0, r; E) depends only on the short-distance behavior of the LO potential, which is determined in perturbative QCD.

Gluonic correlators
The pNRQCD expressions of the NRQCD LDMEs in eqs. (2.4) and (2.8) depend on gluonic correlators that scale with powers of Λ QCD . Also, corrections to the wavefunctions at the origin from the velocity-dependent potential involve V p 2 (0), which in DR, is proportional to the correlator iE 2 . While the gluonic correlators of mass dimension two contribute to the NRQCD LDMEs at relative order Λ 2 QCD /m 2 , the dimensionless correlator E 3 contributes to 0|χ † · σψ|V and 0|χ † ψ|P at relative order v 2 , and the correlator iE 2 contributes to the wavefunctions at the origin at relative order Λ QCD /m.
The dimensionless correlator E 3 in the MS scheme has been determined in ref. [37] from measured decay rates of P -wave charmonia. At the MS scale Λ = 1 GeV, The correlator E 3 depends logarithmically on the scale. We compute E 3 at other scales by using the one-loop renormalization group improved expression [37,39] Reference [37] also provides a determination of iE 2 from measured electromagnetic decay and production rates of P -wave charmonia. However, the determination in ref. [37] has uncertainties that are larger than the typical size of the correlator that is expected from its power counting. For this reason, instead of taking the determination in ref. [37], we consider the effect of V (2) p 2 (0) to the wavefunctions at the origin in the uncertainties by assuming |V (2) p 2 | 500 MeV, which corresponds to the typical size of Λ QCD . Since the gluonic correlators of mass dimension two contribute to the NRQCD LDMEs at relative order Λ 2 QCD /m 2 , we neglect them in calculations of the LDMEs compared to corrections of relative order Λ QCD /m and v 2 .

Numerical results for S-wave charmonia
In this section, we compute the MS-renormalized wavefunctions at the origin for the 1S and 2S charmonium states. We identify the J/ψ and η c as the 1S charmonium states in spin-triplet and spin-singlet states, respectively, while the ψ(2S) and η c (2S) states are the 2S charmonium states in spin-triplet and spin-singlet states, respectively.
As we discussed in previous sections, we solve the Schrödinger equation numerically with the LO potential in eq. (5.12) and the charm quark mass in eq. (5.8a) to determine Ψ LO n (r), E LO n , andĜ n (r , r). Then, we compute the corrections to the wavefunctions at the origin in the finite-r regularization using eq. (3.17). In order to compensate for the use of the RS mass, we add to eq. (3.17) the finite correction from the RS subtraction term in eq. (5.7). We also add to eq. (3.17) the Coulombic correction at second order in QMPT in eq. (5.14). In the calculation of the corrections to the wavefunctions at the origin, we use the 1/m potential given by eq. In computing the finite-r regularized wavefunctions at the origin, the regulator r 0 must be chosen to be as small as possible, as long as the numerical calculation is stable. We determine an optimal choice of r 0 by numerically testing the approximate relation in eq. (4.41). We find that the relation is well reproduced numerically at 1% level for r 0 0.1 GeV −1 . Hence, we choose r 0 = 0.2 GeV −1 , and vary r 0 between 0.1 GeV −1 and 0.3 GeV −1 . We set the MS scale Λ to be the charm quark mass m, and choose the central value of the QCD renormalization scale µ R to be 2.5 GeV, as discussed in sec. 5.1.1.
We list the central values of the LO wavefunctions at the origin and the LO binding energies in table 1. We also list the corrections to the wavefunctions at the origin relative to Ψ LO n (0) in table 1. We classify the corrections by their origins in the following way: the non-Coulombic correction δ NC Ψ comes from the 1/m and 1/m 2 potentials, the Coulombic corrections δ C1 Ψ and δ C2 Ψ come from δV C (r) at first and second order in the QMPT, respectively, and the correction δ RS Ψ comes from the RS subtraction term. The explicit expressions for δ NC Ψ and δ C1 Ψ are given by while δ C2 Ψ is given by dividing eq. (5.14) by Ψ LO n (0), and δ RS Ψ is given by dividing eq. (5.7) by Ψ LO n (0). The MS-renormalized wavefunctions at the origin are then given by  Figure 3. The non-Coulombic corrections δ NC Ψ at finite r 0 for the charmonium 1S (solid lines) and 2S (dashed lines) states, for spin triplet (black) and spin singlet (red). The r 0 dependences are mild for the range 0.1 GeV −1 < r 0 < 0.3 GeV −1 that we consider. We note that the r 0 dependence cancels in δ NC Ψ between δZ and the finite-r regularized integral for small r 0 . We demonstrate this cancellation of the r 0 dependence in fig. 3.
The results for the LO binding energies for the 1S and 2S states in table 1 are roughly compatible with the mass difference between J/ψ and ψ(2S). We see that the non-Coulombic corrections from 1/m and 1/m 2 potentials, given by δ NC Ψ in table 1, are large and positive for both 1S and 2S states. This is in contrast with the order-α 2 s corrections to the NRQCD SDCs in appendix C, which are large and negative at Λ = m. This implies that if we combine the pNRQCD expressions of the LDMEs with the NRQCD SDCs, large cancellations will occur between the order-α 2 s corrections to the SDCs and the corrections to the wavefunctions at the origin. The contribution from the long-distance part of the 1/m potential, which is given by the second term in eq. (5.19), amounts to about +10% of the LO wavefunction at the origin for the 1S state, and about +6% of the LO wavefunction at the origin for the 2S state. We note that while the Coulombic corrections at first order are positive, the Coulombic corrections at second order are small and negative, signaling good convergence of the Coulombic corrections. The corrections from the renormalon subtraction term δ RS Ψ are mild for both 1S and 2S states. We use the results for the MS wavefunctions at the origin in table 1 to compute decay constants and electromagnetic decay rates of S-wave charmonium states. We first compute the decay constants f V of V = J/ψ and ψ(2S). By using the pNRQCD expressions of the LDMEs in eqs. (2.4) and (2.5) and the SDCs in appendix C, and expanding the corrections to the SDCs and to the wavefunctions at the origin, we obtain . This expression is valid up to corrections of relative order α 3 s , v 3 , and Λ 2 QCD /m 2 . We set n f = 3 in the SDCs. Since we assume δ C1 Ψ to be of order α s , we keep the cross term α s c v , while the Λ dependence in the order-α s correction to d v cancels with the scale dependence of the correlator E 3 . Hence, variation of the factorization scale Λ has almost no effect in eq. (5.31). We set the scale Λ = m in the one-loop correction to d v , and compute the correlator E 3 at the same scale using the renormalization group improved expression in eq. (5.28). We take the measured quarkonium masses from ref. [71].
The numerical result for the J/ψ decay constant is

+0.003
−0.000 ± 0.069 ± 0.054 GeV = 0.363 +0.089 −0.088 GeV, (5.32) where the first uncertainty comes from varying µ R between 1.5 GeV and 4 GeV, and the second uncertainty comes from varying r 0 between 0.1 GeV −1 and 0.3 GeV −1 . The third uncertainty comes from the neglect of the correction −V (2) p 2 (0)/(2m) to the wavefunction at the origin, which we take to be ±500 MeV/(2m) times the central value. The last uncertainty comes from the uncalculated corrections of order v 3 , which we take to be 15% of the central value. In the last equality, we add the uncertainties in quadrature.
We note that the central value for f J/ψ that we obtain is very close to the leading-order v , δ NC Ψ , δ C2 Ψ , δ RS Ψ , and the cross term α s c Ψ add up to about +14% of the central value, so that the numerical result for f J/ψ in eq. (5.32) is only about 1% larger than the leading-order value. If we had ignored the corrections to the wavefunctions at the origin, the one-loop correction would have been −23% of the leading order result, while the two-loop correction would have been −39% of the leading order value, so that the loop corrections add up to −62% at two-loop accuracy. The inclusion of the corrections to the wavefunction at the origin in the calculation of the decay constant f J/ψ has substantially improved the convergence of the expansion in α s and v.
The result for f J/ψ agrees within uncertainties with the lattice QCD determination using relativistic charm quarks in ref. [57], which gives f J/ψ = 0.4104 (17) GeV. In order to compare with measurements, we compute the leptonic decay rate of J/ψ from f J/ψ by using eq. (5.4). We obtain Γ(J/ψ → e + e − ) = 4.5 +2.5 −1.9 keV, (5.33) where we used α = 1/133, which is computed at the scale of the J/ψ mass. This result agrees with the experimental value Γ(J/ψ → e + e − ) = 5.53 ± 0.10 keV in ref. [71]  where we used α = 1/133, which is computed at the scale of the ψ(2S) mass. This result agrees with the experimental value Γ(ψ(2S) → e + e − ) = 2.33 ± 0.04 keV in ref. [71] within uncertainties. Now we compute the decay constants f P of P = η c and η c (2S). Although f P cannot be obtained directly from experimental measurements, this decay constant appears in exclusive production cross sections of pseudoscalar quarkonia at high energies [54,55]. We obtain the following expression for f P by expanding the corrections to the SDCs and the corrections to the wavefunctions at the origin: p + O(α 3 s ), and α s = α s (µ R ). We neglect the small imaginary part in c (2) p , which amounts to less than 0.02. This expression is valid up to corrections of relative order α 3 s , v 3 , and Λ 2 QCD /m 2 . We set n f = 3 in the SDCs. We note that the Λ dependence in δ NC Ψ | S 2 =0 cancels exactly with α 2 s c (2) p . Since the order-α s correction to d p is not available, our expression for f P in eq. (5.36) depends mildly on Λ through the scale dependence of E 3 . 4 Nevertheless, variation of the factorization scale Λ has a very small effect in eq. (5.36). We compute E 3 at the scale Λ = m. We take the measured quarkonium masses from ref. [71].
The numerical result for f ηc is where the uncertainties are as in eq. (5.32). We add the uncertainties in quadrature. We neglect the uncertainty from the scale dependence of E 3 , which is small compared to other uncertainties. This result for f ηc agrees with the lattice QCD determination using relativistic charm quarks in refs. [57], which gives f ηc = 0.3981 (10)  where the uncertainties are as in eq. (5.37). For the case of η c (2S), the central value for f ηc(2S) that we obtain is about 18% smaller than the leading-order value f LO ηc(2S) = 0.321 GeV. The difference is larger than the case of η c , because the binding energy of the 2S state is larger than the binding energy of the 1S state, and so, the correction proportional to E LO ηc(2S) is larger compared to the η c case. Finally, we compute the two-photon decay rate of P = η c and η c (2S). The NRQCD factorization formula for the decay rate is given in appendix C. The following expression for the decay rate is obtained by expanding the corrections to the SDCs and the corrections to the wavefunctions at the origin at the amplitude level: γγ + O(α 3 s ), and α s = α s (µ R ). This expression is valid up to corrections of relative order α 3 s , v 3 , and Λ 2 QCD /m 2 . We set n f = 3 in the SDCs, and use e Q = 2/3 for charm. We choose α = 1/137, because the QED coupling constant in eq. (5.39) is associated with on-shell photons in the final state. The dependence on the MS scale Λ in δ NC Ψ | S 2 =0 cancels completely with the Λ dependence in α 2 s c (2) γγ , while the Λ dependence in the order-α s correction to d γγ cancels with the scale dependence of the correlator E 3 . Similarly to the case of decay constants, variation of the factorization scale Λ has almost no effect in eq. (5.39). We set the scale Λ = m in the one-loop correction to d γγ , and compute the correlator E 3 at scale m using the renormalization group improved expression in eq. (5.28).
The numerical result for the two-photon decay rate of η c is where the uncertainties are as in eq. (5.32). This result is compatible within uncertainties with the PDG value of the two-photon decay rate Γ(η c → γγ) = 5.06 ± 0.34 keV in ref. [71]. The central value for the decay rate that we obtain is not very different from the leading-order value Γ(η c → γγ)| LO = 16N c πα 2 e 4 Q |Ψ LO ηc (0)| 2 /m 2 ηc = 6.0 keV. If we had ignored the corrections to the wavefunctions at the origin, the one-loop correction would have been −15% of the leading order result at the amplitude level, while the two-loop correction would have been −46% of the leading order amplitude, so that the effect of the loop corrections would add up to −61% at two-loop accuracy at the amplitude level. In contrast, the order-α s corrections in eq. (5.39) from α s c (1) γγ and δ C1 Ψ amount to about 3%, and the corrections proportional to E LO ηc is about −10% compared to the leadingorder amplitude. The remaining corrections from α 2 s c Ψ add up to about +14% of the leading-order amplitude, so that the central value for the decay rate that we obtain is about 7% larger than the leading-order value at the amplitude level, or about 13% larger at the squared amplitude level. Inclusion of the corrections to the wavefunctions at the origin in eq. (5.39) has improved the convergence of the corrections of higher orders in α s and v.
The numerical result for the two-photon decay rate of η c (2S) is where the uncertainties are as in eq. (5.32). The central value of the decay rate is about 19% smaller than the leading-order value Γ(η c (2S) → γγ)| LO = 16N c πα 2 e 4 Q |Ψ LO ηc(2S) (0)| 2 /m 2 ηc(2S) = 3.7 keV. We note that the result for the decay rate in eq. (5.41) disagrees with existing experimental values of the decay rate in refs. [72,73], which disagree with each other.

Numerical results for S-wave bottomonia
Now we compute the MS-renormalized wavefunctions at the origin for the 1S, 2S, and 3S bottomonium states. We identify the spin-triplet states as Υ(nS), while the spin-singlet states correspond to η b (nS), where n = 1, 2, and 3. The calculations for bottomonia are done similarly as the calculations for charmonium states, except that we use the bottom quark RS mass for m, set n f = 4, and choose the central value of the QCD renormalization scale to be µ R = 5 GeV. The range for r 0 is again determined from numerically checking the approximate relation in eq. (4.41). We choose the central value for r 0 to be r 0 = 0.1 GeV −1 , and vary r 0 between 0.05 GeV −1 and 0.2 GeV −1 . We set the MS scale Λ to be the bottom quark mass m. We use the measured masses of the Υ(nS) and η b (nS) states from ref. [71]. Because the mass of the η b (3S) state has not been measured, we estimate m η b (3S) by m Υ(3S) − (m Υ(2S) − m η b (2S) ), assuming that the hyperfine splitting is same for the 2S and 3S states. The r 0 dependences are mild for the range 0.05 GeV −1 < r 0 < 0.2 GeV −1 that we consider. We list the numerical results for the LO wavefunctions at the origin, the LO binding energies, and the corrections to the wavefunctions at the origin relative to the LO wavefunctions at the origin in table 2. The relative corrections δ NC Ψ , δ C1 Ψ , δ C2 Ψ , and δ RS Ψ are defined in the previous section. We show the r 0 dependence of the non-Coulombic correction δ NC Ψ in fig. 4.
The results for the LO binding energies for the 1S, 2S and 3S states in table 2 are roughly compatible with the mass differences between Υ(1S), Υ(2S), and Υ(3S) states. We see that the non-Coulombic corrections from 1/m and 1/m 2 potentials, given by δ NC Ψ in table 2, are large and positive for the 1S, 2S, and 3S states, although the relative sizes of the corrections are smaller than the case of charmonium 1S and 2S states. As it was in the case of charmonia, the order-α 2 s corrections to the SDCs in appendix C are large and negative at Λ = m, so that if we combine the pNRQCD expressions of the LDMEs with the SDCs, large cancellations will occur between the order-α 2 s corrections to the SDCs and the corrections to the wavefunctions at the origin. The contribution from the long-distance part of the 1/m potential, originating from the second term in the expression for the 1/m potential in eq. (5.19), amounts to about +4%, +2%, and +2% of the LO wavefunction at the origin for the 1S, 2S, and 3S states, respectively, which are less than half of the corresponding corrections to the charmonium wavefunctions at the origin. We note that while the Coulombic corrections at first order are positive, the Coulombic corrections at second order are small and negative, signaling good convergence of the Coulombic corrections. The corrections from the renormalon subtraction term δ RS Ψ are small. Now we compute the decay constants f Υ(nS) , f η b (nS) , and the electromagnetic decay rates of Υ(nS) and η b (nS) based on the bottomonium wavefunctions at the origin that we obtained. We use the same pNRQCD expressions for these quantities in eqs. (5.31), (5.36), and (5.39) that we used in the previous section for the charmonium states, except that we set n f = 4 in the SDCs, and use e Q = −1/3 for bottom. We note that the correction terms in the pNRQCD expressions of the LDMEs in eqs. (2.4) and (2.8) that come from the gluonic correlators may not be valid for 1S bottomonium states, because the assumption mv Λ QCD mv 2 may not hold for these states. Hence, when we make predictions for the bottomonium 1S states, we assume that eqs.
p 2 (0)/(2m) to the wavefunctions at the origin, which we take to be ±500 MeV/(2m) times the central value. For f Υ(1S) , the final uncertainty comes from the uncalculated order-v 2 corrections to the LDME, which we take to be 10% of the central value. For f Υ(2S) and f Υ(3S) , the final uncertainties come from the uncalculated corrections of order v 3 , which we take to be 3% of the central value. We add the uncertainties in quadrature.  [71]. We note that the result for Γ(Υ(1S) → e + e − ) that we obtain also agrees well with the perturbative QCD prediction at third order in ref. [75]. However, the convergence of the perturbative expansion in the perturbative QCD calculation is poor; according to ref. [75], the size of the corrections at first and second order, when combined, exceeds the leadingorder result, while the third order correction is moderate. In contrast, in the calculation of the decay rate Γ(Υ(1S) → e + e − ) in this work, the loop corrections to the SDCs and the corrections to the wavefunction at the origin combine to be 25% of the leading order result at the squared amplitude level. It seems that the improvement of the convergence has mostly to do with the Coulombic corrections, because the non-Coulombic correction δ NC Ψ does not change much from the result in table 2 when we neglect the long-distance, nonperturbative part of the LO potential given by the second term in eq. where the uncertainties are as in f Υ(nS) . We add the uncertainties in quadrature. Compared to the LO calculation of the decay rates, the corrections from loop corrections to the SDCs and the corrections to the wavefunctions at the origin combine to be about 58%, 15%, and 3% for the η b (1S), η b (2S), and η b (3S) states, respectively. At the amplitude level, the corrections amount to about 26%, 7%, and 2% for the η b (1S), η b (2S), and η b (3S) states, respectively. If we had ignored the corrections to the wavefunctions at the origin, the orderα s correction would have been −11%, and the order-α 2 s correction would have been −26% of the central value at the amplitude level, so that the loop corrections to two-loop accuracy would add up to −38% of the leading-order amplitude. By the inclusion of the corrections to the wavefunctions at the origin, the sizes of the corrections are reduced considerably for the η b (2S) and η b (3S) states, while the improvement is moderate for the η b (1S) state.

Summary and discussion
In this paper, we have computed the wavefunctions at the origin of S-wave heavy quarkonia in the MS renormalization scheme. We include the nonperturbative long-distance contributions to the potential, which are neglected in perturbative QCD calculations. We compute corrections to the wavefunctions at the origin at subleading orders in 1/m in position space, where the ultraviolet divergences are regulated by using finite-r regularization. The position-space expressions for the corrections to the wavefunctions at the origin are given in section 3. The wavefunctions at the origin in finite-r regularization is then converted to the MS scheme by computing the scheme conversion in perturbative QCD. The result for the scheme conversion coefficient is given in section 4. We use the results for the wavefunctions at the origin to make first-principles based, model-independent predictions of decay constants and electromagnetic decay rates of S-wave charmonium and bottomonium states in section 5.
The predictions for the leptonic decay rates of J/ψ, ψ(2S), η c , and Υ(nS) states in this work agree with experimental measurements within uncertainties. The predictions for the J/ψ and η c decay constants agree within uncertainties with lattice QCD calculations in ref. [57], which make use of relativistic charm quarks. The predictions for the Υ(1S) and Υ(2S) decay constants agree with the lattice NRQCD determinations in ref. [74], where lattice regularization is used to compute both the short-distance coefficients and the NRQCD matrix elements.
The calculation of the wavefunctions at the origin in this work contains several improvements compared to existing model dependent methods. First of all, in this work, we include potentials at leading and subleading orders in 1/m, which are determined by perturbative QCD at short distances, while their nonperturbative behaviors at long distances are fixed by lattice QCD. Secondly, the ultraviolet divergences that appear in corrections to the wavefunctions at the origin are properly renormalized in the MS scheme, so that the wavefunctions at the origin that we obtain have the correct scale dependences that are expected from perturbative QCD. Finally, the ambiguity in the heavy quark pole mass is removed by the use of the modified renormalon subtracted mass, whose numerical values are accurately known. These improvements are generally not possible in potential-model calculations.
Because the wavefunctions at the origin that we have computed have the correct scale dependence that reproduce the anomalous dimensions of NRQCD LDMEs, the dependences on the NRQCD factorization scale cancel completely through two-loop order in the pN-RQCD expressions for the decay constants and electromagnetic decay rates. Together with the order-by-order cancellation of the dependence on the QCD renormalization scale, this greatly reduces the uncertainty associated with scale dependences.
A novel feature of the calculation of the decay constants and decay rates in this work is that large cancellations occur between the corrections to the wavefunctions at the origin and the perturbative corrections to the short-distance coefficients. These cancellations substantially improve the convergence of the expansion in α s and v. This may have important implications in understanding the appearance of large perturbative corrections in calculations of short-distance coefficients in the MS scheme. A possible explanation of the cancellations is that, due to the confining nature of the nonperturbative potentials, including the long-distance contributions to the potentials in calculating the wavefunctions at the origin may have the effect of introducing an infrared cutoff, so that the renormalon ambiguities associated with the infrared contributions of loop corrections in perturbative QCD are resolved.
The pNRQCD expressions of the wavefunctions at the origin, as well as the decay constants and decay rates, depend on gluonic correlators, whose values are in general not well known. Especially, the correction from the velocity-dependent potential at zero distance, which is given by a gluonic correlator whose size is of order Λ QCD , is the largest source of uncertainties, with the exception of the bottomonium 1S states. Improved determinations of the gluonic correlators, which can in principle be done in lattice QCD, will be necessary in further reducing the uncertainties.
The calculation of the renormalization of the wavefunctions at the origin in this work is accurate to two-loop accuracy. In principle, the calculation in this work can be extended to three-loop accuracy, by computing the divergent corrections to the wavefunctions at the origin to second order in the quantum-mechanical perturbation theory, and computing the scheme conversion from finite-r regularization to the MS scheme to order-α 3 s accuracy. Such a calculation will make possible the inclusion of the long-distance nonperturbative contributions to the potentials of order 1/m 2 , because the second-order correction in the quantum-mechanical perturbation theory is necessary in extending the calculation of the unitary transformation between on-shell matching and Wilson-loop matching in section 4 to order-1/m 2 accuracy.
Finally, we note that the calculation in this paper may be extended to states with nonzero orbital angular momentum. This necessarily involves considering the orbital angular momentum dependent terms in the potential, as well as the angular dependence of the wavefunctions in dimensional regularization, which were not present in this work thanks to the rotational symmetry of the S-wave states. Such calculations will allow us to make accurate predictions of production and decay rates of P -wave heavy quarkonium states. Feynman diagrams for one-loop corrections to the NRQCD LDMEs. Solid lines are heavy quarks and antiquarks, dashed lines are temporal gluons, and curly lines are transverse gluons. Open circles represent insertions of the p·A vertex, and filled squares represent the operator χ † · σψ for spin triplet, and χ † ψ for spin singlet.

Acknowledgments
The author is grateful to Nora Brambilla and Antonio Vairo for fruitful discussions and their encouragement in completing this work. This work is supported by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) cluster of excellence "ORIGINS" under Germany's Excellence Strategy -EXC-2094 -390783311.

A Anomalous dimensions
In this appendix, we compute the anomalous dimensions of the NRQCD LDMEs that are given in eqs. (2.3) and (2.7). Although the results are available in refs. [4,[10][11][12][13], it can be useful to compute them through loop calculations in NRQCD, because such calculations can reveal the origins of the anomalous dimensions which are obscure in perturbative QCD calculations of SDCs.
The anomalous dimensions can be computed as perturbation series in α s by replacing the quarkonium states in the definitions of the LDMEs by perturbative QQ states, and computing loop corrections to the LDMEs, with the vertices coming from the operators in the NRQCD Lagrangian. The Q andQ in the QQ states are on shell, which have nonrelativistic 4-momenta (E, q) and (E, −q), respectively, where E = q 2 /(2m). We work in Coulomb gauge, and use the NRQCD Feynman rules given in ref. [76]. We use DR in d = 4 − 2 spacetime dimensions, where the anomalous dimensions are simply given by the coefficients of the 1/ poles that are associated with UV divergences.
The NRQCD loop integrals are evaluated in the following way. First, we integrate over the temporal components of the loop momenta, using contour integration. Then, we expand the integrand in powers of 1/m, which is necessary in preserving the nonrelativistic power counting in DR. Finally, we integrate over the spatial components of the loop momenta, regulating the resulting divergences in DR. The anomalous dimension is given by the coefficients of the single UV poles, after differentiating and multiplying by g s = √ 4πα s .
A.1 One-loop anomalous dimension at relative order v 2 We first consider the one-loop corrections to the NRQCD LDMEs 0|χ † · σψ|QQ and 0|χ † ψ|QQ , which come from the Feynman diagrams in fig. 5. Because the one-loop integrals that we compute only involve single poles in , we may set = 0 in the loop integrands without affecting the 1/ poles.
The one-loop heavy quark self energy from the first two diagrams in fig. 5 reads where the first and second terms come from exchange of temporal and spatial gluons, respectively. The first term is scaleless power divergent, and therefore can be discarded. Since the second term already has a factor of q 2 /m 2 , we can expand in powers of 1/m and keep only the leading contribution. To find the one-loop correction to the quark field renormalization factor Z NRQCD Q , we differentiate Σ(E, q) by E and set E = q 2 /(2m) to obtain where we only keep the UV pole in the last equality. We use the subscript UV to make clear that the pole is associated with a UV divergence. The vertex correction diagram from exchange of a temporal gluon gives which does not have a UV divergence, and hence does not contribute to the anomalous dimension. The transverse-gluon exchange diagram gives where A = (k + q) 2 /(2m) − q 2 /(2m) − iε. Here, the first term comes from the residue of the pole from the quark propagator, and the second term comes from the residue of the pole from the gluon propagator. Since eq. (A.4) already has a factor of q 2 /m 2 , we can keep only the leading contribution in the 1/m expansion, which gives where we only keep the UV pole. We note that the first term on the left-hand side does not have a UV divergence, and the UV pole comes only from the second term. It can be shown that if we replace the p·A vertices with σ ·B vertices, the transverse-gluon exchange diagram does not produces logarithmic UV divergences. Since the diagrams in fig. 5 give rise to logarithmic UV divergences at relative order α s v 2 already at leading power in the 1/m expansion, it is not necessary to consider vertices from higher dimensional operators in the NRQCD Lagrangian. We combine eqs. (A.2) and (A.5) to find the UV pole in the one-loop correction to the NRQCD LDME 0|χ † ·σψ|QQ , which reads Feynman diagrams for two-loop corrections to the NRQCD LDMEs that produce logarithmic UV divergences. There are additional diagrams that can be obtained from charge conjugation, which we do not show here. Solid lines are heavy quarks and antiquarks, dashed lines are temporal gluons, and curly lines are transverse gluons. Open and filled circles represent insertions of the p · A vertex and the spin-dependent σ · B vertex, respectively. Filled squares represent the operator χ † · σψ for spin triplet, and χ † ψ for spin singlet.
Since the logarithmic UV divergences in the vertex correction diagrams come from the contribution from the gluon pole, there are no contributions to the anomalous dimension at relative order α s v 2 that comes from gluon exchanges between the quark and the antiquark when the virtual quark or the virtual antiquark is on shell. Hence, at one loop level, there is no contribution to the anomalous dimension that comes from exchanges of potentials between the Q andQ. This is consistent with the pNRQCD expressions of the LDMEs in eqs. (2.4) and (2.8), where the one-loop anomalous dimension comes from the gluonic correlator E 3 , and not from corrections to the wavefunctions at the origin.

A.2 Two-loop anomalous dimension at leading order in v
Now let us consider the UV divergences in the two-loop corrections to the NRQCD LDMEs. Since we work at leading order in v, we can set the relative momentum q between the quark and antiquark to zero. The two-loop diagrams that contain logarithmic UV divergences are shown in fig. 6. We neglect diagrams that do not contain logarithmic UV divergences, most of which are scaleless power divergent, and hence vanish in DR. As we will see later, in the Coulomb gauge the diagrams only involve single poles in , and hence, we can set = 0 in the numerators of loop integrands without affecting the 1/ poles. The non-Abelian diagrams yield where we only keep the UV pole. The spin-independent ladder diagrams yield where again we only keep the UV pole. The spin-dependent ladder diagrams give where we use the notation ⊗ to make clear that the Pauli matrix on the left acts on the quark line, while the one on the right acts on the antiquark line. We first integrate over k, and average over the angles of to obtain the following UV-divergent contribution from the spin-dependent ladder diagrams: where again we keep only the UV pole. The spin-dependent factor σ i ⊗ σ i yields, for spin triplet, and for spin singlet, We collect the results in eqs. (A.8), (A.9), and (A.11) to obtain the logarithmic UV divergence in the two-loop correction to the LDME 0|χ † · σψ|QQ , which reads This reproduces the order-α 2 s term of the anomalous dimension in eq. (2.3). Similarly, the logarithmic UV divergence in the two-loop correction to the LDME 0|χ † ψ|QQ is which agrees with the order-α 2 s term of the anomalous dimension in eq. (2.7). We note that in the calculation of the two-loop diagrams in fig. 6, at least one of the integrations over the temporal components of the loop momenta must involve residues of the poles from the quark or antiquark propagators in order to produce logarithmic UV divergences. This is clear in the ladder diagrams, because the temporal-gluon propagator does not have a pole in the temporal components of loop momenta. For the non-Abelian diagrams, it can be shown that if we neglect the pole that comes from the transverse gluon propagator in the integration over 0 , we obtain the same UV pole as in eq. (A.8). This shows that the two-loop anomalous dimensions come solely from exchanges of potentials between the Q andQ. This is consistent with the pNRQCD expressions of the NRQCD LDMEs in eqs. (2.4) and (2.8), which imply that the two-loop anomalous dimension can only come from the wavefunctions at the origin.

B Potentials in perturbative QCD
In this appendix, we list the short-distance behaviors of the potentials, which are obtained from perturbative QCD. In perturbative QCD, the static potential is given through relative order α 2 s by [49,50] V (0) (r) pert = − α s (µ)C F r 1 + where α s = α s (µ) is the MS-renormalized QCD coupling constant at scale µ, and the functions a n (r; µ) are defined by a 1 (r; µ) = 31C A − 20T F n f 9 + 2β 0 log(µe γ E r), (B.2a) , γ E is the Euler-Mascheroni constant, n f is the number of light quark flavors, andā 1 = a 1 (r = e −γ E /µ; µ). We note that the dependence on µ cancels order by order in eq. (B.1).
The forms of the 1/m and 1/m 2 potentials generally depend on the matching scheme in which the potentials are determined. In on-shell matching, where we match on-shell S-matrix elements in NRQCD and pNRQCD in momentum space, we obtain [30,43,44,47,[77][78][79][80]  We use the superscript OS to denote the on-shell matching scheme. In Wilson-loop matching, the potentials at short distances are given by [45] V (1) (r) where m P is the mass of the quarkonium P , and c p and d p read, in the MS scheme [13,87,88], Again, Λ is the scale at which the NRQCD LDMEs are renormalized. To the best of the author's knowledge, the order-α s correction to d p has not been computed yet. Since the QCD renormalization scale dependence cancels order by order in c p , eq. (C.5a) is still valid if we replace α s (m) by α s (µ R ) and add − 3αsC F 2π × αsβ 0 4π log(µ 2 R /m 2 ), which compensates for the running of α s .
Finally, the two-photon decay rate of a pseudoscalar quarkonium P is given by where the SDCs c γγ and d γγ are given in the MS scheme by [4,12,83,[89][90][91][92][93][94] lbl + O(α 3 s ), (C.8a) Here, Λ is the scale at which the NRQCD LDMEs are renormalized. The constants f (2) reg and f (2) lbl have been determined numerically in ref. [92], which read where n H = 1, e q is the fractional charge of the light quark with flavor q, and the sum runs over n f light quark flavors. Since the QCD renormalization scale dependence cancels order by order in c p , eq. (C.8a) is still valid if we replace α s (m) by α s (µ R ) and add − αsC F π π 2 8 − 5 2 × αsβ 0 4π log(µ 2 R /m 2 ).

D Wavefunctions at the origin in perturbative QCD
If we ignore the nonperturbative long-distance behavior of the static potential, so that V LO (r) = −α s C F /r, the Schrödinger equation can be solved exactly, and the S-wave contribution to the Green's function in position space is known analytically: where r < = min(|r|, |r |), r > = max(|r|, |r |), λ = α s C F / −4E/m, and This result can be obtained by solving the differential equation in eq. (5.21) analytically. The bound states can be identified from the poles of Γ(−λ), which are located at λ = n with principal quantum numbers n = 1, 2, 3, . . .. The reduced Green's functions can also be obtained analytically from eq. (D.1). This makes possible analytical calculations of the corrections to the S-wave wavefunctions at the origin. Such a calculation has been done in the context of heavy quark pair production near threshold in perturbative QCD in refs. [25][26][27][28][29][30][31][32][33]35]. We use the known results in perturbative QCD to check the numerical procedure used in this paper for computing the divergent corrections to the wavefunctions at the origin that originate from the 1/m and 1/m 2 potentials. The explicit analytical expressions for the two-loop non-Coulombic corrections to the wavefunctions at the origin from the 1/m and 1/m 2 potentials can be found in ref. [28] for the spin-triplet state. Rather than comparing the wavefunctions at the origin, which depends on the renormalization scheme, it is simpler to compare the corrections to the decay constant f V , which is scale and scheme independent. At leading order in α s and v, f V for the spin-triplet nS state is given in perturbative QCD by where |Ψ LO n (0)| 2 = (α s C F m) 3 /(8πn 3 ). The corrections to the wavefunctions at the origin coming from the 1/m and 1/m 2 potentials are given in ref. [28] by where µ f is a factorization scale, and the ellipsis represent the Coulombic corrections that we neglect. The wavefunctions at the origin in ref. [28] are renormalized in a scheme that is different from the MS scheme that we use in this paper, so it is not possible to compare eq. (D.4) directly with the results in this paper. This scheme dependence cancels in the decay constant against the scheme dependence in the hard matching coefficient given by eq. (2) in ref. [28], which is obtained from the direct matching procedure [51]. Since the non-Coulombic corrections are proportional to α 2 s C 2 F and α 2 s C F C A , we only need to keep the contributions that are proportional to C 2 F and C F C A in the loop corrections to the hard matching coefficient. By combining the corrections to the wavefunctions at the origin and the hard matching coefficient, we obtain the two-loop non-Coulombic correction given by where the last two lines correspond to the C 2 F and the C F C A terms of the two-loop corrections to the hard matching coefficients in ref. [28]. The µ f dependence cancels exactly between the non-Coulombic corrections to the wavefunctions at the origin and the two-loop corrections to the hard matching coefficient. Note that the last two lines of eq. (D.5) differ from the C 2 F and the C F C A terms of the two-loop corrections to the SDC c v in eq. (C.2a). This reflects the difference between the renormalization scheme used in ref. [28] and the MS scheme used in this work. This analytical result can be compared with the numerical calculation in this paper, which is given by where δ NC Ψ pert is equal to eq. (5.29a), except that we set V LO (r) = −α s C F /r, V