Theory for quarkonium: from NRQCD factorization to soft gluon factorization

We demonstrate that the recently proposed soft gluon factorization (SGF) is equivalent to the nonrelativistic QCD (NRQCD) factorization for heavy quarkonium production or decay, which means that for any given process these two factorization theories are either both valid or both violated. We use two methods to achieve this conclusion. In the first method, we apply the two factorization theories to the physical process $J/\psi \to e^+e^-$. Our explicit calculation shows that both SGF and NRQCD can correctly reproduce low energy physics of full QCD, and thus the two factorizations are equivalent. In the second method, by using equations of motion we successfully deduce SGF from NRQCD effective field theory. By identifying SGF with NRQCD factorization, we establish relations between the two factorization theories and prove the generalized Gremm-Kapustin relations as a by product. Comparing with the NRQCD factorization, the advantage of SGF is that it resums the series of relativistic corrections originated from kinematic effects to all powers, which gives rise to a better convergence in relativistic expansion.

II. J/ψ → e + e − IN SOFT GLUON FACTORIZATION APPROACH A. Factorization formula According to Ref. [8], for exclusive decay process, one have the following SGF formula at the amplitude level, where n denote intermediate states. In general, n can contain dynamical soft partons (gluons or light quarks) in addition to a QQ pair. But for simplicity, we only discuss intermediate states without dynamical soft partons in this work, although dynamical soft partons can be discussed similarly. Then nonperturbative matrix elements R n * Q are defined as where Ψ stands for Dirac field of heavy quark, K n is projection operator defining the intermediate state n, and the subscript "S" means that, to evaluate the matrix elements, one only picks up integration regions where off-shellness of all particles is much smaller than heavy-quark mass. From the point view of method of regions [10,11], the effect of "S" keeps only small regions which are everything except the hard region. For process J/ψ → e + e − , symmetries of QCD tell us that only n = 3 S [1] 1 , 3 D [1] 1 are relevant, where we use the spectroscopic notation with superscript " [1]" denoting color singlet. We thus have with projection operators defined explicitly as where d = 4−2 is the space-time dimension, D µ is the gauge covariant derivative with Ψ ← → D µ Ψ = Ψ(D µ Ψ)−(D µ Ψ)Ψ, P is the momentum of J/ψ, µ Sz is a polarization vector with P · Sz = 0, m is the heavy-quark mass, M is the mass of J/ψ, the color projector is defined as C [1] = 1/ √ N c , and the spin projection operator P αβ is defined as The hard partsÂ n can be perturbatively calculated according to the matching procedure discussed in [7,8]. To this end, we first replace the J/ψ with a on-shell color-singlet state cc( 3 S which fix q 0 and |q| in the rest frame of P . The rest of degrees of freedom of q are removed by partial wave expansion, S-wave for this case. After the replacement, the l.h.s. of Eq. (3) becomes where Ω is the solid angle of relative momentum q in the cc rest frame, L µ is the leptonic current and A µ cc is the hadronic part of decay amplitude with spinors of cc removed. The cc pair is projected to state 3 S 1 by replacing spinors of cc pair by Similarly we have (11b) Based on these equations, one can calculate A cc Denoting perturbative expansion of any quantity W as W = W (0) + α s W (1) + α 2 s W (2) + · · · , we have the following orthogonal relations [7]: Based on this, the Eq. (3) results in the following matching relations: and up to order α s one has [9], and 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 the M in Ref. [9] is a free parameter. But to use these expressions for SGF, it needs to be the mass of quarkonium.

C. Perturbative calculation of matrix elements in SGF
Now we describe our method to calculate R . As was pointed out in Refs. [7,8], this quantity is defined to include only small loop momentum region. In the following, we will provide an explicit definition and choose a UV renormalization scheme.
Up to order α s , the corresponding Feynman diagrams are shown in Fig. 1, where the solid circle represents the operator ΨK 3 S [1] 1 Ψ. The calculation at tree level is straightforward, the result is Let's take the vertex correction Fig. 1(d) as an example to explain the calculation at one-loop level. The amplitude reads where T S is an operator that forces the loop momentum k to be in small region [7,8], which will be defined explicitly by using the method of regions [10,11]. Before continue, we note that although the full QCD integral in Eq. (20) is well regularized by dimensional regularization, the manipulation to define T S will generate unregularized integrals. Therefore, other regulator is needed to make our manipulation mathematically rigorous. We propose a new regularization at the full QCD level by multiplying the power of all propagator denominators by 1 + η, which can regularize all possible divergences encountered in the derivation using the method of regions, including both ultraviolet and non-ultraviolet divergences 2 . For the integral that we are interested in, we get where ν is introduced to compensate the mass dimension changed by the new regularization, and we assume η to make sure that the theory is eventually regularized by dimensional regularization.
With the new regularization in hands, we decompose the loop momentum k µ = (k 0 , k) into three domains: where relation " < ∼ " is understood as the negation of " ", D = R d is the complete integration domain, and an implicit cutoff scale exists to rigorously separate the three domains. This division satisfies Then for any original integral F , one can split it into the three domains 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 UV divergences removed by an MS renormalization scheme. T (s) and T (p) are conventionally called soft region and potential region, respectively, while the overlap region T (s,p) removes double counting between T (s) and T (p) .
For the soft region contribution R , we apply the operator T (s) to expand the integrand in Eq.(21) for small quantities, which results in the following integrals where the term k m1 0 (k · q) m2 comes from the numerator in Eq. (21) and with j i being non-negative integers. As scaleless and infrared-finite integrals can be set as zero (this is again a choice of scheme), we find only integrals with m 1 = m 2 = j 1 = j 2 = j 4 = j 5 = 0 are relevant, The above integrals have pinch poles around k 0 = 0, which can be regularized by the new regulator η. As integrals can be obtained by expanding the k 2 0 term in the denominator in Eq. (29), it is clear that they cancel exactly with the contribution from pinch poles around . This cancellation shows that the overlap region is important conceptually, although it contains only scaleless integrals that are usually set to zero in dimensional regularization. By combining the soft region and the overlap region, we only need to consider the contribution from the gluon pole and we can take η → 0 safely, which results in We emphasize that, if one wants to calculate R 1 ) separately, the correct order is to do the integration over k 0 first with fixed k. Or else, poles around k 0 = 0 will be regularized by different regulators between soft region and overlap region which necessarily breaks symmetries of the theory and makes the cancellation between the two regions impossible.
, terms in the expansion are proportional to where we can take η → 0 as all divergences are well regularized by dimensional regularization. Then because k 0 and k · q can be expressed as linear combinations of denominators, if any of m 1 , m 2 , j 1 , j 2 and j 3 is nonzero the integral can be decomposed to either scaleless and infrared-finite integrals or integrals with pure virtual value. By keeping only 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 MS renormalization procedure, which gives renormalized matrix element Similarly we can find the real part of R  (13), we obtain where We find the matrix element defined in SGF reproduces all infrared and Coulomb divergences in full QCD, and the obtained hard part is free of divergences. Therefore we conclude that the SGF factorization holds at least at one-loop level.

E. Validity of SGF at all orders in αs
The correctness of SGF at 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 A cc( 3 S [1] 1 )→e + e − is expressed as Where the first term on r.h.s has nonzero support only in the hard domain and thus is infrared safe. It is straight forward to check that the low energy part of A cc( 3 S [1] 1 )→e + e − has been correctly reproduced by matrix element in SGF, i.e., which leaves the corresponding short-distance hard part defined by high energy part of A cc( 3 S [1] 1 )→e + e − , and thus infrared safe. More precisely, we havê where MS means to remove UV divergences by MS subtraction scheme 3 . The above one-loop argument can be generalized to all orders. By definition, low energy part ("small" region of loop momenta) in full QCD can be reproduced by matrix element in SGF at any order in α s expansion, with a proper definition of T S at multi-loop level as discussed in Ref. [11]. Therefore, the short-distance hard part is perturbatively infrared-safe, which means the SGF formula for the decay width of J/ψ → e + e − holds to all orders in α s .

A. NRQCD result
Ignoring operators involving gauge fields, the NRQCD factorization for the decay amplitude J/ψ → e + e − is given by [9,13,14] A J/ψ→e + e − = n s n ( 3 S where ψ and χ † are the two-component heavy quark fields in NRQCD, 1 ) are short-distance hard parts which can be perturbatively calculated. The first two orders in α s expansion for s n ( 3 S [1] 1 ) are given by [9,13,14] with G , H 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,15,16] q 2n J/ψ = q 2 n J/ψ , with These relations are obtained by computing the matrix element 0|O An |J/ψ in potential-model [15].

B. Equivalence between SGF and NRQCD factorization
The basic reason for the validity of SGF in this process is that SGF matrix elements reproduce low energy part of full QCD. As NRQCD matrix elements also correctly reproduce low energy part of full QCD, the SGF is equivalent to NRQCD factorization. It means that, for any process, the two factorization formulas are either both valid or both broken.
Because the amplitude of J/ψ → e + e − can be factorized in both SGF and NRQCD factorization, we have (D-wave contributions are suppressed) which can generate relation between SGF matrix element and NRQCD matrix elements with O(v 2 ) denoting contributions from operators with explicit gauge fields. Actually, we can generate even more relations by applying the two factorization formulas to any well defined QCD quantity W , For example, if we choose W = 1 To determine corresponding w n , we replace J/ψ by a cc pair with invariant mass M 5 . Using the nonrelativistic expansion formulas given in [17], we have Here we used nonrelativistic normalization for the spinors u and v. Then we obtain Where the extra factor √ 2M in the first line appears because NRQCD matrix elements have nonrelativistic normalization, while SGF matrix elements have relativistic normalization. As both SGF matrix elments and NRQCD matrix elements keep the same low energy physics and renormalized in the same way, the coefficients w n ( 3 S [1] 1 ) should vanish at higher orders in α s , i.e.
Thus we get a relation which provides complete relations between the SGF matrix element and NRQCD matrix elements. Using these relations, we obtain which agrees with Eq. (45). We thus have proved the generalized Gremm-Kapustin relations by using the equivalence between SGF and NRQCD factorization. Based on our proof, it is clear that the O(v 2 ) terms in the relations are contributions from operators with explicit gauge fields.

A. Exclusive processes
Because the equivalence between SGF and NRQCD factorization gives rise to the generalized Gremm-Kapustin relations, it implies that Gremm-Kapustin-like relations are the key to relate NRQCD to SGF. Indeed, due to these relation, the introduction of nonperturbative quantities 0|O An |J/ψ with n ≥ 1 is unnecessary for exclusive processes. The dominant contributions of these quantities are purely kinematic that can be taken into account by coefficients of 0|O A0 |J/ψ . Thus the NRQCD factorization formula can be resummed to obtain which is exactly the SGF formula, noticing the Eq. (53) which relates 0|O A0 |J/ψ to R 1 * J/ψ . Now let us show generally that 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 similar equation for χ field. Because we are not interested in gluon fields, we can replace D 0 by ∇ 0 and D by ∇.
The equations of motion are usually used to replace operators involving ∇ 0 by operators involving ∇ 2 in NRQCD. Then only spacial components ∇ appear in NRQCD operators, which can be decomposed to relative derivative and total derivative when it acts on quark-antiquark bilinear operators. For example, beginning from the bilinear operator χ † ψ, in NRQCD one can construct more operators by adding relative derivative to obtain χ † ← → ∇ 2 ψ, or total derivative to obtain ∇ 2 (χ † ψ), or their combinations.
However, the equations of motion can also be used to replace relative derivatives ← → ∇ 2 and ← → ∇ 0 by total derivatives 6 , which results in the following matrix elements for exclusive processes with non-negative integers for n 1 and n 2 . As we are working in the rest frame of Q, by using integration by parts one finds that ∇ 0 gives rise to quarkonium mass M and ∇ vanishes. Therefore, in a factorization formula the matrix element with n 1 = n 2 = 0 is enough to take care of all contributions in this series of matrix elements, although short-distance hard parts will depend on heavy quark mass m as while as quarkonium mass M . This is nothing but the SGF formula. It is also clear that the SGF resums a series of power corrections originated from kinematic effects in NRQCD.

B. Inclusive processes
For inclusive quarkonium processes, we can also use equations of motion to decompose NRQCD matrix elements by Using integration by parts we can eliminate all matrix elements except n 1 = n 2 = 0, but with short-distance hard parts depending on P 2 , P · P X and P 2 X , 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 originated from kinematic effects in NRQCD.
Because short-distance hard parts depend on P X , nonperturbative matrix elements, soft gluon distributions, must be also functions of P X [7]. Note the difference between soft gluon distributions and shape functions introduced at endpoint region [18][19][20][21]. The purpose of shape functions is to resum large logarithms at endpoint region in NRQCD factorization framework, which are defined at fixed power in relativistic expansion (usually leading power). The SGF with soft gluon distributions are aimed at resumming a power series of relativistic expansion, which can both be applied at endpoint region and regions away from that. If SGF is applied at endpoint region, large logarithms can be naturally resummed by renormalization group equations of soft gluon distribution. Detailed discussion of this topic will be presented in a separate work [22].

V. SUMMARY
In summary, taking Γ(J/ψ → e + e − ) as an example, we demonstrated that the SGF is equivalent to the NRQCD factorization for heavy quarkonium production or decay. We also shown that the SGF can be deduced from NRQCD effective field theory at the operator level by using equations of motion. To achieve the above conclusion, we introduced a new regulator and defined SGF matrix elements rigorously. Based on the equivalence between the two factorizations, we derived explicit relations between SGF matrix elements and NRQCD matrix elements and proved the generalized Gremm-Kapustin relations.
The results obtained in this paper means that, for any given process, the SGF and the NRQCD factorization are either both valid or both violated. Therefore, the SGF is valid to all orders in perturbation theories for many processes where NRQCD factorization have been proved . This provides a solid theoretical foundation for the SGF. Comparing with the NRQCD factorization, the SGF effectively resums a subset of relativistic correction terms originated from kinematic effects, which can reduce theoretical uncertainties and thus may provide a better description of experimental data.
As the first term on the r.h.s can be defined by integrals in hard domain, it is infrared safe. Now we apply the formula Eq. (A3) to calculate the QCD correction to A µ c+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