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 {\sc\small HELAC-Onia}. We have applied it to estimate the last missing leading-twist contribution from the spin-triplet color-singlet S-wave production at $\mathcal{O}(\alpha_s^5)$, which is a NNLO term in the $\alpha_s$ expansion for the quarkonium $P_T$ 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 an 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 leadingorder (LO) in the QCD strong coupling constant α s can be schematically written as (1.1) where n represents a Fock state, dσ(n) is a perturbatively calculable short-distance cofficient (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 enviroment. 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 hence are not physical objects beyond LO. 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.   [1].
The most subtle part is the α s stablity 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 1 hadroproduction is still at LO level, while the NLP piece is indeed NLO accurate. A NLO accuracy of the LP contribution can only be achieved with a next-to-NNLO calculation in α s for the SDC. The situation is slightly better though still similar for the other Fock states listed in Table 1. Like 1 S [1,8]  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 . 1 production is possible to reproduce 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 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 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 P T in the large transervese momentum region.
-3 -color-octet LDMEs based on NLO calculations [6,7,8,9,10]. On contrast, a suspect 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 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 double-parton 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 at high P T , which show the necessity of taking into account both the single-parton (at LP) and the double-parton (at NLP) fragmentation contributions. The factorization theorem for 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 stablize 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.

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 have finite remainders proportional to the S-wave SDCs and P-wave LDMEs. The 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 for large NLO QCD corrections to the J/ψ production at a high-energy hadron collider are mainly due to the appearance 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 light-flavoured 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 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 FastJet [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 6 cc Fock states Fig. 1 at the 13 TeV proton-proton collisions. The complete NLO curves (denoting as NLO) from Refs. [36,7] 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 Fig. 2 with the NLO curves from Refs. [8,37], where the nontrivial spin-density matrix elements are from 2 Fock states. 5 The negative differential cross sections for P-wave Fock states can be attributed to the negative finite remainders of the NRQCD P-wave counterterms dσ C . The SDCs presented in this paper are calculated by the program HELAC-Onia [31,32]. 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. 1 , 2 between our aNLO calculations and the complete NLO calculations. 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 ob--9 -served from the left-panel of Fig. 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 right-panel of Fig. 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 1 , which should be in principle cancelled by the virtual contributions because of the unitarity. Therefore, one must introduce a more general infrared-safe method to avoid these large logarithms, and at meantime one should keep the hard radiations from the real contributions.

A general infrared-safe method 4.1 Infrared-safe cutoffs
Let us assume a LO fragmentation process for a given Fock state O n is accompanying 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 -10 -production, 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.
2. In the onium-jet, there are at least k light-flavoured partons to fulfill the above fragmentation process. Let us say there are m light-flavoured partons inside the onium-jet with m ≥ k. 6 3. If m ≥ 2, each pair of parton 1 and 2 inside the onium-jet should pass the following soft drop condition [39] 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 safty is guaranteed by choosing β < 0. R 0 is the original jet radius, which is an order one number.
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 . In order to avoid large logarithms from infrared cuts, when i + k + o − m ≥ 2, it is necessary to impose an asymmetric cut on these parton jets min (P T (j 1 ), · · · , P T (j i+k+o−m )) max (P T (j 1 ), · · · , P T (j i+k+o−m )) > z cut,a . where z 1 can be energy fraction E 1 E(j On ) , transverse mass fraction , 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. 6 At LO one should have m = k since the configuration of i + k − m < i recoiling partons is zero by definition for the given observable when m > k. 7 For example, O n = QQ( 3 S -11 -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 begining. 8 In other words, we do not need to pay a special attention to which kind of fragmentation process F allowed for a given Fock state. We call such cuts as STOP cuts, where "STOP" is an acronym of "STablize 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 reminanders 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.

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 Fig. 4 and Fig. 5 for the spinsummed and spin-dependent differential cross sections respectively. In the left panel of Fig. 4 and the upper panels of Fig. 5, 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. In contrast to the aNLO results in the right panel of Fig. 3, it indeed shows that the STOP cuts improve the perturbative calculations, and the transverse-momentum 8 When k = 0, the cut in item 2 will not be applied.
-12 -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 tunning of STOP cut parameters can reproduce the NLO results at high precision. In the right panel of Fig. 4 and the lower panels of Fig. 5, we have calculated the

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 Fock states do not show LP behaviour at LO. The comparisons of nLO calculations to NLO calculations for spin-summed and spin-dependent differential cross sections are displayed in Fig. 6 and Fig. 7 respectively. With the scale variations shown by the grey bands, nLO results in general are able to successfully reproduce the NLO calculations in both cases. The only exception is the differential cross section of 1 S [8] 0 at very large P T (i.e. P T > 90 GeV). Such a discrepancy in 1 S [8] 0 can be better understood from the LO fragmentation function g → 1 S [8] 0 + g [40], which -14 -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 Fig. 8. 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, we think it is not straightforward to judge which is the better one between the two different choices z cut,s = 0.1 m and z cut,s = 10 −2 m . In fact, we believe z cut,s = 0.1 m is a compromising choice in order to avoid spoiling the perturbative convergence in the fixed-order calculations by a large logarithm log z cut,s .
We have compared the recent CMS measurement [41] to our nLO calculations (with and without STOP cut tunning on 1 ) for ψ(2S) production at 13 TeV LHC in Fig. 9, where the nonperturbative LDMEs are taken from Eqs.(2.17) and (2.18) in Ref. [42]. A factor 10 −1 has been multiplied to the nLO results with tunned 3 S [8] 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 -19 -calculations will be questionable. This is one of the reasons why the importance of color-octet contributions in J/ψ hadroproduction is still under debat. Although the accomplishment of NNLO calculations for 3 S [1] 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 safty.
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 ) in the 13 TeV proton-proton collisions. The spin-summed P T differential distributions are shown in Fig. 10, 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 part 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 only 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 spin-dependent differential distributions from Fig. 11. We believe a complete NLO calculations of

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 fragmention appears at O(α 4 s ), -21 -one should not overlook the associated production processes of a quarkonium plus a heavy quark pair. They were first studied in Ref. [44]  2 . 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 Figs. 12 and Fig. 14, we compared the tree-level 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 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 productions as clearly shown in Fig. 13 and Fig. 15 for the spin-summed and spindependent distributions.

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 -26 -with tree-level amplitudes only. These spin-denpendent results can be used to predict the polarization observables. We are also firstly able to obtain the spin-summed NLO results for 1 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 behaviour 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 channel. Our approach stablizes 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.
The color-octet P-wave Fock states 3 P [8] 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 -29 -counterterms dσ C .

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: -30 -# 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 HELAC-Onia commands to calculate the O(α 4 s ) real terms dσ R STOP are