Boosting perturbative QCD stability in quarkonium production

The aim of this paper is to introduce a general way to stabilize the perturbative QCD computations of heavy quarkonium production in the boosted or high-momentum transferring region with tree-level generators only. Such an approach is possible by properly taking into account the power-enhanced perturbative contributions in a soft and collinear safe manner without requiring any complete higher-order computations. The complicated NLO results for inclusive quarkonium hadroproduction can be well reproduced within our approach based on a tree-level generator HELAC-Onia. We have applied it to estimate the last missing leading-twist contribution from the spin-triplet color-singlet S-wave production at Oαs5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}_s^5\right) $$\end{document}, which is a NNLO term in the αs expansion for the quarkonium PT spectrum. We conclude that the missing NNLO contribution will not change the order of the magnitude of the short-distance coefficient. Such an approach is also quite appealing as it foresees broad applications in quarkonium-associated production processes, which are mostly absent of complete higher-order computations and fragmentation functions.


Introduction and motivations
As a class of the simplest hadrons, heavy quarkonium is usually viewed as the "hydrogen atom" in the strong-interaction theory QCD. While the knowledge of the nonperturbative aspect in QCD is still quite limited, heavy quarkonium provides a unique opportunity to probe the quark confinement in QCD by exploring the physics at the scale around the nonperturbative and perturbative boundary. The intrinsic scales of the heavy quark mass m Q and their binding energy m Q v 2 lie in the perturbative and nonperturbative regimes respectively, where v is the relative velocity between the heavy quark pair in the rest frame of the quarkonium. Due to the smallness of the relative velocity v (e.g. v 2 0.3 and v 2 0.1 for the charmonium and bottomonium), the relativistic QCD can be reorganized via the operator product expansion in the power counting of v. The effective theory was dubbed as non-relativistic QCD (NRQCD) [1]. The reformulation of QCD provides a factorization conjecture for calculating the rates of the quarkonium production and decay. In the case of the quarkonium H production, the (differential) cross section at leading-order (LO) in the QCD strong coupling constant α s can be schematically written as where n represents a Fock state, dσ(n) is a perturbatively calculable short-distance coefficient (SDC) with the heavy quark pair in the Fock state n and O H (n) is the vacuum expectation number of an operator O H (n). If the factorization formula eq. (1.1) holds, the nonperturbative long-distance matrix element (LDME) O H (n) is independent of quarkonium production process as well as the production environment. The universal LDMEs, which are analogous to the parton-distribution functions (PDFs) in the perturbative QCD factorization, are to be determined from a subset of the experimental data and to predict all of the rest experimental measurements. They have the probability explanations at LO, while LDMEs depend on the renormalization scheme and they are not physical objects. The prediction power of eq. (1.1) heavily relies on the perturbative convergences of v 2 and α s in dσ(H). The leading power counting of various Fock states up to O(v 7 ) for S-wave and P-wave quarkonia is listed in table 1 according to the NRQCD velocity scaling rule [1]. The convergence in v 2 can be improved by including the relativistic corrections. However, the prices to pay are that one has to introduce more nonperturbative LDMEs that can not be determined from the first principle, and the good relations like heavy-quark spin symmetry holding at LO in v will be violated too.
The most subtle part is the α s stability in the SDCs dσ(n), which is the main point to be discussed in this paper. For the high-transverse-momentum (P T ) quarkonium production at a high-energy hadron collider, it was found that 3 S [1] 1 receives a giant K factor from QCD corrections to its SDC a decade ago [2], which was understood by the fact that due to the quantum number conservation, there is a factor α s s ) (next-toleading order, NLO) compared to O(α 3 s ) (LO). This enhancement spoils the perturbative convergence in α s , shedding light on another possible enhancement from O(α 5 s ) (next-tonext-to-leading order, NNLO) corrections, while the accomplishment of the full NNLO calculation is even lacking today. The sole reason is the partonic cross sections dσ dP 2 T , before convoluting PDFs, are asymptotically scaling as  Table 2. The first α s orders needed in the SDCs for both LP and NLP in P T of various Fock states in their hadroproduction in order to achieve the LO and NLO QCD accuracies.
NNLO respectively. 1 Therefore, even with a full NNLO calculation at O(α 5 s ), the accuracy for the LP part of 1 has the leading P T behaviour as the jet, which means the LP channel is already opened at LO O(α 3 s ). In table 2, we have collected the first α s powers in order to achieve the LO and NLO QCD accuracies for various Fock states at both LP and NLP in P T .
Following this observation, the complete NLO result for 1 production is possible to be reproduced by the tree-level matrix element alone at O(α 4 s ) after introducing an ad hoc infrared cutoff. A first attempt was given in ref. [3] to introduce an invariant-mass cut on any final-final and initial-final massless parton pairs, which was called NLO . It can successfully reproduce the high-P T NLO calculation for 1 production. 2 The same infrared cut can be imposed in the phase-space integration of the O(α 5 s ) tree-level matrix element. Another giant K factor was observed compared to the NLO calculation at high P T , which may question on the extractions of color-octet LDMEs based on NLO calculations [6][7][8][9][10]. In contrast, a suspicion in ref. [11] on the size of O(α 5 s ) was given from their P T scaling reanalysis of the NNLO curves. Instead of the P T power enhancement, the observed giant K factor dσ NNLO dσ NLO is mainly attributed to the introduction of the infrared cutoff. Therefore, a reliable estimate of the size of O(α 5 s ) is still missing. It is necessary to clarify the situation before drawing a solid conclusion.
The aim of this paper is to introduce an infrared-safe method to cure the problematic giant K factors appearing in the SDC calculations in particular for high-P T quarkonium production without performing complete higher-order calculations. 3 In contrast to the NLO 1 Rigorously speaking, the associated production of 1 production, NLO cut was also applied to 1 hadroproduction in ref. [4]. NLO calculation is able to well reproduce the complete NLO result [5] in the double charmonium/bottomonium production. Its good performance may rely on the fact that like the single 1 is also NNLP in PT in the large transverse momentum region. 3 In the processes of elementary particle production, a few proposals to cure the giant K factors, which are mainly from logarithmic terms in perturbative calculations, are present [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. Unfortunately, none of them is straightforwardly applicable to the power-enhanced contributions in quarkonium production.

JHEP01(2019)112
calculations, the new method will not introduce the logarithmic dependence from the infrared cutoff. The estimate of the missing higher orders is to use the conventional renormalization and factorization scale variations. It is complemented with the fragmentation function approach, which requires the analytical calculations of different single-and doubleparton fragmentation functions for single and multiple quarkonium production. Another nontrivial task to use the fragmentation function approach is to solve the corresponding coupled evolution equations. It has been shown in ref. [29] that the fragmentation function approach without scale evolution can reproduce the spin-summed NLO cross sections of J at high P T , which shows the necessity of taking into account both the single-parton (at LP) and the double-parton (at NLP) fragmentation contributions. The factorization theorem for the single-inclusive quarkonium production cross sections in terms of single-and double-parton fragmentation functions was first proven in ref. [30] under the assumption of perturbative QCD factorization.
There are also other appealing reasons to introduce such a method. First of all, it can be used to stabilize the higher-order QCD corrections in quarkonium associated production processes, where most of them are absent of complete NLO calculations. The possible cancellations between S-wave and P-wave are guaranteed in our approach. For instance, in the double J/ψ at the LHC, it requires a NNLO calculation to have the full cancellations between S-wave and P-wave Fock states. As we will see later in this paper, the good reproduction of the NLO results both in the spin-summed and spin-dependent cross sections for single quarkonium production at high P T can serve as a fast way to the future phenomenology studies. In practice, the phenomenology from a complicated calculation scales as an inverse power of the computation time.
The outline of the remaining context is following. After introducing the remainders of P-wave counterterms in section 2, we will show that one can reproduce the NLO results for most of the Fock states (except 1 ) with fairly simple cuts based on tree-level matrix elements in section 3. These simple cuts are not sufficient to remove large logarithms introduced by the phase space cut parameters. Hence, a general infrared-safe method is introduced to obtain the giant K factors for all the Fock states relevant for J/ψ and χ cJ production in section 4. Finally, we draw our conclusions in section 5. An instruction on how to use HELAC-Onia [31,32] to perform the calculations done in this paper is given in appendix A. The appendix B contains supplemental figures.

Remainders of P-wave counterterms
It is well-known that the remaining infrared divergences in the SDC computations for the productions and decays of P-wave Fock states should be cancelled by the P-wave counterterms arising from the renormalization group running of S-wave LDMEs beyond LO in α s , which is analogous to the remaining collinear divergences absorbed by the PDF counterterms in a peturbative QCD calculation. The renormalization of NRQCD operators links the S-wave LDMEs with the P-wave LDMEs as shown in eq. (150) of ref. [33]. Such counterterms, after cancelling infrared divergences with the real and virtual matrix elements, will leave finite remainders proportional to the S-wave SDCs and P-wave LDMEs. The

JHEP01(2019)112
introduction of the P-wave counterterms is crucial especially in the case that the S-wave SDCs are much larger than the P-wave SDCs. In particular, the negative P-wave SDCs for heavy quarkonium hadroproduction at high P T could be attributed to these negative remainders. We have implemented the following finite remainders of P-wave counterterms: where m Q is the mass of the heavy quark and µ Λ is the NRQCD scale. In the following, we will set µ Λ = m Q as usually done in the complete NLO calculations. These remainders have already been implemented in the HELAC-Onia [31,32]. They are necessary ingredients to reproduce the NLO results, which we will show in the following two sections.

A first step towards NLO
From the discussion in the section 1, it is known that large NLO QCD corrections to the J/ψ production at a high-energy hadron collider are mainly due to the emergence of new P T power-enhanced fragmentation contributions. Hence, all S-and P-wave Fock states except 1 receive giant K factors from NLO QCD calculations. We first introduce the following basic phase space cuts in order to take into account the hard radiations without using virtual amplitudes. In real part at O(α 4 s ), exact 2 lightflavoured jets 4 satisfying P T (j) > P min T and |y(j)| < y max are required, which is denoted as dσ R 0 . The phase space integrations of Born dσ B (O(α 3 s )) and the remainders of the NRQCD P-wave counterterms dσ C (O(α 4 s )) are infrared safe with P T (onium) larger than J ) /(2J + 1) 1.16 GeV 3 9.03 · 10 −3 GeV 3 1.46 · 10 −2 GeV 3 3.43 · 10 −2 GeV 5 0.107 GeV 5 Table 3. The values of LDMEs used in the differential distributions of various Fock states.
a given positive value P min T (onium). We call the summed results of dσ B + dσ R 0 + dσ C as approximated NLO (aNLO).
In the following, we take P min T (onium) = 5 GeV, and light-flavoured jets are clustered with anti-k T algorithm [34] using radius R = 0.5 and |y(j)| < 5, P T (j) > P min T by Fast-Jet [35]. We will vary P min T from 3 GeV to 6 GeV as a way to estimate the infrared-cut dependence. We have shown the spin-summed double differential distributions for the cc Fock state 2 are displayed in figure 18 as our supplemental material. The complete NLO curves (denoting as NLO) from refs. [7,36] are also shown in order to have a comparison. The red-hatched bands represent the infrared cut variations P min T ∈ [3, 6] GeV, and the grey bands are the uncertainty from the independent variations of renormalization and factorization scales µ R , µ F around the central value µ 0 = P 2 T (onium) + 4m 2 c by a factor of 2. It is interesting to notice that the scale uncertainty in general captures the missing virtual and soft/collinear pieces. The agreements between NLO and aNLO are improved as P T (onium) increases. A similar behaviour can be observed for the spin-dependent differential cross sections shown in figure 2 for 0 are trivial. 5 We have utilized CTEQ6M PDF [38] to be consistent with the NLO results. For the reproducible purpose, the values of LDMEs for the distributions of the Fock states are listed in table 3.
Because the LP in P T for 3 S [8] 1 already exists at Born dσ B (i.e. O(α 3 s )) from the gluon fragmentation, it is expected that the scale uncertainty at LO would already give a reliable estimate of the missing NLO QCD corrections, which is indeed observed from the leftpanel of figure 3. In such a case, a request of 2 light-flavoured jets in the computation of dσ R 0 is insufficient to obtain an infrared-safe differential cross section. From the rightpanel of figure 3, the aNLO P T spectra are too hard compared to the complete NLO ones. The reason is because of the large logarithms arising from the very asymmetric dijet system P T (j 1 ) P T (j 2 ). Such a configuration is suppressed in other Fock states, because the leading fragmentation topologies require at least one light-flavoured parton along with the quarkonium direction at high P T . The weights of the asymmetric dijet events will be enhanced due to the unphysical logarithm log P T (j 1 ) P T (j 2 ) in the aNLO calculations of

Infrared-safe cutoffs
Let us assume a LO fragmentation process for a given Fock state O n is accompanying with k final massless partons: For a given observable, one needs to consider O n plus i recoiling partons. For example, in the case of the transverse-momentum distribution for a single quarkonium production (schematically depicted in figure 4), the quarkonium at least recoils against one parton at the lowest order (bar the zero transverse momentum bin). The soft-and collinear-safe calculations can be achieved based on pure tree-level matrix elements via the following conditions: 1. The number of jets is larger than i + 1 with the transverse momentum of jet P T (j) > P min T and the rapidity |y(j)| < y max . O n is also taken into account in the jet-clustering procedure. One should make sure that there is exactly one jet containing O n passing the above P T and rapidity cuts. Such a jet is called an onium-jet here.  3. If m ≥ 2, each pair of parton 1 and 2 inside the onium-jet should pass the following soft drop condition [39] min where p T,i is the transverse momentum of parton i and ∆R 12 = ∆φ 2 12 + ∆y 2 12 . The above cut already excludes the soft singularity as long as z cut > 0, while the requirement of the collinear safety is guaranteed by choosing β < 0. R 0 is the original jet radius, which is an order one number.
The condition eq. (4.2) in item 3 is chosen to kill the infrared unsafe configurations (a) and (b) given in figure 5. Either when partons 1 and 2 are close to be collinear ∆R 12 R 0 or if one parton is soft p T,2 p T,1 , eq. (4.2) cannot be fulfilled when z cut > 0, β < 0. In practice, the absolute value of β is at order one and z cut is at the order of v 2 .
If one goes to extra o radiations (i.e. O n plus i+k+o final light-flavoured QCD partons), one should impose the following additional cuts: 4. There are i + k + o − m partons outside the onium-jet. Each parton should form a single jet within P T (j) > P min T and |y(j)| < y max to avoid the collinear divergences. In order to get rid of large logarithms from infrared cuts arising from the soft largeangle radiations illustrated in the case (d) in figure 5, when i + k + o − m ≥ 2, it is necessary to impose an asymmetric cut on these parton jets The value of z cut,a should be a positive number smaller than 1 but not close to 0. It is important to vary its value in order to assess this cut dependence. 5. If k = 0, 7 and m > 0, each parton 1 in the onium jet j On should pass the soft cut

JHEP01(2019)112
where z 1 can be the energy fraction , the transverse momentum fraction p T,1 P T (j On ) or other similar fractions corresponding to z 1 → 0 when the parton 1 is soft. This condition is needed in order to kill the case (c) in figure 5, where all light-flavoured partons in the onium-jet j On can be soft and the condition eq. (4.2) is still satisfied. Similar to the value of z cut in eq. (4.2), the proper value of z cut,s should be O(v 2 ) as the effect of the soft radiations should be absorbed into the long-distance part of the quarkonium.
In fact, the combination of items 1-5 introduces a general infrared-safe method for any Fock state production if k = 0 is assumed at the beginning. 8 In other words, we do not need to pay a special attention to which kind of fragmentation process F is allowed for a given Fock state. We call such cuts as STOP cuts, where "STOP" is an acronym of "STabilize quarkOnium Production".
In the case of the P T spectrum of a quarkonium O n production at a hadron collider, i is equal to 1 and the LO process is O n plus one parton. For a real emission process O n plus o+1 partons with o > 0, we should impose the cuts listed in items 1-5 with k = 0, where the condition in item 2 is fulfilled automatically. Same as the previous section, we will denote the Born contribution at O(α 3 s ) as dσ B and the remainders of the P-wave counterterms at O(α 4 s ) as dσ C . dσ R STOP (dσ R 2 STOP ) stands for the contribution from O n plus two (three) partons within the STOP cuts.
1 ) fragmentation at LO. 8 When k = 0, the cut in item 2 will not be applied.

Reproducing NLO results
In this section, we will present the results up to NLO QCD corrections (i.e. O(α 4 s )). In order to differentiate our partial NLO calculations with the complete NLO results, we will denote our partial NLO calculations by imposing STOP cuts as "nLO", i.e. dσ nLO ≡ dσ B + dσ C + dσ R STOP . In the following, we will illustrate that the complete NLO results can be reproduced with the tree-level generators under the following setup of the STOP cuts: where m is the number of light-flavoured partons inside the onium jet. Jets are reconstructed with the anti-k T clustering algorithm using FastJet. Since there is no infrared divergence in the Born after imposing P T (onium) > 0 cut, the STOP cuts will not be applied to the Born and Born-like counterterm events.

Reproducing
After imposing the STOP cuts on 1 , we can reproduce the complete NLO curves within the theoretical uncertainties. They are shown in figure 6 and figure 7 for the spin-summed and spin-dependent differential cross sections respectively. In the left panel of figure 6 and the upper panels of figure 7, we estimate infrared cutoff dependence (the red-hatched bands) via the combined variations of P min T ∈ [3, 6] GeV and z cut,a ∈ [0.1, 0.7]. The grey-shadowed bands represent the scale uncertainties. Opposed to the aNLO results in the right panel of figure 3, it indeed shows that the STOP cuts improve the perturbative calculations, and the transverse-momentum dependence in dσ nLO dP T is the same as the NLO distributions dσ NLO dP T . It demonstrates that the large logarithmic dependence from the simple cuts in section 3 disappears after imposing the STOP cuts. The STOP-cut dependence (the red-hatched bands) is not reduced by increasing the P T of the quarkonium. It is expected since the LP contribution is already present at LO. The variations of the STOP cut variables only alter the fractions of hard radiations in the real matrix elements, which are not logarithmically enhanced. In fact, a careful tuning of STOP cut parameters can reproduce the NLO results at high precision. In the right panel of figure 6 and the lower panels of figure 7, we have calculated the 1 differential distributions after using z cut,a = 0.6 and z cut,s = 0.2 m . The comparisons to the full NLO calculations imply that the P T spectra of

Reproducing other Fock states
We are now in the position to check the calculations for the other Fock states with the STOP cuts. Like the case of the simple cuts in section 3, the general infrared-safe STOP cuts can reproduce the complete NLO results within theoretical uncertainties well. These can be better understood from the LO fragmentation function g → 1 S [8] 0 + g [40], which has the functional form where z onium is the momentum fraction of 1 S [8] 0 . The function peaks at z onium = 1. A finite value of z cut,s in the STOP cuts will remove a non-negligible fraction of radiations in the LP contributions. In fact, we have explicitly checked that if we set z cut,s = 10 −2 m instead of z cut,s = 0.1 m , the agreement between nLO and NLO results are significantly improved at large P T , which can be found in figure 10. In the spirit of the NRQCD factorization, the soft gluons from the heavy quark pair with the momentum fraction smaller than v 2 should be absorbed into the LDMEs as well as their energy evolutions, where v 2 is around 0.3 for the charmonium. Therefore, without taking into account the relativistic corrections, the resolution of NRQCD in describing the heavy quarkonium production should be not better than v 2 . Hence, it is not straightforward to judge which is a better choice between the two different values z cut,s = 0. We have compared the recent CMS measurement [41] to our nLO calculations (with and without STOP cut tuning on 1 in order to improve the visibility between the two theoretical bands. Without surprising, the CMS data agree very well with our nLO calculations, because nLO does a similarly good job as NLO.

Going beyond NLO
It is usually believed that the color-octet states for J/ψ hadroproduction will not receive giant K factors beyond NLO as the LP topologies in P T appear at NLO. On the other hand, the color-singlet Fock state   1 from NLO to NNLO might be possible in J/ψ production, though the NLO calculation shows that the 1 contribution to J/ψ hadroproduction seems to be negligible compared to the color-octet contributions. If it is true, the extractions of color-octet NRQCD LDMEs solely based on NLO calculations will be questionable. This is one of the reasons why the importance of color-octet contributions in J/ψ hadroproduction is still under debate. Although the accomplishment of NNLO calculations for 1 is still beyond state-of-the art, it was indeed suggested in ref. [3] that the partial calculation shows a giant K factor dσ NNLO dσ NLO . Later on, it was pointed out in ref. [11] that the giant K factor observed in ref. [3] is in fact due to the logarithmic enhancement induced by the infrared cutoff. Such a logarithm is expected to be absent in a full NNLO calculation because of the infrared safety.
We have the opportunity to clarify the situation with our infrared-safe STOP cut method. With the same setup done in section 4.2, we have performed the calculations for 1 plus three light-flavored jets production at O(α 5 s ) at the 13 TeV. The spin-summed P T differential distributions are shown in figure 12, where we have used "nnLO" and "nNLO" for the O(α n s ), n ≤ 4 parts being nLO and NLO cross sections respectively. In other words, we have used dσ nnLO ≡ dσ nLO + dσ R 2 STOP and dσ nNLO ≡ dσ NLO + dσ R 2 STOP . In the nNLO results, no theoretical uncertainties are taking into account from the NLO piece dσ NLO . In contrast to the finding made in ref. [3], we do not observe any giant K factor up to P T 100 GeV. In fact, the P T spectra of nnLO and nNLO are not harder than NLO ones. Such an observation can be explained if the coefficient of the LP P T part arising from the single-gluon fragmentation is much smaller than the coefficient of the NLP P T part and/or if the average momentum fraction of 1 taking from the original gluon is significantly smaller than 1. The calculation based on the gluon fragmentation function shows a similar behaviour, and the normalization of 1 is significantly smaller than the color-octet contributions [43]. In our calculation, the K factor dσ nnLO dσ NLO is ranging from 1 to 3 depending on the infrared cutoff choices. A similar conclusion can be drawn for the spindependent differential distributions from figure 13. We believe a complete NLO calculations of 1 plus two jets will help to reduce the remaining large infrared cutoff as well as the renormalization/factorization scale dependence.

Reassessing the charm fragmentation
So far, we have only considered the light-flavoured jet(s) accompanying with the quarkonium, which is usually thought to be dominant because the gluons are more often produced than the heavy quarks at high-energy hadron colliders. However, since the LP P T contribution from the charm quark fragmentation appears at O(α 4 s ), one should not overlook the associated production processes of a quarkonium plus a heavy quark pair. They were first studied in ref. [ To the best of our knowledge, the existing calculations only focus on the spin-summed differential cross sections, while we will also present the spin-dependent results in this section. In fact, one has to examine the relevance of these contributions if large cancellations between various Fock states happen.
In figure 14 and figure 15, we compared the tree-level 1 + cc (tagged as "cc") production with the nnLO calculations of 1 + cc contribution has a harder P T spectrum than the nnLO contribution. The former one exceeds the latter one above P T 55 GeV in the spin-summed case, while such a kind of crossover happens earlier for the spin transverse component dσ 11 dP T around P T 20 GeV. On the other hand, the charm quark associated contributions are orders of magnitude smaller than the light-flavoured jet contributions for

JHEP01(2019)112 5 Summary and outlooks
After implementing the remainders of P-wave counterterms in section 2, we have introduced a general infrared-safe method to estimate the giant K factors in quarkonium production in high P T region. As a proof of concept, we have validated our approach with the existing complete NLO QCD calculations of the Fock states in both spin-summed and spin-dependent cases. They are relevant for J/ψ and χ cJ hadroproduction up to O(v 7 ). Our approach only requires the tree-level amplitudes provided by HELAC-Onia. To the best of our knowledge, it is the first time to be able to reproduce the complete NLO spin-dependent results with tree-level amplitudes only. These spin-dependent results can be used to predict the polarization observables. We are also firstly able to obtain the spin-summed NLO results for production. It is believed to be at LP in P T scaling starting at this order, and is the last missing piece for the heavy quarkonium P T spectrum up to O(v 7 ). In contrast to the NNLO calculations based on the simple invariant-mass cuts [3], we do not observe the similar orders of magnitude enhancement compared to the NLO calculations, while an enhancement factor of 1 to 3 is still possible up to P T 100 GeV depending on the infrared cutoff choices. We believe the complete NLO calculations of 1 plus 2 jets will reduce this uncertainty. Finally, we have also calculated the charmonium plus a charm quark pair production, where the spin-dependent differential cross sections presented here are new. Their contributions to the inclusive P T distributions of charmonium are only relevant in the 1 channel. Our approach stabilizes the QCD corrections in the heavy quarkonium production rate calculations at high P T . It is quite appealing not only because it provides a fast way to perform the phenomenology studies of inclusive quarkonium production but also it can be used to improve the predictions in the associated quarkonium production processes. Together with the controlled perturbative SDCs, it is feasible to study various nonperturbative effects in the heavy quarkonium hadroproduction in an acceptable amount of computation time. Last but not least, with a similar method, we believe that we are able to promote the accuracy of both LP and NLP pieces to NLO level simultaneously with the full one-loop calculations.
Acknowledgments I thank Jean-Philippe Lansberg for useful discussions. This work is supported by the ILP Labex (ANR-11-IDEX-0004-02, ANR-10-LABX-63). The computations in this paper were performed with the help of the computing facilities at IPN Orsay.

A Calculations with HELAC-Onia
In this section, we will give an instruction on how to use HELAC-Onia to perform nLO and nnLO calculations. The implementations are available from version 2.3.6 and onwards,

A.1 Born and counterterms
The Born dσ B at O(α 3 s ) for S-wave Fock states can be achieved via the following commands: HO > define ppsi = cc~(1 s08 ) cc~(3 s11 ) cc~(3 s18 ) HO > decay ppsi > m + m -@ 1.0 HO > generate p g > ppsi j HO > generate g p > ppsi j HO > launch where we have always excluded the quark-quark initial states due to their very small parton luminosity from PDFs. No extra kinematical cuts are needed for the Born-like events. The color-octet P-wave Fock states J can be grouped together as we will always use the relation from heavy-quark spin symmetry O( HO > set generate_CT = T HO > decay cc~(3 p08 ) > m + m -@ 1.0 HO > generate p g > cc~(3 p08 ) j HO > generate g p > cc~(3 p08 ) j HO > launch The command sets generate_CT to be T in order to get the contributions from counterterms dσ C .
The color-singlet P-wave Fock states J will be calculated separately since they contribute to χ cJ , J = 0, 1, 2 respectively. The commands are: HO > set exp3pjQ = T

A.2 Real terms
One should apply the STOP cuts to the real terms at O(α n s ), n ≥ 4. It requires us to implement the following additional entries in user.inp: # STOP cuts use_stop_cut T # whether use stop cuts (following cuts will be ignored if it is F) stop_minptjet 3d0 # minimum jet pt cut (include onium in the jet clustering) stop_maxrapjet 5d0 # max jet rapidity, negative no such a cut stop_zcut 0.1d0 # zcut in the soft drop stop_beta -1d0 # beta in the soft drop (negative to make it collinear safe) stop_jet_alg -1 # 1: kt; 0: C/A; -1: anti-kt stop_jet_radius 1.0d0 # jet radius R stop_jet_dyn_radius -1d0 # if it is > 0, it will use dynamical jet radius R=max[stop_jet_radius,stop_jet_dyn_radius*M_onium/P_T,onium] stop_min_n_jet 2 # min number of jet (should be n hard jet + 1 onium jet), negative no such a cut stop_max_n_jet -1 # max number of jet, negative no such a cut stop_n_frag_gluon 0 # minimal number of final gluon in the LO fragmentation process stop_n_frag_quark 0 # minimal number of final bare quark in the LO fragmentation process stop_zsoftcut 0.1d0 # z_s,cut for the soft cut applied to (stop_n_frag_gluon+stop_n_frag_quark)==0 # It will be divided by number of partons inside the onium jet stop_zasymcut 0.1d0 # asymmetric cut for the light-flavoured jets if the number of light jets are >= 2 The

B Supplemental plots
We provided the supplemental plots in this appendix for the sake of completeness. The comparisons of spin-summed differential cross sections for the 5 Fock states 2 between aNLO and NLO are shown in figure 18, while the spin-dependent ones for

JHEP01(2019)112
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.