Fragmentation functions of polarized heavy quarkonium

Study of the polarized heavy quarkonium production in recently proposed QCD factorization formalism requires knowledge of a large number of input fragmentation functions (FFs) from a single parton or a heavy quark-antiquark pair to a polarized heavy quarkonium. In this paper, we calculate these FFs at the input scale in terms of nonrelativistic QCD (NRQCD) factorization. We derive all relevant polarized NRQCD long-distance matrix elements based symmetries and propose a self-consistent scheme to define them in arbitrary $d$ dimensions. We then calculate polarized input FFs contributed from all $S$-wave and $P$-wave NRQCD intermediate states. With our calculation of the polarized input FFs, and the partonic hard part available in literatures, the QCD factorization formalism is ready to be applied to polarized heavy quarkonium production.


Introduction
Since the discovery of the J/ψ, heavy quarkonia have been serving as ideal systems to test our understanding of QCD bound states and the hadronization processes. Unfortunately, a theoretically and phenomenologically satisfying framework is still lacking for heavy quarkonium production. The problem is more acute with the recently discovered XY Z-mesons, since the production of these exotic mesons requires good understanding of the production of conventional heavy quarkonium [1,2]. A phenomenologically successful model for heavy quarkonium production is the non relativistic QCD (NRQCD) factorization [3,4], which expands the cross section by powers of α s and v, the relative velocity of the heavy quark pair inside of a heavy quarkonium. All nonperturbative contributions are included in a set of NRQCD long-distance matrix elements (LDMEs). If the factorization is correct to all order in α s , these LDMEs will be universal, i.e. process independent. For heavy quarkonium production, although a formal proof of NRQCD factorization to all orders in α s is still lacking, it has been shown to be valid at the next-to-leading order (NLO) in many processes. With certain modification, it also works at the next-to-next-to-leading order (NNLO) in some specific cases [5][6][7]. Phenomenologically, NLO NRQCD factorization calculations successfully explain the p T distribution of many heavy quarkonium states with a few NRQCD LDMEs fitted from the data [1,2].
Nevertheless, the current NRQCD factorization formalism is far from perfect in the high transverse momentum, p T , heavy quarkonium production. It is understood now that fixedorder NRQCD calculation suffers from large high-order corrections, due to the large power enhancement in forms of p 2 T /m 2 Q and large logarithm log(p 2 T /m 2 Q ) at high orders [8], where m Q is the mass of heavy quark. For example, the NLO NRQCD calculation for the yield of J/ψ is orders of magnitude larger than the leading order (LO) calculation for 3 S [1] 1 and 3 P [8] J channels at p T 15 GeV [9][10][11][12]. Although the most phenomenologically important leading power contribution is claimed to be included in the NLO NRQCD calculation [11,13], the existence of large logarithm may potentially undermine the convergence of α s expansion at large p T m Q . Recently, a new QCD factorization formalism has been proposed to study heavy quarkonium production at large p T [8,[14][15][16][17][18]. In this formalism, the cross section is expanded by powers of m 2 Q /p 2 T , with both leading-power (LP) and next-to-leading power (NLP) terms considered. It is proved to all orders in α s that the dominant LP (NLP) term can be factorized into the production of a single parton (a heavy quark pair) in the hard part, convoluted with the single parton (the heavy quark pair) fragmentation functions (FFs) to heavy quarkonium. All nonperturbative contributions are included in these FFs, whose scale dependence is determined by a closed set of evolution equations. By solving the evolution equations, large logarithms log(p 2 T /m 2 Q ) are resummed to all orders in α s . Because of the systematic treatment of the powers of p 2 T /m 2 Q and log(p 2 T /m 2 Q ), the QCD factorization is expected to converge faster in α s expansion than the NRQCD.
With partonic hard part calculated in Ref. [18], and evolution kernels calculated in Ref. [8,16,17], a set of FFs at some input scale µ 0 is still required for phenomenological application of the QCD factorization formalism. These input FFs are the only components in the QCD factorization framework which are sensitive to the observed heavy quarkonium states and their polarizations. In principle, these input FFs should be extracted from data, like the case for pion and kaon. Nonetheless, this task is extremely hard in practice, especially with the inclusion of NLP contributions. For example, polarization-summed J/ψ production needs at least 4 one-variable single parton FFs and 6 three-variable heavy quark pair FFs [19,20]. The number of FFs doubles for polarized J/ψ production.
However, different from pion and kaon FFs, the heavy quarkonium FFs have an intrinsic large scale m Q Λ QCD , thus they are partially perturbative. Since 1993, Braaten, Cheung and Yuan proposed to apply the Color-singlet model and the NRQCD factorization to further separate the perturbative and nonperturbative contributions in the input FFs [21,22]. By choosing the input QCD factorization scale µ 0 2m Q and NRQCD factorization scale µ Λ ∼ m Q , neither µ 2 0 /m 2 Q nor m 2 Q /µ 2 Λ is large, consequently the NRQCD expansion of the input FFs is expected to have fast convergence. In Refs. [19,20], we calculated polarization-summed LP and NLP input FFs for all S-wave and P -wave channels at the NLO in α s with NRQCD factorization. 1 With our calculated input FFs, a first test of QCD factorization for heavy quarkonium production is shown in Ref. [23]. It was found that without the resummation of large logarithms, our simple and analytical LO QCD factorization calculation successfully reproduces the complicated and numerical NLO NRQCD calculation for all phenomenologically important channels for the yield of J/ψ at p T 15 GeV. Surprisingly, the NLP contribution, although suppressed by m 2 Q /p 2 T compared to LP contribution, are crucial for 1 S [8] 0 and 3 S [1] 1 channels even if p T approaching 100 GeV. A similar calculation with only the LP contribution cannot reproduce these two channels [24]. This demonstrates the importance of NLP contribution in the heavy quarkonium production. Since polarization is an important observable to explore the production mechanism of heavy quarkonium [25][26][27], it will be interesting to study polarized quarkonium production in the QCD factorization formalism. To this end, a set of polarized input FFs is required. Like the polarization-summed input FFs, we can use the NRQCD factorization to calculate these polarized input FFs. We will show in this paper that, if the polarization of heavy quarkonium is observed, many more NRQCD channels are needed in the factorization formula. Especially, offdiagonal channels are needed for the production of some heavy quarkonia even at the leading power in v. Conservation laws, angular momentum addition rules, and velocity scaling rules are very important to constrain various channels.
As in the case of the polarization-summed input FFs [19,20], the intermediate steps in the calculation of polarized input FFs involve ultraviolet (UV) and infrared (IR) divergences. Conventional dimensional regularization (CDR) is the most convenient tool to regularize all these divergences. However the application of CDR on the polarized input FFs is not straightforward, since definitions of polarized NRQCD LDMEs in arbitrary d-dimension are very non-trivial. In principle, one needs to consider the coupling of two d-dimensional heavy quark spinors, and L-S coupling if orbital angular momentum is nonzero, such as for the 3 P [1] J channels. In this paper, by requiring similar symmetries to be preserved when generalizing the polarized NRQCD LDMEs from 4 dimensions to d dimensions, we provide a simple scheme to separate the contributions with different J and |J z |. Since our scheme is quite general, it can also be applied on other circumstances.
The rest of this paper is organized as follows. In Section 2, within the NRQCD factorization formalism, we derive all nonvanishing NRQCD LDMEs for the production of polarized heavy quarkonium, by using conservation laws and angular momentum addition rules. Especially, we single out nonvanishing NRQCD LDMEs that contribute to polar angular distribution of heavy quarkonium production, which is of the most phenomenological interest. We also discuss how velocity scaling rules help to simplify problems in practice. Then in Section 3, we introduce a scheme to define these polarized NRQCD LDMEs in arbitrary d dimensions. In Section 4, we outline the application of NRQCD factorization on the single-parton and heavy quark-antiquark pair (QQ pair) FFs. Finally in Section 5, we briefly discuss our results for the polarized J/ψ production. Since the calculation is very similar to the polarization-summed case, which has already been explained in great details in Refs. [19,20], we do not repeat it in this paper, but providing all mathematical tools specific to the calculation of polarized input FFs in Appendix A. All our results for polarized FFs are listed in Appendixes B and C. As a few of these single parton FFs have been calculated in literature, we compare our results with them in Appendix B.5.

Heavy quarkonium polarization in NRQCD factorization
In this section, we start from experimental observables and find all allowed NRQCD channels for polarized heavy quarkonium production. In Subsection 2.1, the most general selection rules for nonvanishing NRQCD channels are obtained, by employing only the conservation laws and angular momentum addition rules. These selection rules are valid for any observables related to the heavy quarkonium polarization. Then in Subsection 2.2, we focus on the polar angular distribution and derive more selection rules for this specific observable. Finally in Subsection 2.3, we adopt the NRQCD power counting rules to find the relative importance of the survived channels for polar angular distribution.

General selection rules
Experimentally, the polarization of a produced heavy quarkonium is determined by measuring the angular distribution of its decay products in its rest frame. Take J/ψ as an example, the distribution of the decay products l + l − can be parameterized as (see Ref. [28] for a thorough discussion, and [29] for a recent review) dσ J/ψ(→l + l − ) dΩ ∝ 1 + λ θ cos 2 θ + λ ϕ sin 2 θ cos 2ϕ + λ θϕ sin 2θ cos ϕ + λ ⊥ ϕ sin 2 θ sin 2ϕ + λ ⊥ θϕ sin 2θ sin ϕ, where θ and ϕ are the polar angle and the azimuthal angle of l + respect to a chosen frame. Similarly, one can do the decomposition for χ cJ (see [30][31][32]) and other heavy quarkonium states. In the narrow-width approximation for heavy quarkonium, all the coefficients λ's in Therefore for the yield of H, only the diagonal entries in J and J z are needed.

Nonvanishing Channels for polar angular distribution
In this paper, we focus on polar angular distribution, which is of the most phenomenological interest. For J/ψ, it means that we consider only λ θ in Eq. (2.1), which can be separated from other λ's by integrating the azimuthal angle ϕ from 0 to 2π on both sides of Eq. (2.1).
Then we obtain . Thus we can use Eqs. (2.3) and (2.6) to find out all relevant QQ-pair channels (n,ñ). In addition to all diagonal channels with n =ñ, there are also many offdiagonal channels. For example, if n is S-wave, there are S-D mixing channels; if n is P -wave, there are P -P mixing channels as well as P -F mixing channels: J,Jz , Based on these results, we have a simpler factorization formalism for polar angular distribution. The time-reversal invariance of QCD implies 0|O J,−Jz and similarly for −ñ. Then, the Eq. (2.2) gives n,ñ +ρ 11) Table 1. Essential channels for various polarized heavy quarkonium production, with relative power-counting of each channel explicitly. λ = L, T, T T is the polarization of the heavy quarkonium. Q can be either charm or bottom.

Quarkonium
Essential channels

NRQCD power counting
The above general discussions are based only on symmetry of QCD. In principle, one needs all channels obtained there for any quarkonium production. In practice, however, different channels are not equally important. The relative importance in power of v for different channels is given by the velocity scaling rule [4].
To derive the scaling rule, one first expands the wave function of a physical heavy quarkonium H into different Fock states in NRQCD effective field theory. For a conventional heavy quarkonium, its dominant Fock state is a QQ-pair Fock state with definite quantum number (this quantum number is usually used to denote the quarkonium), while all other Fock states are suppressed by powers of v. The relative suppression can be estimated by color multipole expansion: If two Fock states are related by a color electric dipole transition (E1), the suppression is at O(v); While if the two Fock states are related by a color magnetic dipole transition (M1), the suppression is at O(v 3/2 ) [4]. Therefore, one has the expansion for the wave function of a heavy quarkonium H with quantum number 2S+1 L J , (2.14) In addition to suppression caused by the Fock states expansion, there is another suppression at O(v L+L ) for the QQ-pair channel (n,ñ) because of derivative operators. Based on these power counting rules, one can easily work out the relative importance in v of different channels for a definite heavy quarkonium production, e.g., the relative power counting of phenomenologically important channels for S-wave and P -wave quarkonium production is given in Table 1.
Three conclusions are in order: • First, a channel mixing between different waves is always suppressed by at least O(v 2 ). For example, the 1 S channel for any quarkonium state.
• Second, a channel mixing between different total angular momentum but the same orbital angular momentum can be expressed in terms of "diagonal" channel up to O(v 2 ) correction. For example, the mixing channels 3 P in Eq. (2.10) contributes to the production of H = J/ψ, ψ(nS) or Υ(ns). If we assume the heavy quark spin symmetry, the spin of the intermediate QQ-pair in the amplitude, S z , and in the complex-conjugate amplitude,S z , should equal to the spin of the hadron J H z . Further more, the corresponding mixing NRQCD LDMEs equal to each other due to heavy quark spin symmetry. As a result, we can define a "diagonal" channel 3,Sz P [8] , 3,Sz P [8] to study heavy quarkonium production, in which the spin of the QQ-pair is S z , and the orbital angular momentum is summed over. As heavy quark spin symmetry holds up to O(v 2 ) correction [4], this treatment is valid at the same precision order.
• Third, for any heavy quarkonium state, there is at least one S-wave channel which contributes at leading power in v; for any P -wave and higher-wave heavy quarkonium state, there is at least one P -wave channel which contributes at leading power in v; and so on. Thus lower-wave channels are usually more important.
In this consideration, we will calculate only S-wave and P -wave channels in this paper, including "diagonal" channels 3,Sz P [b] , 3,Sz P [b] reduced from P -P mixing channels in Eq. (2.10), and leave complete treatment of offdiagonal channels and higher-wave channels for future study. As it involves only diagonal channels in our treatment [up to O(v 2 ) corrections], the factorization formula in Eq. (2.11) is further simplified to where, similar to the notation in Eq. (2.11), λ = L, T, T T, · · · represent |J z | = 0, 1, 2, · · · , respectively. For example, and To relate our results in this paper to polarization-summed results, we sum over the polarization λ in Eq. (2.15), which gives is the usual polarization-summed NRQCD LDME. To obtain the second equal sign in Eq. (2.17), we have used the factor that the summation λ 0|O H λ [QQ(n λ )] |0 is independent of λ because of rotation invariance. Our calculations in this paper will be enough for the current purpose of studying Swave and P -wave heavy quarkonia production phenomenologically. In the next section, we provide a simple scheme to define all these operators in d dimensions.

Definition of polarized NRQCD LDMEs in d dimensions
In Refs. [19,20], we calculated both the single parton and heavy quark pair FFs for an polarization-summed heavy quarkonium. We used CDR to regularize all the UV and IR divergences. By doing this, we implicitly generalized the polarization-summed NRQCD LDMEs to arbitrary d dimensions. This generalization is simple, in the sense that with J z of the heavy quarkonium summed over, there is no special direction in the heavy quarkonium rest frame. In this case, one possible generalization can be achieved by extend the O(3) symmetry to the O(d − 1) symmetry [19,20], i.e., we demand all QQ-pair channels to be covariant under the action of O(d − 1) rotation group, similar to the 4 dimensions case.
The situation is more complicated with the polarization, where a specific direction z-axis needs to be specified. For example, in the helicity frame, z-axis is chosen to be along the moving direction of the heavy quarkonium in the Laboratory frame. To separate contributions with the same J but different |J z |, we need to know more details of the O(d − 1) rotation group. The situation is more severe with angular momentum couplings, such as the L-S coupling, which is exactly what we have in the NRQCD factorization.
In this section, we provide a simple scheme to separate the contributions with different J and |J z |, which is required in the calculation of λ θ . In our scheme, we only require that, first, the wave function of the heavy quark pair preserves all the symmetries about the z-axis when it is generalized to d dimensions, and second, the d-dimensional polarizationsummed NRQCD LDMEs in Ref. [19,20] are recovered when adding up the polarized ones. For S-wave and P -wave channels, the following rules can be applied, • Case 1: for 3 S 1 and 1 P 1 , the wave functions of the heavy quark pair with |J z | = 1 (or J z = 0) in its rest frame are antisymmetric (or symmetric) when flipping the direction of all axes exceptẑ axe. Note that, it is hard to define the rotation operation in d dimensions, but the flipping operation described here is always well-defined 2 .
• Case 2: for 3 P J , the wave function with J = 0 is a scalar, i.e. it has SO(d − 1) symmetry. Wave functions with J = 2 (or J = 1) are symmetric (or anti-symmetric) tensor in its orbital and spin indices and are constructed to be traceless.
• Case 3: for 3 P 2 , the wave functions with |J z | = 1 are separated from those with |J z | = 0, 2 by requiring that they are antisymmetric under flipping the direction of all axes except z, similar to Case 1. The wave functions with J z = 0 and |J z | = 2 are further separated by requiring the wave function with J z = 0 has SO(d−2) symmetry in the space perpendicular to z-axis.
• Case 4: for 3 P 1 , wave functions with |J z | = 1 are singled out by requiring that they are antisymmetric under flipping the direction of all axes except z, similar to Case 1. We define the symmetric part under this operation as J z = 0 3 .
Notice that in d dimensions, there is no unique way to group states into categories with different J and |J z |. Different grouping methods are equally good as long as they are consistent and give the correct decomposition at d = 4. They serve as different schemes in dimensional regularization.
As an example, we consider the polarized LDMEs for J/ψ production below, which is of great phenomenological interest, while giving definitions of all S-wave and P -wave NRQCD LDMEs in Appendix A. The most important [QQ(n)] channels in Eq. (4.2) for J/ψ production are 3 S J (see Table. 1). For 3 S [1] 1 channel, the wave function for polarization-summed QQ-pair in its rest frame can be chosen as a where ψ (χ) is the annihilation (creation) operator for a heavy quark (heavy antiquark) at the center of mass frame of the pair, base , and c j are the corresponding coordinates. In this Hilbert space with d = 4, the 3 × 3 matrix representation of the flipping operation described in Case 1 is where the summation of j ⊥ is suppressed, a H λ is the annihilation operator of heavy quarkonium H with polarization λ, and 1 (d−2)(2Nc) is a normalization factor. The corresponding NRQCD LDMEs for 3 S For J/ψ production, the analysis of the 3 P

[8]
J channel is greatly simplified if we assume the heavy quark spin symmetry, as explained at the end of Subsection 2.3. Since the spin of the produced J/ψ must be the same as the spin of the intermediate QQ-pair in the vanish at d = 4 and they have different parity with |Jz| = 1 wave functions, so it is natural to consider |Jz| = 2 wave functions as longitudinal. 4 The coupling of two spinors in d dimensions also have high-spin representations. On the one hand, because of heavy quark spin symmetry, their contribution to spin-0 and spin-1 representation is at least suppressed by v 3 [35]. On the other hand, one can alway ignore them by using a different dimensional regularization scheme, where spinors are in 4 dimensions.
Here we only consider the components at leading powers in v, relativistic corrections can be constructed by inserting covariant derivatives, as in Ref. [4].
leading power of v, the orbital angular momentum of the heavy quark pair can be effectively summed over. In this way, we can apply the same flipping operation as in Eq. (3.2), but this time only to the spin index. Then we obtain A complete treatment of the 3 P

NRQCD factorization on input FFs
The NRQCD factorization formula derived in Eq. (2.15) can also be applied to FFs at the initial scale. Specifically, input single parton FFs and QQ-pair FFs to polarized heavy quarkonium can be written in the factorized form 5 where µ 0 2m Q is the QCD factorization scale and µ Λ ∼ m Q is the NRQCD factorization scale. In the above equation, we single out the polarization parameters, with λ (λ ) denotes the polarization of the intermediate NRQCD QQ pair (observed heavy quarkonium H). In Eq. (4.1a), f = Q,Q, q,q, g, and the variable z is the light-cone momentum fraction taken by H; in Eq. (4.1b), κ = v [1,8] , a [1,8] or t [1,8] represents the heavy quark pair in spinor state vector, axial-vector or tensors, respectively, with the superscript labeling color-singlet (1) or color-octet (8). In Eq. (4.1b), we adopt the variables ζ 1 and ζ 2 defined in Ref. [19,20], with ζ 1 (ζ 2 ) labeling the relative momentum of the heavy quark pair in the amplitude (complex conjugate of the amplitude).
To calculate the short-distance coefficients (SDCs)d's in Eq. (4.1), one replaces the final heavy quarkonium state by an asymptotic heavy quark pair [QQ(n λ )] [QQ(n λ )] , (4.2b) where we suppress the factorization scales. In Eq. (4.2), the LHS can be calculated with perturbative QCD, while the LDMEs on the RHS can be calculated with perturbative NRQCD. By matching the LHS and the RHS, the SDCs can be obtained order by order in α s . In this paper, we follow the convention in Ref. [19,20] and expand the SDCs as where λ and λ represent the polarization of the NRQCD heavy quark pair and the observed heavy quarkonium, respectively. Now with the d-dimensional polarized NRQCD LDMEs and the projection operators derived from these LDMEs, we can apply CDR to regularize all the divergences in the calculation of the LHS of Eq. (4.2), with a particular polarized [QQ(n λ )] state. In the NLO calculation of the LHS, there exists (1) UV divergence, which are removed by QCD renormalization and the renormalization of composite operators; (2) rapidity divergence, which are canceled when adding up all Feynman diagrams; (3) Coulomb divergence and IR divergence, which are canceled between the LHS and the NRQCD LDMEs on the RHS, if NRQCD factorization is valid. The cancellation of the Coulomb divergence and IR divergence needs NLO calculation of polarized NRQCD LDMEs, which are also given in Appendix A.
These polarized SDCs and the polarization-summed SDCs are not independent. With the definitions of the polarized NRQCD LDMEs and the projection operators in Appendix A, one can obtain the relation [QQ(κ)]→[QQ(n λ )] for one polarization state (two polarization states) for 3 S 1 , 1 P 1 , and 3 P 1 ( 3 P 2 ) channel. The rest one can be calculated with Eq. (4.4).

Discussion and Summary
In this paper, we derived all nonvanishing channels for polarized heavy quarkonium production in NRQCD factorization. We then defined a scheme to generalize these polarized NRQCD LDMEs to arbitrary d dimensions. With these d-dimensional NRQCD LDMEs, we apply CDR to calculate the FFs for all polarized S-wave and P -wave heavy quarkonium with input scale µ 0 2m Q , to the NLO of α s . We find that all divergences are canceled and the SDCs are finite as expected. With our results, as well as the evolution kernels and the hard parts calculated in Refs. [8,18], the QCD factorization formalism can be used to give predictions to the polarizations of heavy quarkonia in hadron helicity frame, which have been measured in the Tevatron and the LHC.
In addition to the surprisingly large NLP contribution to the yield of J/ψ [23], our results in this paper further demonstrate the potential importance of the NLP to the J/ψ polarization. In Table 2, we show the LO contributions from different channels to the J/ψ polarization. The single-quark fragmentations are not important since they are suppressed at large z [19]. The contribution with a fragmenting QQ pair in a tensor state is suppressed by the hard part [18]. These suppressed channels are not showed in Table 2. For completeness, we also list the 1 S [8] 0 channel, which contributes to unpolarized J/ψ production. Table 2. Contributions of LO FFs to the J/ψ polarization. The labels "T", "L", and "Un" represent transversely polarized, longitudinally polarized, and unpolarized, respectively.
a [8] T Un Contrary to the LP contribution, the NLP mainly contributes to longitudinally polarized J/ψ. Recall that the NLP contributes both directly, by the NLP term in the factorized cross section, and indirectly, through the mixed kernel in the LP evolution equation (see Ref. [8] for more details). Especially, the indirect contribution increases the longitudinal component to the gluon fragmentation function when it is evolved from the input scale µ 0 2m Q to the hard scale µ ∼ p T . Although definite conclusion needs more studies, we believe our results in this paper will be very helpful to understand the polarization of heavy quarkonium production.

Acknowledgments
We thank E. Braaten for helpful discussions and G.T. Bodwin for useful communication regarding comparisons between our work and the results in Ref. [36]. HZ would like to thank the hospitality of the Peking University. This work was supported in part by the U.

A Polarized NRQCD LDMEs
This appendix is organized as follows. In Subsection A.1, the definitions of all essential d-dimensional polarized NRQCD LDMEs with S-wave and P -wave heavy quark pair are given. In Subsection A.2, we derive the projection operators for the calculation of SDCs. In Subsection A.3, we expand the polarized NRQCD LDMEs to NLO, which are necessary for the cancellation of IR divergences in the NLO calculation of SDCs.

A.1 NRQCD LDMEs
With the method expained in Section 3, we give our definitions of NRQCD four-fermion operators for polarized heavy quarkonium production in arbitrary dimension d.
where the summations of j ⊥ and k ⊥ run over all directions perpendicular toẑ-axis, and Subscripts T T , T , and L represent the non-relativistic QQ pair with |J z | = 2, |J z | = 1 and J z = 0, respectively. Polarization of the heavy quarkonium H is labelled by λ, which can be T T , T and L, depending on the specific state of H. For completeness, we also list the operators for states 3 P 0 and 1 S 0 , which are unpolarized. As mentioned at the end of Section 3, the 3 P

[8]
J channel is greatly simplified with heavy quark spin symmetry, where we only need the two NRQCD LDMEs in Eq. (3.4): The color-singlet operators can be obtained by removing the two T a 's and multiplying the factor 1/(2N c ) in their color-octet counterparts in Eqs. (A.1) and (A.3). With the definitions above, it is straightforward to check that by adding up operators with the same J but different |J z | and choose d = 4, the definitions of unpolarized NRQCD LDMEs in Ref. [4] can be retrieved 6 .

A.2 Projection Operators
With the definitions of d-dimensional polarized NRQCD LDMEs, we derive the projection operators in the same way as in the polarization-summed case.
α and β (α and β ) are the indices for the orbital angular momentum and spin of the heavy quark pair in the amplitude (complex conjugate of the amplitude), respectively. By adding up the projection operators with different |J z | but the same J, we can get the same results in Ref. [19,20]. Choosing d = 4 in the P NR 3 P 2,λ with λ = L, T, T T , we retrieve the results in Ref. [37].

A.3 Expand polarized LDMEs to NLO with perturbative NRQCD
To cancel the IR divergence from the NLO full QCD calculation, we need to expand the polarized NRQCD LDMEs to NLO, as for the polarization-summed case. We calculate the NLO NRQCD corrections for four-fermion operators using similar method as described in Ref. [38]. For our purpose, we only give the NLO correction of S-wave 4-fermion operators: 2,T T ) where ellipsis denotes terms that are irrelevant for our calculation.
To derive the expressions in Eqs. (A.6), we replace d by (4 − 2 ) only for the dimension of loop momentum in the calculation while keeping d's in other places untouched, and we then define a renormalization scheme for LDMEs to subtract all terms that are proportional to C U V , which is the same as C by replacing IR by U V . Our renormalization scheme is similar to the conventional MS scheme, but not exactly the same. The advantage of our renormalization scheme is that the calculated finite short-distance coefficients are unaltered if LDMEs are defined differently by a factor of 1 + O( ).

B Single-Parton Fragmentation Functions to Polarized Heavy Quarkonium
In this appendix we list the SDCs for all single-parton fragmentation functions to S-wave and P -wave polarized QQ-pair up to order O(α 2 s ) in Eq. (4.3a). The polarized FFs and unpolarized FFs are related by Eq. (4.4). Therefore for outgoing QQ with n+1 polarizations, we only give n polarized FFs below. The other one can then be calculated with Eq. (4.4) and the unpolarized FFs calculated in Ref. [19]. Those channels, which are zero for all polarizations, are not listed.
FFs in g → QQ( 3 P 2 ) channels in different polarization states have been calculated in Ref. [41] with cutoff regularization scheme. Comparing their results with ours in Eqs. (B.2), (B.3) and (B.4), we find that if we change log(µ 2 Λ /m 2 Q ) by log(µ 2 Λ /m 2 Q )+ 5/3 in our results, and take into account the different normalization of the LDMEs, we can reproduce their results in Ref. [41]. The 5/3 difference is due to different regularization scheme and was also realized in the calculation of polarization-summed FFs by authors of Ref. [42]. Therefore, our results of FFs in g → QQ( 3 P [1] 1 ) and g → QQ( 3 P [1] 2 ) channels in different polarization states are consistent with the results in Ref. [41].
Polarized FFs in q → QQ( 3 S 1,L ) channels have been reported in Ref. [36] when we are preparing our paper. 7 The authors of Ref. [36] present their results with NRQCD LDMEs defined slightly different from ours. After taking into account the difference of the LDMEs, our results in Eqs. (B.9), (B.10) and (B.15) are consistent with theirs.

C Double-Parton Fragmentation Functions to Polarized Heavy Quarkonium
In this appendix we list the SDCs for all heavy quark pair fragmentation functions to S-wave and P -wave polarized QQ-pair up to order O(α s ) in Eq. (4.3b). The polarized 7 Our results were published previously in Ref. [43], where there is a typo for q → QQ( 3 S [8] 1,L ) channel.
FFs and unpolarized FFs are related by Eq. (4.4). Therefore for outgoing QQ with n + 1 polarizations, we only give n polarized FFs below. The other one can then be calculated with Eq. (4.4) and the unpolarized FFs calculated in Ref. [19,20].Those channels, which are zero for all polarizations, are not listed.

C.1 ∆-functions
In our results below, we use the same definitions of ∆-functions as in Ref. [20]. We repeat these definitions below for readers' convenience.