New factorization theory for heavy quarkonium production and decay

The widely used nonrelativistic QCD (NRQCD) factorization theory now encounters some notable difficulties in describing quarkonium production. This may be due to the inadequate treatment of soft hadrons emitted in the hadronization process, which causes bad convergence of velocity expansion in NRQCD. In this paper, starting from QCD we propose a rigorously defined factorization approach, soft gluon factorization (SGF), to better deal with the effects of soft hadrons. After a careful velocity expansion, the SGF can be as simple as the NRQCD factorization in phenomenological studies, but has a much better convergence. The SGF may provide a new insight to understand the mechanisms of quarkonium production and decay.

The widely used nonrelativistic QCD (NRQCD) factorization theory now encounters some notable difficulties in describing quarkonium production. This may be due to the inadequate treatment of soft hadrons emitted in the hadronization process, which causes bad convergence of velocity expansion in NRQCD. In this paper, starting from QCD we propose a rigorously defined factorization approach, soft gluon factorization (SGF), to better deal with the effects of soft hadrons. After a careful velocity expansion, the SGF can be as simple as the NRQCD factorization in phenomenological studies, but has a much better convergence. The SGF may provide a new insight to understand the mechanisms of quarkonium production and decay.

I. INTRODUCTION
As the simplest bound state of strong interactions, heavy quarkonium is a perfect system to study both perturbative and nonperturbative physics of QCD. Ever since the discovery of the first heavy quarkonium J/ψ in 1974, a lot of efforts have been devoted to interpret the production and decay mechanisms of heavy quarkonium. Among them, the most notable theories include the color-singlet model [1][2][3] and the nonrelativistic QCD (NRQCD) factorization theory [4] 1 .
best fit of J/ψ yield data at high transverse momentum in hadronic collisions determines two linearly combined LDMEs, M 0 = 0.074GeV 3 and M 1 = 0.0005GeV 3 [22] (the J/ψ polarization data requires almost the same two combined LDMEs [38]). There is a two-orders difference between the two combined LDMEs. However, velocity scaling rules in NRQCD [4] expect these LDMEs to be at the same order of magnitude. (3) The universality problem: A necessary condition for NRQCD factorization is that LDMEs are universal, i.e. process independent. Yet phenomenology studies show that M 0 extracted from hadron colliders [22,23,39] is much larger than the upper bound set by e + e − collisions M 0 < 0.02GeV 3 [45].
In this paper, we show that the velocity expansion in the present NRQCD framework suffers from large high order relativistic corrections due to ignoring the momentum of soft hadrons, which are mainly soft gluons perturbatively, emitted in the hadronization process. Thus, by including only a few low order contributions in relativistic expansion, NRQCD is hard to provide good descriptions for quarkonium production, which may be the reason for the above mentioned difficulties. Starting from QCD, we propose a new factorization approach, called soft gluon factorization (SGF), to better deal with the effects of soft hadrons. The SGF has a much better convergence in the velocity expansion, and it may provide a new insight to fully understand the mechanisms of quarkonium production and decay.
The rest of the paper is organized as follows. In Sec. II, we explain why NRQCD factorization may have bad convergence in velocity expansion, and then define the SGF formula with rigorous operator definition. Relations to other approaches are also discussed. In Sec. III, we do a careful simplification to the SGF so that it will be easy to use in practice, but at the same time it can capture main physics. Following this, we compare the SGF and NRQCD by studying the gluon fragmentation to yield of J/ψ in Sec. IV. As expected, we find that the lowest order NRQCD approximation is not good, which may be the reason why NRQCD encounters many difficulties. The summary and the possibility of solving these difficulties in the SGF framework are given in Sec. V.

II. SOFT GLUON FACTORIZATION
A. Bad convergence of v 2 expansion in NRQCD In NRQCD factorization [4], the differential cross section of a heavy quarkonium H production can be factorized as where P H is the momentum of H, n denotes intermediate QQ states, whose quantum numbers are usually expressed in terms of spectroscopic notation 2S+1 L [c] J , with c = 1, 8 denotes color singlet or color octet of the pair, H n (P H ) are perturbative calculable short-distance coefficients which can be expanded order by order in α s , O H n are gauge-completed [25][26][27] nonperturbative longdistance matrix elements (LDMEs) which can be classified according to the power counting rules of v, and the ellipses denote other relativistic correction terms. It is needed to emphasize that, in NRQCD factorization, one expands the mass of quarkonium M H around twice of heavy quark mass 2m, which results in P 2 H ≈ 4m 2 in Eq. (1).
If convergence of velocity expansion in NRQCD factorization is very good, one only needs to consider a very limited number of LDMEs to describe experimental data. Unfortunately, it may not be the case. To see this, for an example let us consider the differential cross section of quarkonium production in hadron colliders with high transverse momentum P HT . In NRQCD factorization, one calculates the partonic QQ production differential cross sections with transverse momentum P T , and then expands P T around P HT . On average, let us set P HT ∼ (1 − λ)P T , with λ being the average transverse momentum fraction carried by soft gluons which can be at the order of v or v 2 . At a specific range of large P HT , differential cross section behaves as with c being usually larger than 4. Expansion of this result with respect to λ mimics the relativistic expansion in NRQCD. For a typical choice λ = 0.2 and c = 4 [46], we have (1 − λ) c = 0.4096, but the expansion of (1 − λ) c to O(λ 0 ) gives 1, and expansion to O(λ 1 ) gives 0.2. The above argument shows a bad convergence of relativistic expansion. We note that the convergence still cannot be improved if one instead is to calculate, e.g., P 4 HT dσH dP 2

HT
. Furthermore, as we will see in Sec. IV, there are other relativistic correction terms that make the convergence of relativistic expansion even worse.
Because of the bad convergence of relativistic corrections, theoretic calculations based on only a limited number of LDMEs are sometimes hard to describe experimental data. In literature, e.g. Refs. [47][48][49][50], for each problem one can resum a specific subset of LDMEs to improve the theoretic results. The resummation results in many so-called "shape functions".
A factorization method aiming to improve convergence for an arbitrary problem was proposed in Ref. [51]. However, this method can only be thought of as a model because there is no operator definition for nonperturbative functions. Without an operator definition, the method is not well defined and it is not possible to do rigorous calculation beyond tree level.

B. Soft gluon factorization formula
The aim of soft gluon factorization is to resum a subset of relativistic correction terms in NRQCD factorization that are important for a phenomenological purpose. It is convenient to demand the subset to be Lorentz invariant. The formula of SGF for a quarkonium H production is must be produced in the hadronization process. In this case, the corresponding F n→H (P, P H ) has short-distance effects and can be refactorized, which enables us to eliminate this contribution from the factorization formula. Furthermore, if the intermediate state contains other energetic partons in addition to a QQ pair, the corresponding F n→H (P, P H ) can also be refactorized because there is no energetic gluon or light quark in the dominant Fock state of conventional quarkonium. In this consideration, we only need to sum over intermediate states which contain a QQ pair and some soft partons. According to physical scales existing in the quarkonium system, energy of these soft partons can be at the order of mv or mv 2 , which are assumed to be much smaller than m.
As each conventional quarkonium has a QQ state with specific quantum numbers as its leading Fock state, the QQ intermediate state with the same quantum numbers has leading contribution. If the intermediate state contains soft partons in addition to a QQ pair, which means that soft partons are produced in the hard parts, the contribution is at least suppressed by v 2 comparing with the contribution from the intermediate state with the same QQ but without soft partons. In other words, power counting of the corresponding F n→H (P, P H ) in the factorization formula is at high order in v 2 . Unlike v 2 corrections due to kinematic effects of soft gluons emission, contribution of intermediate states with soft partons should have a mild v 2 corrections. Thus, as the first approximation, we will ignore these contributions and consider only intermediate states constituted by a QQ pair in the rest of the paper.
Then n in Eq. (3) denotes intermediate QQ states, whose quantum numbers are usually expressed in terms of spectroscopic notation n = QQ( 2S+1 L

[c]
J,Jz ), with c = 1, 8 denoting color singlet or color octet of the pair. We note that the state in the amplitude n can be in general different from that in the complex-conjugate amplitudẽ J,Jz ). Color charge, C-parity, P -parity and angular momentum conservations provide general selection rules [52], which gives that the following relations always hold: c =c, S =S, and |L −L| = 0, 2, 4, · · · . For the production of a polarization-summed quarkonium, we have further constraints J =J, J z =J z . Thus, even in the polarization-summed case, there is interference between the 3 S  1,Jz state. Selection rules for polarization distributions can be found in Ref. [52]. For simplicity, in the rest of the paper we will only discuss polarization-summed quarkonium production and choose the quantum numbers of the intermediate QQ pair in the amplitude to be the same as that in the complex-conjugate amplitude, while polarization and interference contributions can be studied similarly.

D. Nonperturbative distributions
To make SGDs Lorentz invariant, it is convenient to define them using QCD fields instead of NRQCD fields. We define SGDs for polarization-summed quarkonium H production as vacuum expectation values of bilocal operators constructed from QCD fields, where projects final state to include a polarization-summed H with relativistic normalization H( The subscript "S" means that, to evaluate the matrix element, one only picks up integration regions where off-shellness of all particles is much smaller than heavy quark mass 2 . K n (rb) are projection operators defining the intermediate state n, where angular momentum operators Γ n are defined explicitly in Appendix A, color operators C [c] will be defined later, MH ± / P H 2MH project out large components in velocity expansion, and, as we will see later, the introduction of prefactor √ MH MH +2m makes our SGDs more closely related to NRQCD LDMEs.
If the intermediate QQ pair is in color singlet, we define color operator as which is nothing but SU (3) color Clebsch-Gordan coefficient 3i;3ī|00 and 1 c is the identity matrix in color space. If the intermediate QQ pair is in color octet, we define color operator by the multiplication of a color matrix in fundamental representation with a gauge link in adjoint representation, where √ 2tā is SU(3) color Clebsch-Gordan coefficient 3i;3ī|8ā and the introduction of gauge link Φ The next question is what direction of the gauge link should we choose. A natural choice is to define the direction along a light cone, because a gauge link can be obtained by eikonal approximation for soft gluon interaction between the QQ pair and jets. Indeed, this is the choice to define gauge-completed NRQCD LDMEs in Refs. [25][26][27]. In the SGF framework, however, this choice will result in uncanceled gauge-link-collinear rapidity divergence because the position b in Eq. (4) is in general not along a light cone. This is the same problem encountered in transverse-momentum-dependent (TMD) factorization where there is a mismatch of rapidity divergence between virtual correction and real correction [54]. To cancel the rapidity divergence, soft factor as while as jet functions may need to be introduced.
If we do not want to introduce additional soft factors 3 , the gauge links can be defined along the b direction, where P denotes path ordering, A (A) is gluon field with color matrix in adjoint representation, and b µ ℓ = b µ + εℓ µ . We choose 0 < ε ≪ 1 and a fixed light-like vector ℓ µ so that: when b µ is finite, b µ ℓ is the same as b µ ; while as b µ → 0, the gauge link is well defined (along the ℓ µ direction).
In the above, we have defined SGDs for n with fixed J z . SGDs for polarization-summed intermediate states can be defined accordingly, Similar to LDMEs in NRQCD, each SGD has a definite power counting in v 2 , which will be explained later. Thus, for a certain accuracy, we can truncate to use only a few SGDs for each H production.

E. Perturbative matching
With the definition of SGDs, short-distance hard parts can be obtained order by order in perturbation theory by matching both sides of Eq. (3). As H n (P ) are independent of nonperturbative physics, they are unchanged if we replace the quarkonium H in Eq. (3) by an on-shell QQ pair with specific quantum numbers. To make the following equations more convenient to use, we first relabel n and P in Eq. (3) by n ′ and P ′ , respectively, and we then replace H in both sides of the equation by a state J,Jz ) with momentum P , which results in where Note that, although dσ H and F n ′ →H are changed to dσ n and F n ′ →n , respectively, H n ′ (P ′ ) stays unchanged, which is guaranteed by the fact that H n ′ (P ′ ) is perturbatively calculable and thus it is independent of nonperturbative external states. The reason for choosing the on-shell QQ pair is that dσ n and F n ′ →n are then gauge invariant and can be expanded order by order in α s in any gauge.
The key to determine all H n ′ is to project H to a complete set of QQ states, and then to match coefficients of α s in both sides of Eq. (3). Although any complete set is doable for this purpose, a good choice should satisfy as much as possible that leading order in α s expansions of F n ′ →n are delta functions.
To define a good complete set, we begin with a QQ with momenta We project it to color state c by color projectors: where N 2 c − 1 in the denominator is to average over color-octet states. We then project it to state with total spin S and S z by replacing spinors of QQ by 4 where M 2 = P 2 and On-shell conditions p 2 There are still 2 degrees of freedom of the relative momentum q, which can be chosen as spatial angles in the rest frame of the pair. We further do a partial wave expansion in this frame, which removes the q dependence although introducing orbital angular momentum L and L z . Finally, we project spin and orbital angular momentum to total momentum J and J z . More precisely, in the rest frame of P , dσ n and F n ′ →n with quantum J,Jz ) can be obtained from the corresponding production of QQ with fixed c, S and S z by the following operation: 4 This expression can be interpreted as covariant spin projectors [55,56] multiplied by a normalization factor 2/M [57]. But for our purpose, this interpretation is unnecessary. The key in our definition is to guarantee the orthogonal relations Eq. (19).
In Appendix B, we show that the above definition can indeed result in orthogonal relations at lowest order where the superscript "(0)" denotes leading order in α s expansion. By inserting perturbative expansions n ′ →n + · · · , (20b) into Eq. (11) and using the orthogonal relations Eq. (19), coefficients of different powers of α s give the following relations: and so on. Based on these relations, H n can be obtained by perturbative calculation of dσ n and F n ′ →n for on-shell QQ pair production. This perturbative calculation relies on the cancellation of IR divergences between dσ n and F n ′ →n [like in Eq. (21)] to all orders in perturbation theory, which is a very difficult problem and will be discussed in the next subsection. It is worth emphasizing that the above perturbative calculation of on-shell QQ pair production can only directly determine H n (P ) with P 2 > 4m 2 . Although on average we may have P 2 > 4m 2 [7], there is a nonvanishing contribution from the P 2 < 4m 2 region in Eq. (3). The value of H n (P ) in the later region can be obtained by analytical continuation of its value in the former region. Analytical continuation here is almost trivial, which means that we use the same functional form of H n (P ) in all regions of P 2 .

F. Justification
The above proposed SGF is only a conjecture, and we cannot provide an all order proof at present. The most dangerous interaction that may ruin the factorization is the elastic scattering between the QQ pair and hard jets by exchanging gluons, the kinematic region of these exchanged gluons are called "Glauber" region in literature. We find that IR divergence from the "Glauber" region at one-loop order is purely imaginary, which cancels with the contribution from its complex conjugate diagram. Therefore, SGF holds at one-loop order. Thanks to the two-loop study of infrared divergences of quarkonium production [25][26][27], it is hopeful that SGF may also hold to two-loop order, which however still needs further examination.

G. Relation to NRQCD factorization
The difference between NRQCD and SGF can be partially understood from the treatment of intermediate momenta. Based on relations in Eq. (17), we can express q ′ 0 and q ′2 by P ′ 0 and P ′2 , which corresponds to the SGF strategy. Alternatively, we can also use Eq. (17) to express q ′ 0 and P ′ 0 by q ′2 and P ′2 , which corresponds to the strategy of NRQCD factorization.
To further discuss the relation between SGF and NRQCD, by ignoring velocity corrections we can approximate H n (P ) by H n (P H ) in Eq. (3). Then the integral over P is applied on SGDs defined in Eq. (4), and thus Eq.
which has the same form as the NRQCD factorization Eq.(1). In the above formula, are closely related to gauge-completed LDMEs O H n although they are defined by different fields. Considering the differences between normalization in color space, up to lowest order in v 2 approximation, we have 5 Based on NRQCD velocity scaling rules, these relations also tell us v 2 power counting rules for SGDs. Beyond the above approximation, the complete SGF formula resums a series of velocity corrections in 5 Let us first check color-octet states. In NRQCD one usually chooses nonrelativistic normalization for the H state, which means that our definition has an extra factor 2E H due to normalization of state. For a color operator, in NRQCD one uses tā instead of √ 2tā, due to which we have an extra factor 2. Finally, the prefactor in Eq. (6) introduces an extra factor M H /(M H + 2m) 2 . The multiplication of these three factors equals 1 if we choose the H rest frame and take v 2 → 0. Except these differences, all other parts of our definition are identical to the NRQCD definition at the lowest order in v 2 approximation. For color-singlet states, in NRQCD one uses 1c as color operator, which results in a difference of factor 2Nc comparing with the color-octet case. NRQCD, which are important for phenomenological study. As we know, the TMD factorization is a generalization of collinear factorization with a series of higher twist corrections resummed [54]. In this sense, we can say that SGF is a "TMD" version of NRQCD. We expect that there are a lot of similarities between SGF and TMD factorization.

H. Relation to QCD collinear factorization
It is also needed to explain the relation between the SGF method and the QCD collinear factorization method [8][9][10]. In QCD collinear factorization, one assumes the transverse momentum of the produced quarkonium to be much larger than the heavy quark mass, P HT ≫ m. Therefore, one can factorize out hard physics at the scale of P HT , and leave the physics at the scale of m to be described by input functions, which are called fragmentation functions (FFs). Unlike NRQCD factorization and SGF, QCD collinear factorization has been proved to all orders in perturbation theory at both leading power and next-to-leading power levels. Phenomenologically, the advantage of this method is that large logarithms log(P HT /m) can be resummed to all orders in perturbation theory using renormalization group equations for single-parton and double-parton FFs [10], and therefore there can be smaller theoretical uncertainties.
Due to the above advantages, when studying high transverse momentum quarkonium production, one should better first use QCD collinear factorization to express cross sections in terms of the convolution of perturbative calculable hard parts [58] with FFs. Based on renormalization group equations, these FFs are fully determined by their values at an initial scale. Because of the heavy quark mass, these FFs at the initial scale still have perturbative physics that can be further separated from nonperturbative physics. This can be achieved by using NRQCD or SGF. Using NRQCD factorization, these FFs at an initial scale have been studied extensively [59][60][61][62][63][64][65][66][67][68][69][70]. In Sec. IV, we will see that FFs can also be calculated by using the SGF method, which can improve the convergence of v 2 expansion comparing with the using of NRQCD factorization.
Note, however, that both NRQCD factorization and SGF can also be used to directly describe low transverse momentum quarkonium production as well as quarkonium decay, where the QCD collinear factorization does not apply. where M 2 = P 2 and λ ∼ a few hundreds MeV << M is the typical energy of emitted gluons. Thus we can expand O(λ) and O(λ 2 /M ) terms in hard parts. In this expansion, although space components of P µ are O(λ), the convergence should be much better than usual velocity expansion in NRQCD because P is integrated out symmetrically around the origin. The leading term in this expansion gives where In this way, we significantly reduce the complication of input functions from four dimensional to one dimensional. We refer to Eq. (25) as SGF-1d and Eq. (3) as SGF-4d. Note that a similar formula that relates momentum of χ cJ to its decaying particle J/ψ, p J/ψ ≈ m J/ψ mχ cJ p χcJ , was first introduced by us in Ref. [21], and later also used in other publications [7,39,71]. The approximation was found to be very good, with only 8% deviation from full results for both yield and polarization measurements [71]. We thus expect the approximation Eq.(25) to be also very good.
In the rest of the paper, we will only take Eq. (25) as an example to show the effect of collinear approximation.

B. Statical approximation
As hard gluon emissions are excluded from F n→H (z) and very soft gluon emissions are suppressed by phase space, F n→H (z) should peak around z = z n ∼ 1 − O(λ/m H ). If the distribution of F n→H (z) is narrow enough, we can approximate it by δ(z − z n ) O H n . Then the SGF becomes which we call SGF-0d. While this formula is as simple as NRQCD factorization, we will see later that the lowest order SGF-0d is a much better approximation of real physics than that of the lowest order NRQCD.

C. Expansion of m
In hard parts of all the above SGFs, there are at least two independent hard scales, 2m and M . Their difference is, however, much smaller than 2m. We thus can further simplify hard parts by expanding m around M/2, which defines a relativistic correction series in SGFs. We will see that there are no large corrections in this expansion. To compare SGFs with NRQCD factorization, we apply them to the J/ψ hadron production via gluon fragmentation. A charmonium H production cross section via gluon fragmentation is given by where dσ g (p T /x) is the well-known hard part that produces a gluon with transverse momentum p T /x, and D g→H (x) is the FF of a gluon into an H that will be determined in the following.

C. Model assumptions
For the nonperturbative function F n→H (P, P H ), we simply employ an existing model in Ref. [51], 1 →H (P, P H ) = a k 2 exp(− where k = P − P H is a timelike momentum with positive energy k 0 . Λ is an energy cutoff for emitted soft gluons, for which we will choose 500MeV 6 . As we are only interested in the cross section ratio in the following, the value of a is irrelevant. We further set m = 1.55GeV in studying J/ψ production. We note that, our conclusions are in fact independent of these choices.

D. Numerical results
With these inputs in hand, FFs calculated by using different methods are shown in Fig. 1. In this figure, D 4d g→H (x) and D 1d g→H (x) are calculated by using the exact value of m, while D 1d (0) g→H (x) is similar to D 1d g→H (x) but calculated by setting m = M/2 in Eq. (34). We find that FFs calculated by using SGF-4d and SGF-1d have different shapes, but they have the same accumulated values, This is because SGF-4d and SGF-1d are equivalent for the integrated FF. The FF obtained by expanding m to the lowest order, D 1d (0) g→H (x), has only a small difference from the complete result D 1d g→H (x), which implies that relativistic correction due to the expansion of m is very small. We also find that F 3 S [8] 1 →H (z) has almost the same shape as that of D 1d g→H (x), which is δ-function-like. For cross section, we will find that SGF-0d with z 0 = 0.86 can well reproduce SGF-4d. Using the above FFs, we can calculate the J/ψ cross section based on Eq. (30). Let us assume that the cross section calculated by SGF-4d is "exact", and examine how good are the SGF-1d, SGF-0d, and NRQCD expansions. We show ratios of J/ψ differential cross section calculated by different methods over that calculated by SGF-4d in Fig. 2. We find that the SGF-1d result is close to the exact value, with largest error less than 6%. This implies that, as expected, the SGF-1d expansion should be very good. By expanding m to the lowest order, SGF-1d (0) also provides a good approximation, with deviation smaller than 10%. The SGF-0d can well reproduce SGF-4d if we choose z 0 = 0.86, as shown in Fig. 2. So, we expect that the convergence of velocity expansion in SGF is good in general.
On the contrary, the lowest order NRQCD result is larger than the exact value by more than a factor of 4, where NRQCD LDMEs are determined by the approximation Eq. (24). There are two main sources for the large deviation. One comes from the hard part of Eq. (30), which has an approximate scaling behavior dσ g (p T /x) ∼ (p T /x) −4 with average value of x being around 0.86. Yet, by ignoring soft gluon emission, NRQCD approximates x ≈ 1, which enhances the total result by almost a factor of 2. This effect has also been pointed out in Sec.II A. The other one comes from Eq. (34). By ignoring soft gluon emission and then expanding M around 2m, NRQCD approximates ∆ ≈ m J/ψ /0.86 2m ≈ 1.16 by 1, which enhances the total result by another factor of 2. Roughly speaking, the lowest order NRQCD approximates 0.86 9 ≈ 1/4 ∼ (1 − v 2 /2) 9 by 1 in this problem, which is hard to be recovered by traditional relativistic correction in NRQCD [73].
For bottomonia production, because v 2 is smaller, NRQCD approximation can be a little bit better than that for charmonia. According to the above discussion, we expect the lowest order NRQCD approximation overestimates bottomonia production rate by a factor of (1 − v 2 /2) −9 ≈ 1.6.

V. SUMMARY AND OUTLOOK
In this paper, we propose a soft gluon factorization (SGF) approach to describe quarkonium production and decay, which keeps the momentum difference between intermediate QQ pair and physical quarkonium. From the point view of NRQCD, the SGF effectively resums a subset of Lorentz-invariant relativistic correction terms in NRQCD factorization. In this sense, the SGF is a generalization of the NRQCD factorization.
By construction, the SGF has a good convergence in velocity expansion, so the lowest order approximation may already capture most physics; whereas for NRQCD factorization, we find that the lowest order result in velocity expansion can deviate from the full result by more than a factor of 4, mainly due to ignoring the momentum taken away by soft gluon emission. With so large deviation, it is not surprising that the NRQCD calculation faces many difficulties.
Hopefully, these difficulties may be resolved or relieved in the SGF framework with well controlled relativistic corrections. Specifically, the universality problem may be due to the fact that, for quarkonium production in different processes, e.g., e + e − collision or pp collision, lowest order NRQCD calculations suffer from large but different relativistic corrections. While in SGF, we do not expect large relativistic corrections. Moreover, considering the large relativistic effect on yield, we may also expect a significant effect on polarization, like color-magnetic dipole transition effect, which may alter the transverse polarization of the 3 S As total spin has only two choices, S = 0 or 1, the corresponding spin operators are where ǫ µ Sz are polarization vectors. Orbital operators are defined as where D µ is the gauge covariant derivative withΨ ← → D µ Ψ = Ψ (D µ Ψ) − D µΨ Ψ and ǫ µ1···µL Lz are L-rank polarization tensors. The above polarization tensors ǫ µ Sz and ǫ µ1···µL Lz have only spatial components in the rest frame of P Hµ , which is equivalent to the following relations: Furthermore, polarization tensors are chosen to be orthonormal, e.g., ǫ Sz · ǫ * S ′ z = −δ SzS ′ z . The combination of the above two properties results in Although orbital operators in Eq. (A3) are similar to that in NRQCD, their meanings are actually very different. In fact, in addition to these operators, there are an infinite number of relativistic-correction operators in NRQCD, like (− i 2 ) 2 ← → D 2 , which do not show up in SGF. The reason is that, instead of relativistic expansion, the SGF is performing a partial-wave expansion. This can be seen more clearly in momentum space. Then in the rest frame of P H , Eq. (A3) becomes where q µ is half of the relative momentum between intermediate QQ, and θ and φ are polar angle and azimuthal angle of q, respectively. Equation (A6) can be thought of as the definition of Γ o LLz in momentum space. As the partial-wave expansion already forms a complete set of operators, we do not need further relativistic-correction operators. Relativistic corrections in SGF are encoded in perturbative calculable hard parts. Our definition of SGDs is to pick up the minimal number of |q| in each partial wave.
where (2π) 4 δ 4 (p Q − p ′ Q ) presents because the diagram is disconnected and the terms inside of "Tr" project both initial QQ and final QQ to specific spin states. By setting P = P ′ and using the delta function to integrate out q, we get (B2) where we have used the last line of Eq. (B2) gives where we have used the fact that P H = P at this order and thus P · ǫ Sz = 0. Therefore, Eventually, we have