Charm-quark mass effects in NRQCD matching coefficients and the leptonic decay of the $\Upsilon(1S)$ meson

We compute two-loop corrections to the vector current matching coefficient involving two heavy quark masses. The result is applied to the computation of the $\Upsilon(1S)$ decay width into an electron or muon pair. We complement the next-to-next-to-next-to-leading order corrections of Ref. arXiv:1401.3005 by charm quark mass effects up to second order in perturbation theory. Furthermore, we apply the formalism to $\Gamma(J/\Psi\to \ell^+\ell^-)$ and compare to the experimental data.


I. INTRODUCTION
Bottomonia, the bound states of a bottom and an antibottom quark, are excellent systems to investigate the dynamics of bound states in QCD. On the experimental side, there exist precise measurements of their properties. And on the theoretical side, the large mass of the bottom quark means that perturbation theory can be applied. This is in particular the case for the Υ(1S) meson. Nevertheless, its description is complicated by the fact that aside from the bottom-quark mass m b (the hard scale), there are two more relevant scales: the typical momentum and energy of the quarks, which are of order m b v (the soft scale) and m b v 2 (the ultrasoft scale), respectively. The Υ(1S) is a non-relativistic bound state, where the relative velocity v of the quark and antiquark is small, which means that these scales are well separated. It is then convenient to use an effective theory for the description of this multiscale problem. Starting from QCD, we first integrate out the hard modes to arrive at non-relativistic QCD (NRQCD). In a second step we integrate out the soft modes and potential gluons with ultrasoft energies and soft momenta to arrive at potential NRQCD (PNRQCD). At each step one has to determine the Wilson or matching coefficients of the corresponding effective theory, which are the couplings of the effective operators. For comprehensive reviews on this topic we refer to [2,3] The main focus of this paper is the matching coefficient c v of the vector current in NRQCD. Among other observables, it contributes to the decay rate of an Υ(1S) to a lepton-antilepton pair. In PNRQCD and to next-to-next-tonext-to-leading order (N 3 LO) accuracy, the decay rate is given by the formula [4] where α is the fine structure constant and m b the bottom-quark pole mass. E 1 and ψ 1 (0) are the binding energy and the wave function at the origin of the (bb) system. For convenience we provide the leading order results which are given by The matching coefficient c v of the leading current is known at the three-loop level [5][6][7] for the case of one massive quark and n l massless quarks. d v is the matching coefficient of the sub-leading bb current in NRQCD. Since it is multiplied by E 1 , it is only required at the one-loop level. This result can be found in Ref. [3]. Together with the N 3 LO results for the energy levels and the wave function at the origin [3,4,8], this made it possible to evaluate the decay rate at N 3 LO in Ref. [1]. One approximation that was made in Ref. [1] was to treat the charm quark as massless. The aim of this paper is to go beyond this approximation and include the corrections due to the charm-quark mass at next-to-next-to-leading order (NNLO). If we consider the charm-quark mass m c to be formally of the order of the hard scale m b , the charm quark has to be integrated out of QCD, leading to NRQCD with two heavy quarks with different masses. All NRQCD matching coefficients will then receive contributions due to m c . However, at NNLO only c v is affected. Thus we have to compute the fermionic contribution to the two-loop corrections to c v for a second non-zero quark mass. The analytic result for this contribution completes the two-loop evaluation of c v and together with its application to Γ(Υ(1S) → + − ) constitutes the main result of our paper.
Another possibility to include the charm-quark mass effects is to consider m c to be soft. In this case the charm quark is integrated out of NRQCD. Then there is no contribution to c v , but instead to the matching coefficients of PNRQCD, which are the potentials in the Schrödinger equation describing the (bb) system. At NNLO, only the Coulomb potential is affected (see Section 3.5 of Ref. [9]). Thus, the m c dependence then enters in the wave function and binding energy. We will compare the results of these two approaches.
The remainder of the paper is organized as follows: In the next section we describe the calculation of c v and in Section III the discussion is extended to external axial-vector, scalar and pseudo-scalar currents, where in addition to the non-singlet also the singlet contributions have to be considered. In Section IV we provide updated predictions for Γ(Υ(1S) → + − ) and in Section V we consider the decay of the J/Ψ and provide predictions of Γ(J/Ψ → + − ) up to N 3 LO. Our conclusions are presented in Section VI. In the Appendix analytic results for all matching coefficients up to two loops, which are not presented in the main part of the paper, are provided. The supplementary material to this paper [10] contains computer-readable expressions of all matching coefficients and all master integrals, which we compute in this paper.

II. TWO-LOOP MATCHING COEFFICIENT FOR THE VECTOR CURRENT WITH TWO MASSES
The matching coefficient for the vector current is defined via where m q is the heavy quark mass and j i v andj i v are currents defined in the full (QCD) and effective (NRQCD) theory. They are given by where φ and χ are two-component spinors. Note that in the heavy quark limit the 0 th component of j µ v is of order 1/m 2 q . A convenient approach to compute c v is based on the so-called threshold expansion [11,12] which is applied to the vertex corrections of a vector current and a heavy quark-antiquark pair, Γ v . Denoting by Z 2 the on-shell quark wave function renormalization constant one obtains the equation [6] Note that the vector current in QCD has a vanishing anomalous dimension whereasZ v deviates form 1 at order α 2 s . It gets contributions from the colour factors C 2 F and C A C F which are not considered in this paper. The momenta q 1 and q 2 in Eq. (5) correspond to the outgoing momenta of the quark and antiquark which are considered on-shell. Furthermore, we have (q 1 + q 2 ) 2 = 4m 2 q , a consequence of the threshold expansion. The quantity Γ v is conveniently obtained with the help of projectors applied to the vertex function Γ µ . It is straightforward to show that one gets with where q = q 1 + q 2 .
In Fig. 1 we show sample diagrams contributing to c v up to two-loop order. The main focus of this work is the Feynman diagram in Fig. 1(a) where the quark in the closed loop has mass m 2 . For the computation of this diagram we proceed as follows. • In order to fix the boundary conditions we consider the limits x → 0 and x → 1. This is necessary since in each individual limit some of the integration constants drop out. Alternatively it would be possible to solve the differential equation to higher order in . The values of I for x = 0 and x = 1 are used to determine the integration constants in J.
We want to remark that we use a general QCD gauge parameter ξ for our computation. The independence of our final result from ξ is a welcome cross check. Our results for the master integrals agree with those given in Appendix B of Ref. [17].
For the renormalization of our two-loop contribution we need the two-loop corrections to the on-shell wave-function renormalization constant. The analytic results for the contribution involving m 2 can be found in Refs. [17][18][19][20][21]. Furthermore the one-loop counterterm for α s is needed. At this point we use Eq. (5) in order to extract c v .
For the fermionic contributions to c v it is convenient to introduce n f = n h + n m + n l , where n h = 1 and n m = 1 label the contributions with closed massive quark loops with mass m q and m 2 , respectively. n l counts the massless quarks. Using this notation we can cast the result for c v in the form where α (n l +nm) s (µ) is the strong coupling constant where the heavy quark with mass m q is decoupled from the running of α s . The m 2 -independent contributions to c (2) v can be found in Refs. [5,6,22]. The new contribution proportional to n m reads where H i... = H i... (x) and H i... (x) are harmonic polylogarithms (HPLs) [16]. Note that for x → 1 we reproduce the known result for m 2 = m q which is given by However, in the limit x → 0 we do not obtain the massless fermion contribution but recover the well-known Coulomb singularity, which is regulated by the mass m 2 . For small m 2 we have In the application we discuss in Section IV we need c v expressed in term of α (nl=3) s which means that we have to decouple the charm quark from the running of α s . As a consequence µ 2 is effectively replaced by m 2 2 in Eqs. (14), (15) and (16) disappear.
Let us finally investigate the numerical effect of the new contribution. We specify to the bottom-charm system and use m 2 = m c = 1.65 GeV and m q = m b = 5.1 GeV for the pole masses of the charm and bottom quarks. This leads to where the contributions originating from the closed massless, bottom and charm quark loops are marked by n l = 3, n h = 1 and n m = 1, respectively. One observes that the coefficient of n m is more than a factor four times larger than the coefficients of n l and n h . Thus, the contribution of the heavy quark with mass m 2 is larger than the contributions of the heavy quark with mass m q and all three massless quarks combined.
III. TWO-LOOP TWO-MASS MATCHING COEFFICIENTS FOR AXIAL-VECTOR, SCALAR AND PSEUDO-SCALAR CURRENTS.
In this Section we consider further external currents, which are of phenomenological relevance, and compute the corresponding matching coefficients between QCD and NRQCD to two-loop order. Such currents have, in contrast to the vector case, both non-singlet and singlet contributions. The latter are characterized by the fact that the external current does not directly couple to the quarks in the final state but only through the exchange of two gluons. A sample Feynman diagram is shown in Fig. 1 We write the two-loop corrections in the form x,non−sing + c (2) x,sing , where x ∈ {a, s, p} stands for an axial-vector, scalar or pseudo-scalar. All two-loop corrections which involve only one mass scale have been computed in Ref. [22]. In this paper we concentrate on the diagrams where a second massive quark in a closed loop is present which concerns both the non-singlet and the singlet contribution.
In analogy to Eq. (3) we define the additional currents in QCD via The anomalous dimension of j µ a is zero. For the scalar and pseudo-scalar current we have for the corresponding renormalization constant Z s = Z p = Z m , where Z m is the (on-shell) mass renormalization constant.
In NRQCD the currents read [22] where k = 1, 2, 3. Furthermore we have j 0 a = ij p which constitutes an alternative way to compute the matching coefficient c p . Note the presence of the momentum p, which is the relative momentum of the external quark and antiquark, in the definition of the axial-vector and scalar current. Thus, an expansion in p has to be performed in order to obtain the loop corrections to the corresponding matching coefficients.
The matching equation in (5) also holds for the other currents after the obvious replacements of Γ v ,Z v andΓ v . For the pseudo-scalar current and the zero component of the axial-vector the momentum p is zero and the calculation proceeds in close analogy to the vector case. In fact, we have [22]. with For the axial-vector and scalar cases there are similar equations to Eq. (21). The expansion in p (up to linear order) is conveniently realized by choosing q 1 = q/2 + p and q 2 = q/2 − p, which implies q · p = 0. Thus the projectors are more complicated and are given by After the application of the projectors and the expansion in p we can set p = 0 and q 2 = 4m 2 q . The calculation of the non-singlet contribution is in close analogy to the vector case, see Section II. In particular, it is possible to use anticommuting γ 5 . Furthermore, we can map the scalar integrals contributing to Γ x after the application of the projector to the same integral families and we thus end up with the same master integrals.
The singlet contribution is more involved and a few comments are in order. Let us first mention that a nonzero contribution for the scalar and pseudo-scalar currents is only obtained for massive quarks in the closed fermion loop. Furthermore, for the axial-vector current an effective current formed by the difference of the upper and lower component of a given quark doublet should be considered in order to guarantee the cancellation of anomaly-like contributions. For example, for the top-bottom doublet we have In practice, this means that we have a quark with mass m q in the final state and we consider both a massless quark and a quark with mass m 2 in the closed quark loop and take the difference.
In the singlet diagrams we treat γ 5 according to the prescription of Ref. [23]. In the Feynman diagrams we apply for the axial-vector and pseudo-scalar couplings the replacements We perform the same substitution also in the corresponding projectors. Afterwards we strip off the two ε tensors and interpret the product in d dimensions. This allows us to perform the calculation in close analogy to the scalar current. The remaining calculation of the singlet diagrams proceeds as outlined in the previous section. After applying a partial fraction decomposition, we can map all integrals in our amplitude to two integral families which are given by The reduction to master integrals using FIRE [13] and LiteRed [14] leads to 12 master integrals which are shown in Fig. 3. In a next step we establish differential equations in the variable t defined by x = 2t/(1 + t 2 ). In this new variable the differential equation can be brought into -form with the help of CANONICA [15]. We expand the solution including terms of order since some of the master integrals have 1/ poles in the prefactor. The differential equations are integrated with the help of HarmonicSums [24] in terms of cyclotomic harmonic polylogarithms over the alphabet Alternatively one could factorize the denominators over the complex numbers and arrive at Goncharov polylogarithms. The boundary conditions are obtained from the single-scale master integrals needed for the two-loop calculation of Refs. [22,25]. For the master integrals with dots the naive m 2 = 0 limit is not enough and we have to consider the asymptotic expansion around m 2 = 0 which can be obtained easily from one-dimensional Mellin-Barnes representations or a diagrammatic large momentum expansion of the corresponding Feynman integrals. In a second approach we use the algorithm described in [26] to solve the differential equations in the variable x without going into an -form first.
For the implementation we additionally make use of Sigma [27] and OreSys [28]. This approach introduces the squareroot valued letter √ 1 − τ 2 /τ . Both results agree after the above mentioned variable transformation. We compute the expansion of all master integrals up to the order which is needed to obtain the O( ) terms of the matching coefficients.
After inserting the master integrals into the integration-by-parts-reduced amplitude we obtain for the two-loop singlet contribution to the matching coefficient of the scalar current the following expression c (2) s,sing m2 = n m C F T F 4t 3 1 + t 2 + π 2 t 1 − 28t 2 − 35t 4 18 1 + t 2 3 with t = (1 − √ 1 − x 2 )/x and H a = H a (t). The imaginary part of the matching coefficient is displayed in the last three lines of Eq. (29). For the expansions around x = 0 we find c (2) s,sing m2,x→0 As expected, this contribution to the matching coefficient is zero for vanishing quark mass in the closed triangle. Note that mass corrections are linear in m 2 . On the other hand, for x → 1 c The expressions for the pseudo-scalar and axial-vector currents can be found in Appendix A, where we also show the non-singlet terms.

IV. Γ(Υ(1S) → + − ) AND FINITE CHARM QUARK MASS
In Ref. [1] the charm quark has been treated as massless and the decay rate has been expressed in terms of α (n l ) s (µ) with n l = 4. In the following we discuss the additional ingredients needed for the finite charm quark mass terms. As mentioned in the Introduction we consider two scenarios: A. m c is hard and the charm quark is integrated out when matching QCD to NRQCD. In this approach we express Γ(Υ(1S) → + − ) in terms of α B. m c is soft and thus the charm quark is a dynamical scale within NRQCD. We express Γ(Υ(1S) → + − ) in terms of α (4) s (µ). In this approach charm mass effects to bound-state energies and wave functions are needed. They are known at NLO [29] and NNLO [9,30]. We use the expressions given in Ref. [9].
In case the decay rate shall be expressed in terms of the potential subtracted mass the charm quark mass effects are needed to NNLO [9].
All necessary expressions for this scenario are available in the program QQbar threshold [31].
In both scenarios charm mass effects to the relation between the MS (which we use as input) and on-shell bottom quark mass are taken into account. They are known to three-loop order [19,21].
In scenario A we assume that m c is parametrically of the order of m b . In such a situation both m b and m c have to be decoupled from the running of α s and α (3) s is used as an expansion parameter. In fact it has been observed (see, e.g., Ref. [32]) that, e.g., the finite-m c terms to the MS-on-shell relation of the bottom quark are quite sizeable and do not converge in case α (4) s is used as parameter. On the other hand, charm quark mass corrections are small and well convergent for α (3) s . To arrive at the new result for Γ(Υ(1S) → + − ) we proceed as follows. Our starting point is the expression derived in Ref. [1] where α (4) s has been used as expansion parameter. For the number of massless quarks we have n l = 4. We restore the dependence on (massless) charm quarks and write n l = n l + n m with n l = 3 and n m = 1. In scenario B we can simply add the finite-m c terms from the binding energy and wave function. This modifies the coefficient of n m such that in the limit m c → 0 the coefficient of n l is recovered.
In scenario A we interpret the result of Ref. [1] in the n l = 3-flavour PNRQCD with an expansion parameter α s . Finite-m c effects enter in Eq. (1) only via the matching coefficient c v (cf. Section II) which also has to be expressed in terms of α (3) s . We are now in the position to provide numerical results for the decay rate. For the numerical evaluation we use α(2m b ) = 1/132.3 [33], α (5) s (M Z ) = 0.1179(10) [34] and the renormalization scale µ = 3.5 GeV. We use the program RunDec [35] to evolve the coupling with five-loop accuracy and obtain α where the four terms in the first lines of the two equations refer to the LO, NLO, NNLO and N 3 LO results. At NNLO and in scenario B also at NLO we display the contributions from a finite charm quark mass separately. We remark that the finite-m c terms of c (2) v , which are computed in Section II, amount to about 2% of the NNLO coefficient and they are of the same order of magnitude as the N 3 LO contribution. In scenario B the m c effects at NLO and NNLO are of the same order of magnitude and amount to about 15% and 5% of the corresponding m c -independent coefficient. The scale uncertainty in the last line of Eq. (32) is computed from the variation of µ in the range µ ∈ [3, 10] GeV. We  It is interesting to note that both scenario A and scenario B lead to the same final prediction at N 3 LO although the contributions from the various orders is different. We oberserve a notable discrepancy to the experimental result which is given by Γ(Υ(1S) → e + e − )| exp = 1.340 (18) keV. We also want to mention a recent lattice evaluation [37] where the value Γ(Υ(1S) → e + e − ) = 1.292(37)(3) keV has been reported.
In Fig. 4 we show the dependence of Γ(Υ(1S) → e + e − )| pole,A on µ successively including higher order corrections. The solid black line corresponds to the N 3 LO prediction. We observe that the inclusion of higher order corrections clearly stabilizes the perturbative predictions for µ ∼ > 3 GeV. Furthermore, it is interesting to note that the thirdorder corrections vanishes close to the value of µ where the N 3 LO curve has a maximum. The dashed black curve corresponds to the N 3 LO prediction of scenario B. The overall shape is very similar to the corresponding curve of scenario A. However, it is remarkable that the two N 3 LO lines cross the NNLO curve for the same value of µ.
In Fig. 5 we show Γ(Υ(1S) → e + e − )| pole,A as a function of α (5) s (M Z ). One observes that the third-order band is embedded by the NNLO band which can be interpreted as good convergence of the perturbative corrections. Note that we do not recompute the bottom pole mass when varying α s .
It is well-known that the pole mass suffers from so-called renormalon ambiguities. They are avoided by choosing a properly defined so-called threshold mass. Such masses have the advantages that they have nice convergence properties (as the MS mass) and that they can also be used for the description of bound-state properties. In the following we want to consider the potential-subtracted (PS) mass scheme [38] as an example and discuss the perturbative corrections to Γ(Υ(1S) → + − ).
Explicit results for the relation between the pole mass and the PS mass to n-th order can be derived from the n-loop expression for the Coulomb potential (see for example Ref. [39]). For scenario A, we use this relation for n = 3 and n l = 3, since in this scenario finite charm-quark mass effects are only included in the relation between the MS mass and pole mass. In scenario B, however, we also have to include charm-mass effects in the relation between the pole mass and PS mass for n = 1 and n = 2. The latter can be found in Appendix B of Ref. [9].
The numerical (input) value for the PS mass is conveniently obtained fromm b (m b ) = 4.163 GeV. Using N 3 LO accuracy we obtain for the two scenarios The final predictions for the decay rate are close to those in the on-shell scheme (cf. Eqs. (32)) and agree well within the uncertainties. However, the transition from the pole to the PS mass leads to a significant redistribution among the various perturbative orders. For example, in scenario A the N 3 LO term in the PS scheme is about four times larger as compared to the on-shell scheme but has a different sign. Similarly, in scenario B the N 3 LO coefficient gets reduced by a factor six.  The inclusion of the finite-m c effects leads to the same conclusions as in Ref. [1]: The perturbative predictions for Γ(Υ(1S) → + − ) are well under control but there is a discrepancy with respect to the experimental result. In [1] one can find an extensive discussion on possible non-perturbative effects. However, no clear conclusion can be drawn and it remains an open question whether a full quantitative understanding of the decay rate based on perturbative and non-perturbative QCD is possible.   In this section we apply the formalism of Ref. [1] to the decay of the J/Ψ to massless leptons. In general, the application to charm bound states is questionable, since the ultra-soft scale in PNRQCD is smaller then Λ QCD . Furthermore, even the hard scale (m c ) is below 2 GeV. Nevertheless, it is interesting to study the perturbative behaviour and to compare with the experimental result.
We work with an on-shell charm quark mass value m c = 1.65 GeV and choose µ = 2 GeV for the renormalization scale. This leads to α where the scale uncertainty is computed from the variation of µ in the range µ ∈ [1.5, 6] GeV. Although the perturbative series does not converge it is instructive to compare to the experimental result. This is done in Fig. 8 where Γ(J/Ψ → + − ) is shown as a function of µ. It is interesting to notice that there is agreement between the N 3 LO prediction and the experimental result Γ(J/Ψ → + − )| exp = 5.53 ± 0.10 keV [34] close to the value of µ where the N 3 LO curve has a maximum and thus the derivative with respect to µ vanishes. Furthermore, for this value of µ the third-order corrections are are quite small, as can be seen from Eq. (34). Note that the N 3 LO corrections vanish for µ = 1.724 GeV. For this value of the renormalization scale we have Γ(J/Ψ → + − )| pole = 5.03 keV. From Fig. 8 we also observe that even the N 3 LO curve shows a sizable dependence on µ. Furthermore, one notices that below µ ≈ 1.5 GeV ≈ m c perturbation theory breaks down. We want to remark that a similar feature has been observed in Ref. [40] where next-to-leading logarithmic (NLL) corrections to the hyperfine splitting of heavy quark-antiquark bound states have been considered. The application to the 1S charmonium states shows good agreement for values of the renormalization scale where the NNL prediction has a maximum. The perturbative uncertainties are sizeable, as for the J/Ψ decay rate.

VI. CONCLUSIONS
In this paper we consider the matching coefficients between QCD and NRQCD of external vector, axial-vector, scalar and pseudo-scalar currents. We compute all two-loop contributions which involve two mass scales, one from the external quarks and one present in a closed fermion loop. Whereas for the vector current only non-singlet contributions have to be considered there are also singlet contributions for the other three currents. We present analytic results including terms of order , which are of relevance for a future three-loop calculation.
In Sections IV and V we apply our results for the vector current to the leptonic decay rates of the lowest spin-1 heavy-quark-anti-quark mesons, Υ(1S) and J/Ψ, and provide update numerical predictions. We discuss the decay rate Γ(Υ(1S) → + − ), including charm quark mass effects, both in the three-and four-flavour scheme and for the heavy quark masses defined both in the on-shell and PS scheme. Although there is a good convergence of the perturbative corrections, we observe a discrepancy with respect to the experimental result which to date is not understood.
where H a = H a (x).
In the case of the singlet axial-vector current we explicitly specify the flavour of the quark in the final state to bottom quark. Furthermore, we split the matching coefficient into the contributions from the strange and charm quarks c = n m C F T F π 2 − t 2 1 + 25t 2 + 19t 4 +  = n m C F T F x 2 c −