Theory for quarkonium: from NRQCD factorization to soft gluon factorization

The widely used nonrelativistic QCD (NRQCD) factorization [1] has encountered notable difficulties in describing heavy quarkonium data. As the NRQCD factorization is based on the NRQCD effective field theory [2], it is likely rigorous, although for inclusive quarkonium production, only two-loop verification is available at present [3-5]. The main known problem of NRQCD factorization is the bad convergence of its relativistic expansion [6], which may be responsible for its deficiencies. Recently, a new factorization approach called soft gluon factorization (SGF) [7,8] was proposed to describe quarkonium production and decay. The aim of SGF is to resum the series of relativistic corrections originating from the kinematic effects in NRQCD, which are the main cause of the poor convergence of the relativistic expansion.

Nevertheless, SGF has not been well-established. In SGF, the hadronization of intermediate quark-antiquark pairs to physical quarkonium is described by nonperturbative soft gluon distributions (SGDs), which are only formally defined by QCD fields in a small loop momentum region [7]. Without an explicit definition of a small region, it is hard to prove the validity of SGF for physical processes. Furthermore, the unclear relation between SGF and NRQCD factorization makes it impossible to verify whether the kinematic effects have been correctly resummed.

In this paper, with the help of a new regulator, we provide a rigorous definition of a small region in SGF. We then present two strategies for exploring the relationship between SGF and NRQCD factorization. In the first strategy, we apply the two factorization theories to the physical process of $ J/\psi \to e^+e^- $ and show that both SGF and NRQCD factorization can correctly reproduce all the low energy physics of full QCD in this process. In the second strategy, we argue that the SGF formula can be deduced from NRQCD at the operator level by using equations of motion. Both strategies demonstrate that SGF and NRQCD factorization are equivalent, which means that, for any process, the two factorizations theories are either both valid or both violated. By identifying the two theories, we generate complete relations between the nonperturbative matrix elements in SGF and NRQCD; this proves the generalized Gremm-Kapustin relation [9] as a byproduct.

The rest of this paper is organized as follows. In Sec. II, we study the exclusive process $ J/\psi \to e^+e^- $ in SGF. We give a rigorous definition of nonperturbative matrix elements in SGF and show that the low energy physics of full QCD can be reproduced by SGF. In Sec. III, we discuss the equivalence between SGF and NRQCD factorization and establish relations between the nonperturbative matrix elements in SGF and NRQCD. In Sec. IV, we show that SGF can be deduced from NRQCD at the operator level. We summarize our results in Sec. V. Some of the technical details of our calculations are provided in Appendix A.

According to Ref. [8], for the exclusive decay process, we have the following SGF formula at the amplitude level:

where $ n $ denotes intermediate states. In general, $ n $ can contain dynamical soft partons (gluons or light quarks) in addition to a $ Q\bar Q $ pair. However, for simplicity, in this work we only discuss intermediate states without dynamical soft partons, which can nevertheless be treated similarly. The nonperturbative matrix elements $ \overline{R}^{n*}_{{ {\cal Q}}} $ are then defined as

where $ \Psi $ stands for the Dirac field of a heavy quark, $ {\cal{K}}_{n} $ is the projection operator defining the intermediate state $ n $, and the subscript "S" indicates that, to evaluate the matrix elements, one only includes integration regions where the off-shellness of all particles is much smaller than the heavy-quark mass. From the point view of the method of regions [10,11], the effect of "S" is to keep only small regions, which include everything except the hard regions.

For the process $ J/\psi \to e^+e^- $, the symmetries of QCD tell us that only $ n = {^{{3}}{S}_{{1}}^{[{1}]}}, {^{{3}}{D}_{{1}}^{[{1}]}} $ are relevant, where the spectroscopic notation, with the superscript "[1]", denotes a color singlet. We thus have

with the projection operators defined explicitly as

where $ d = 4-2\epsilon $ is the space-time dimension, $ D_\mu $ is the gauge covariant derivative with $ \overline\Psi \overleftrightarrow{D}_\mu \Psi = \overline\Psi (D_\mu \Psi) - $$(D_\mu \overline\Psi)\Psi $, $ P $ is the momentum of $ J/\psi $, $ {\epsilon}_{S_z}^{\mu} $ is a polarization vector with $ P\cdot {\epsilon}_{S_z} = 0 $, $ m $ is the heavy-quark mass, $ M $ is the mass of $ J/\psi $, the color projector is defined as $ {\cal C}^{[1]} = 1/\sqrt{N_c} $, and the spin projection operator $ \mathbb{P}_{\alpha\beta} $ is defined as

The hard parts $ \hat{{\cal{A}}}^{n} $ can be perturbatively calculated according to the matching procedure discussed in [7,8]. To this end, we first replace $ J/\psi $ with an on-shell color-singlet state $ c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}}) $ with momenta^①

on both sides of Eq. (3). The on-shell conditions $ p_c^2 = p_{\overline{c}}^2 = m^2 $ result in

which fixes $ q_0 $ and $ |{{q}}| $ in the rest frame of $ P $. The remaining degrees of freedom of $ q $ are removed by partial wave expansion, i.e., $ S $-wave expansion in this case. After the replacement, the l.h.s. of Eq. (3) becomes

where $ \Omega $ is the solid angle of the relative momentum $ {{q}} $ in the $ c\bar{c} $ rest frame, $ L_\mu $ is the leptonic current with

and $ {\cal{A}}^{\mu }_{c\bar c} $ is the hadronic part of the decay amplitude with the spinors of $ c\bar{c} $ removed. The $ c\bar{c} $ pair is projected to state $ {^{{3}}{S}_{{1}}^{[{1}]}} $ by replacing the spinors of the $ c\bar{c} $ pair by

Similarly, we have

Based on these equations, one can calculate $ {\cal{A}}^{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})\to e^+e^-} $, $ \overline{R}^{ {^{{3}}{S}_{{1}}^{[{1}]}} *}_{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})} $ , and $ \overline{R}^{ {^{{3}}{D}_{{1}}^{[{1}]}} *}_{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})} $ perturbatively.

Denoting the perturbative expansion of any quantity $ W $ as $ W = W^{(0)}+\alpha_s W^{(1)}+\alpha_s^2 W^{(2)}+\cdots $, we have the following orthogonal relations [7]:

Based on this, Eq. (3) results in the following matching relations:

By replacing $ J/\psi $ with $ c\bar{c}( {^{{3}}{D}_{{1}}^{[{1}]}}) $, we can obtain similar relations for $ \hat{{\cal{A}}}^{ {^{{3}}{D}_{{1}}^{[{1}]}}} $. These relations enable us to calculate $ \hat{{\cal{A}}}^{ {^{{3}}{S}_{{1}}^{[{1}]}}} $ and $ \hat{{\cal{A}}}^{ {^{{3}}{D}_{{1}}^{[{1}]}}} $ perturbatively. For simplicity, we concentrate on the S-wave contribution $ \hat{{\cal{A}}}^{ {^{{3}}{S}_{{1}}^{[{1}]}}} $ in the rest of the paper.

We first calculate $ {\cal{A}}^{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}}) \to e^+e^-} $ according to Eq. (8). The amplitude $ {\cal{A}}^\mu_{c\bar c} $ in full QCD can be decomposed as

and up to order $ \alpha_s $, one has [9]

with

where $ {\rm{Li}}_2 $ is the Spence function:

In the above results, we have dropped imaginary parts that are irrelevant for our purpose. By inserting Eqs. (14) and (10) into Eq. (8), one gets

Note that in Ref. [9], $ M $ is a free parameter, but to use these expressions for SGF, it needs to be the quarkonium mass.

We now describe our method for calculating $ \overline{R}^{ {^{{3}}{S}_{{1}}^{[{1}]}} *}_{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})} $. As pointed out in Refs. [7,8], this quantity has been defined to include only small loop momentum regions. In what follows, we provide an explicit definition and choose a UV renormalization scheme.

Up to order $ \alpha_s $, the corresponding Feynman diagrams are shown in Fig. 1, where the solid circle represents the operator $ \overline\Psi {\cal{K}}_{ {^{{3}}{S}_{{1}}^{[{1}]}}}\Psi $. The calculation at the tree level is straightforward; the result is

Figure 1.  Feynman diagrams for the $ {^{{3}}{S}_{{1}}^{[{1}]}} $ matrix element

Consider the vertex correction in Fig. 1(d) as an example to explain the calculation at the one-loop level. The amplitude reads

where $ T_S $ is an operator that forces the loop momentum $ k $ to be in a small region [7,8], which will be defined explicitly by using the method of regions [10,11].

Before continuing, we note that, although the full QCD integral in Eq. (20) is well regularized by dimensional regularization, the manipulation to define $ T_S $ generates unregularized integrals. Therefore, another regulator is needed to make our manipulation mathematically rigorous. We thus propose a new regularization at the full QCD level by multiplying the power of all propagator denominators by $ 1+\eta $; this can regularize all possible divergences encountered in the derivation using the method of regions, including both ultraviolet and non-ultraviolet divergences^②. A similar regularization method has been used in the light-cone ordered perturbation theory in Ref. [12] to regularize rapidity divergences. For the integral of interest, we obtain

where $ \nu $ is introduced to compensate for the change in mass dimension caused by the new regularization, and we assume $ \eta\ll \epsilon $ to ensure that the theory can eventually be regularized by dimensional regularization.

With the new regularization at hand, we decompose the loop momentum $ k_\mu = (k_0,{{k}}) $ into three domains:

where the relation "$ \lesssim $" is understood as the negation of "$ \gg $," $ D = \mathbb{R}^{d} $ is the complete integration domain, and an implicit cutoff scale exists to rigorously separate the three domains. This division satisfies

Then, one can split any original integral $ F $ into the three domains:

A possible definition of $ T_S $ is $ \{k\in D_s\cup D_p\} $. However, this definition involves a hard cutoff, which makes high order perturbative calculations inconvenient.

To obtain a more convenient definition, we introduce the operators $ T^{(i)}\; (i \in \{h,s,p \}) $ , which expand the integ

rand to a convergent power series of small quantities in each domain. We also define $ T^{(i,j,\cdots)}\equiv T^{(i)}T^{(j,\cdots)} $. We then have

where the property $ T^{(i,j)} = T^{(j,i)} $ has been used [11]. It is clear that the l.h.s. and the first term on the r.h.s. of this equation are equivalent in the low energy domain and that their difference is an integration in the hard domain, which is infrared safe but may be ultraviolet divergent. The difference can be interpreted as a different choice of renormalization scheme. Therefore, we arrive at our final definition:

with the UV divergences removed by an $ \overline{ {{\rm{MS}}}} $ renormalization scheme. $ T^{(s)} $ and $ T^{(p)} $ are conventionally respectively called the soft region and potential region, while the overlap region $ T^{(s,p)} $ compensates for the double counting between $ T^{(s)} $ and $ T^{(p)} $.

For the soft region contribution $ \overline{R}^{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})*,{(s)}}_{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})} $, we apply the operator $ T^{(s)} $ to expand the integrand in Eq. (21) for small quantities. This results in the following integrals:

where the term $ k_0^{m_1} ( {{k}}\cdot {{q}})^{m_2} $ comes from the numerator in Eq. (21), and

where $ j_i $ are non-negative integers. As scaleless and infrared-finite integrals can be set to zero (which is another choice of scheme), we find that only integrals with $ m_1 = m_2 = j_1 = j_2 = j_4 = j_5 = 0 $ are relevant:

These integrals have pinch poles around $ k_0 = 0 $ that can be regularized by the new regulator $ \eta $. As integrals in $ \overline{R}^{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})*,{(s,p)}}_{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})} $ can be obtained by expanding the $ k_0^2 $ term in the denominator in Eq. (29), it is clear that they cancel exactly with the contribution of the pinch poles around $ k_0 = 0 $ in $ \overline{R}^{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})*,{(s)}}_{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})} $. This cancellation demonstrates that the overlap region is conceptually important even though it contains only scaleless integrals that are usually set to zero in dimensional regularization. In fact, such overlap contributions have also been introduced as "zero-bin subtraction" in effective-theory treatments [13,14]; this resolves a number of problems in effective theories (NRQCD and soft-collinear effective theory) [13], including the pinch singularity problem discussed above. By combining the soft region and the overlap region, we only need to consider the contribution from the gluon pole, and we can safely take the limit $ \eta\to 0 $ , which results in

We emphasize that, if one wants to calculate $ \overline{R}^{ {^{{3}}{S}_{{1}}^{[{1}]}}*,(s)}_{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})} $ separately, the correct order is to first integrate over $ k_0 $ with fixed $ {{k}} $. Otherwise, the poles around $ k_0 = 0 $ will be regularized by different regulators between the soft region and overlap region, which necessarily breaks the symmetries of the theory and makes the cancellation between the two regions impossible.

For $ \overline{R}^{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})*,{(p)}}_{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})} $, the terms in the expansion are proportional to

where we can take $ \eta\to 0 $ as all divergences are well regularized by dimensional regularization. Then because $ k_0 $ and $ {{k}}\cdot {{q}} $ can be expressed as linear combinations of the denominators, if any of $ m_1, m_2, j_1, j_2 $, or $ j_3 $ is nonzero, the integral can be decomposed to either scaleless and infrared-finite integrals or to integrals with pure virtual value. By keeping only the real part, we get

Similarly, for the self-energy diagrams, one can derive

Summing these contributions, we obtain the matrix element at NLO before renormalization:

Ultraviolet divergences in the above result can be removed by the $ \overline{ {{\rm{MS}}}} $ renormalization procedure, yielding the renormalized matrix element

Similarly, because we find that the real part of $ \overline{R}^{ {^{{3}}{D}_{{1}}^{[{1}]}} *,(1)}_{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})} $ is proportional to $ \overline{R}^{ {^{{3}}{D}_{{1}}^{[{1}]}} *,(0)}_{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})} $, we have

Substituting Eqs. (18), (35), and (36) into Eq. (13), we obtain

where

We find that the matrix element defined in SGF reproduces all infrared and Coulomb divergences in full QCD and that the obtained hard part is free of divergences. We therefore conclude that the SGF factorization holds at least at the one-loop level.

The correctness of SGF at the one-loop order can be understood in the following way. The full QCD results in Eq. (18) can also be reproduced by the method of regions (see Appendix A for details), in which $ {\cal{A}}^{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})\to e^+e^-} $ is expressed as

Here, the first term on the r.h.s has nonzero support only in the hard domain and is thus infrared safe. It is straightforward to check that the low energy part of $ {\cal{A}}^{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})\to e^+e^-} $ has been correctly reproduced by the matrix element in SGF, i.e.,

which leaves the corresponding short-distance hard part defined by the high energy part of $ {\cal{A}}^{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})\to e^+e^-} $, and thus, it is infrared safe. More precisely, we have

where $ \overline{ {\rm{MS}}} $ denotes the removal of UV divergences by the $ \overline{ {\rm{MS}}} $ subtraction scheme^③.

The above one-loop argument can be generalized to all orders. By definition, the low energy part (the "small" region of loop momenta) in full QCD can be reproduced by the matrix element in SGF at any order in the $ \alpha_s $ expansion, with a proper definition of $ T_S $ at the multi-loop level as discussed in Ref. [11]. Therefore, the short-distance hard part is perturbatively infrared-safe, which means that the SGF formula for the decay width of $ J/\psi \to e^+e^- $ holds for all orders in $ \alpha_s $.

Ignoring operators involving gauge fields, the NRQCD factorization for the decay amplitude $ J/\psi \to e^+e^- $ is given by [9,15,16]

with^④

where $ \psi $ and $ \chi^\dagger $ are the two-component heavy quark fields in NRQCD, $ \overleftrightarrow{\nabla} $ is defined as $\chi^\dagger \overleftrightarrow{\nabla} \psi = \chi^\dagger (\nabla \psi) - $$ (\nabla\chi^\dagger)\psi $, and $ s_n( {^{{3}}{S}_{{1}}^{[{1}]}}) $ and $ b_n( {^{{3}}{D}_{{1}}^{[{1}]}}) $ are short-distance hard parts that can be perturbatively calculated. The first two orders in the $ \alpha_s $ expansion for $ s_n( {^{{3}}{S}_{{1}}^{[{1}]}}) $ are given by [9,15,16]

with $ {\cal{G}}^\prime $, $ {\cal{H}}^\prime $ given in Eq. (38).

In Ref. [9], the S-wave contributions in Eqs. (42) and (44) are further resummed by using the generalized Gremm-Kapustin relation [9,17,18]

with

These relations are obtained by computing the matrix element $ \langle 0 \vert{\cal{O}}_{An} \vert J/\psi \rangle $ in the potential-model [17].

The basic justification for the validity of SGF in this process is that the SGF matrix elements reproduce the low energy part of full QCD. As the NRQCD matrix elements also correctly reproduce the low energy part of full QCD, SGF is equivalent to NRQCD factorization. This means that for any process, the two factorization formulas are either both valid or both broken.

Because the amplitude of $ J/\psi \to e^+e^- $ can be factorized in both SGF and NRQCD factorization, we have ($ D $-wave contributions are suppressed)

which can generate relations between the SGF matrix elements and NRQCD matrix elements, with $ {\cal{O}}(v^2) $ denoting contributions from operators with explicit gauge fields. In fact, we can generate more relations by applying the two factorization formulas to any well defined QCD quantity $ W $:

For example, if we choose $ W = \frac{1}{4}(M+2m) (M-2m)$$\overline{R}^{ {^{{3}}{S}_{{1}}^{[{1}]}} *}_{J/\psi} $, we have $ \hat{W} = \frac{1}{4}(M+2m)(M-2m) $. To determine the corresponding $ w_n $, we replace $ J/\psi $ by a $ c\bar c $ pair with invariant mass $ M $^⑤. Using the nonrelativistic expansion formulas given in [19], we have

Here, we use nonrelativistic normalization for the spinors $ u $ and $ v $and thus obtain

where the extra factor $ \sqrt{2M} $ in the first line appears because the NRQCD matrix elements have nonrelativistic normalization, whereas SGF matrix elements have relativistic normalization. As both the SGF matrix elements and NRQCD matrix elements maintain the same low energy physics and are renormalized in the same way, the coefficients $ w_n( {^{{3}}{S}_{{1}}^{[{1}]}}) $ should vanish at higher orders in $ \alpha_s $, i.e.,

Thus, we get the relation

Similarly, by choosing $ W = [\frac{1}{4}(M+2m)(M-2m)]^n \times \overline{R}^{ {^{{3}}{S}_{{1}}^{[{1}]}} *}_{J/\psi} $, we can obtain

which provides complete relations between the SGF matrix elements and NRQCD matrix elements. Using these relations, we obtain

which agrees with Eq. (45). We have thus proved the generalized Gremm-Kapustin relation by using the equivalence between SGF and NRQCD factorization. Based on our proof, it is clear that the $ {\cal{O}}(v^2) $ terms in the relations are contributions from operators with explicit gauge fields.

We see that the equivalence between SGF and NRQCD factorization gives rise to the generalized Gremm-Kapustin relation. Moreover, by assuming these relations, the introduction of nonperturbative quantities $ \langle 0 \vert{\cal{O}}_{An} \vert J/\psi \rangle $ with $ n\geq 1 $ is unnecessary for exclusive processes. The dominant contributions of these quantities are purely kinematic and can be accounted for by the coefficients of $ \langle 0 \vert{\cal{O}}_{A0} \vert J/\psi \rangle $. Thus, the NRQCD factorization formula can be resummed to obtain

which is exactly the SGF formula: note Eq. (53), which relates $ \langle 0 \vert{\cal{O}}_{A0} \vert J/\psi \rangle $ to $ \overline{R}^{ {^{{3}}{S}_{{1}}^{[{1}]}} *}_{J/\psi} $.

We now show that, in general, SGF can be deduced from NRQCD at the operator level. The equations of motion of heavy quark fields in NRQCD are given by [2]

with a similar equation for $ \chi $ fields. Because we are not interested in gluon fields, we can replace $ D_0 $ by $ \nabla_0 $ and $ {{D}} $ by $ {{\nabla}} $. The equations of motion are usually used to replace operators involving $ \nabla_0 $ by operators involving $ {{\nabla}}^2 $ in NRQCD. Then, only spatial components $ {{\nabla}} $ appear in the NRQCD operators, which can be decomposed to a relative derivative and total derivative when acting on quark-antiquark bilinear operators. For example, beginning from the bilinear operator $ \chi^\dagger\psi $, in NRQCD, one can construct more operators by adding the relative derivative to obtain $ \chi^\dagger \overleftrightarrow{{\nabla}}^2 \psi $, the total derivative to obtain $ {{\nabla}}^2(\chi^\dagger \psi) $, or a combination of the derivatives.

However, the equations of motion can also be used to replace the relative derivatives $ \overleftrightarrow{{\nabla}}^2 $ and $ \overleftrightarrow{ \nabla}_0 $ by total derivatives^⑥, which results in the following matrix elements for exclusive processes:

where $ n_1 $ and $ n_2 $ are non-negative integers. As we are working in the rest frame of $ {\cal{Q}} $, by using integration by parts, we find that $ \nabla_0 $ gives rise to quarkonium mass $ M $, and $ {{\nabla}} $ vanishes. Therefore, in the factorization formula, the matrix element with $ n_1 = n_2 = 0 $ is enough to handle all the contributions in this series of matrix elements, although short-distance hard parts depend on the heavy quark mass$ m $as well as the quarkonium mass $ M $. This is simply the SGF formula. It is also clear that the SGF resums a series of power corrections originating from the kinematic effects in NRQCD.

For inclusive quarkonium processes, we can also use equations of motion to decompose the NRQCD matrix elements by

Using integration by parts, we can eliminate all the matrix elements except $ n_1 = n_2 = 0 $, with the short-distance hard parts depending on $ P^2 $, $ P\cdot P_X $, and $ P_X^2 $, where $ P_X $ is the momentum of unobserved particles $ X $. The final result is the SGF formula for inclusive quarkonium processes. It again resums a series of power corrections originating from the kinematic effects in NRQCD.

Because the short-distance hard parts depend on $ P_X $, nonperturbative matrix elements and soft gluon distributions must be also functions of $ P_X $ [7]. Note the difference between the soft gluon distributions and shape functions introduced in the endpoint region [20-23]. The purpose of the shape functions is to resum large logarithms in the endpoint region in the NRQCD factorization framework, which are defined at a fixed power in the relativistic expansion (usually leading power). The SGF with soft gluon distributions is designed to resum a power series of the relativistic expansion and can be applied at both the endpoint region and other regions. If SGF is applied at the endpoint region, the large logarithms can be naturally resummed by the renormalization group equations of the soft gluon distribution. A detailed discussion of this topic will be presented in a separate work [24].

In summary, taking $ \Gamma(J/\psi \to e^+e^-) $ as an example, we have demonstrated that the SGF is equivalent to the NRQCD factorization for heavy quarkonium production or decay. We have also shown that SGF can be deduced from NRQCD effective field theory at the operator level by using equations of motion. To arrive at this conclusion, we introduced a new regulator and defined the SGF matrix elements rigorously. Based on the equivalence between the two factorizations, we derived explicit relations between the SGF matrix elements and NRQCD matrix elements and proved the generalized Gremm-Kapustin relation.

The results obtained in this paper imply that, for any given process, SGF and NRQCD factorization are either both valid or both violated. Therefore, SGF is valid for all orders in perturbation theories for many processes in which NRQCD factorization has been proved. This provides a solid theoretical foundation for SGF. Compared with NRQCD factorization, SGF effectively resums a subset of relativistic correction terms originating from kinematic effects; this can reduce theoretical uncertainties and may thus provide a better description of experimental data.

We would like to thank K.T. Chao and C. Meng for useful discussions.

In this appendix, we use the method of regions to calculate the amplitude $ {\cal{A}}^{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}}) \to e^+e^-} $ at the one-loop order. According to Eqs. (24) and (25), for the original integral $ F $, we have

Defining

we have

As the first term on the r.h.s can be defined by integrals in the hard domain, it is infrared safe.

We now apply Eq. (A3) to calculate the QCD correction to $ {\cal{A}}^\mu_{c+\bar c} $. We first consider the vertex correction. As showed in Ref. [9], the vertex correction can be expressed in terms of elementary integrals $ I_{111}, I_{011}, I_{-111}, I_{010} $ , and $ I_{110} $, with $ I_{abc} $ defined as

According to Eq. (A3), $ I_{abc} $ can be expressed as

The overlap contributions $ I_{abc}^{(h,s)}, I_{abc}^{(h,p)}, I_{abc}^{(h,s,p)}, \text{ and } I_{abc}^{(s,p)} $ and the soft region contribution $ I_{abc}^{(s)} $ are scaleless integrals and can be set to zero if there is no infrared divergence. The calculation of $ I_{abc}^{(s)}-I_{abc}^{(s,p)} $ and $ I_{abc}^{(h,s)}-I_{abc}^{(h,s,p)} $ is similar to $ \overline{R}^{ {^{{3}}{S}_{{1}}^{[{1}]}}*,(s)}_{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})} -\overline{R}^{ {^{{3}}{S}_{{1}}^{[{1}]}}*,(s,p)}_{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})} $; we have

Then, we consider the contributions $ I_{abc}^{(h,p)} $. Consider $ I_{111}^{(h,p)} $as an example. The expanded integral reads

These integrals are scaleless and infrared safe and thus vanish. Similarly, we have

The remaining contributions $ I_{abc}^{(h)} $ and $ I_{abc}^{(p)} $ are given by

It should be noted that the ultraviolet and infrared divergences in $ I_{abc}^{(h)} $ are not well distinguished. However, when the overlap contributions $ I_{abc}^{(h,s)}, I_{abc}^{(h,p)}, \text{ and } I_{abc}^{(h,p,s)} $ are subtracted from $ I_{abc}^{(h)} $, the infrared divergent parts in the hard region can be removed. Then, the $ 1/\epsilon $ divergences in $ I_{abc}^{(h)}-I_{abc}^{(h,s)}- $$I_{abc}^{(h,p)} + I_{abc}^{(h,s,p)} $ are converted to $ 1/\epsilon_{UV} $ divergences. Inserting Eqs. (A6), (A8), and (A9) into (A5), we find that the results for the elementary integrals are consistent with those in Refs. [9,25].

We also need to calculate the heavy-quark wave-function renormalization $ Z_Q $, which is given by

where

Here, we have introduced integrals $ T_{ab}, T_{ab}^0, \text{ and } T_{ab}^1 $, which are defined by

Similarly, we can derive

For the remaining contributions, only $ T_{12}^{(s)}-T_{12}^{(s,p)} $ and $ T_{12}^{(h,s)}-T_{12}^{(h,s,p)} $ are non-vanishing; they read

Making use of Eqs. (A10), (A11), (A13), and (A14), we obtain

This agrees with the result in Ref. [26].

Based on the above results, one can reproduce the full QCD result in Eq. (15) by summing the contributions from the different domains:

with

Substituting Eq. (A16) into Eq. (18), $ {\cal{A}}^{c\bar{c}( {^{{3}}{S}_{{1}}^{[{1}]}})\to e^+e^-,(1)} $ can be expressed as

with $ {\cal{A}}^{(i,j,\cdots)} $ determined by $ {\cal{G}}^{(1),(i,j,\cdots)} $ and $ {\cal{H}}^{(1),(i,j,\cdots)} $.