1 Introduction↩︎

The theoretical interpretations of the analysed LHC data nowadays heavily rely on precision calculations of short-distance cross sections, as well as precise Monte Carlo event simulations in the context of both the Standard Model (SM) and its extensions. As of today, next-to-leading order (NLO) QCD calculations and their interface to general-purpose parton-shower Monte Carlo programmes have been automated, as seen, e.g., in the MadGraph5_aMC@NLO [1] framework, for elementary-particle production processes ¹ in the SM and in a large class of new physics models. NLO electroweak corrections have also been automated in recent years by several collaborations, such as MadGraph5_aMC@NLO [2] and Sherpa [3] along with external one-loop matrix element providers like Recola [4], OpenLoops [5] or GoSam [6]. Some low-particle-multiplicity processes, solved using customised methods, have been extended to next-to-NLO (NNLO) and even next-to-NNLO (N\(^3\)LO) accuracies. However, these significant theoretical developments are currently restricted to point-like elementary particles, and they cannot be directly applied to non-relativistic bound states like heavy quarkonia. This limitation can be roughly understood as the latter case intrinsically involving multiple scales and thus requiring simultaneous consideration in relativistic quantum field theories (QFT), such as QCD, and their non-relativistic low-energy effective field theories (EFT), e.g., non-relativistic QCD (NRQCD) [7]. This introduces additional conceptual and technical challenges on the theory side.

Therefore, theoretical progress in perturbative calculations for heavy quarkonium is far less advanced. Fairly speaking, automation has been achieved only for tree-level quarkonium processes (single quarkonium in MadOnia [8] and one or more quarkonia in HELAC-Onia [9], [10]). NLO and even higher-order calculations, in many cases, are indispensable not only for precision or accuracy but also for a qualitative understanding. Due to conservation laws at the quantum level, short-distance cross sections of quarkonium production often receive giant \(K\) factors ² from high-order radiative corrections (see, e.g., ref. [12] and references therein). This places the theoretical interpretations of measured quarkonium data on shaky ground if higher-order radiative corrections are not well under control.

The physics that we can learn from quarkonium is, however, no less interesting. In fact, quarkonia provide powerful and sometimes even unique tools that allow us to conduct rich particle and nuclear physics studies [13]. For instance, they can be used to determine the structures of free nucleons [14]–[25] and nuclei [26]–[29]. Due to their sequential binding energies, quarkonia are widely used as a thermometer of quark-gluon plasma produced in heavy-ion collisions [30], [31] to probe the hot-and-dense QCD. They signify the presence of a QCD phase transition by either disappearing [30], [31] or being abundantly produced, hinting at collective heavy-quark effects [32], [33]. A golden channel of searching for QCD instantons, arising from the non-trivial topological structure of QCD vacuum which is believed to be crucial in understanding quark confinement, was suggested to study charmonium decays [34], [35]. Quarkonia were also proposed as a good system to investigate the non-linear dynamics of QCD, also known as parton saturation [36], in addition to the well-known DGLAP and BFKL dynamics. They have been readily used to extract the fundamental SM parameters, e.g., the strong coupling constant \(\alpha_s\) [37], the Higgs-charm Yukawa coupling [38]–[40], the CKM matrix elements [41], [42], as well as the masses of the charm and bottom quarks [43]. Some exotic QCD hadrons, such as the fully-charmed tetraquark \(X(6900)\) [44]–[46] and the first-observed pentaquark states \(P_c^+(4380)\) and \(P_c^+(4450)\) [47], were also discovered in final states with quarkonia.

It has been well understood that the perturbative calculations at NLO and beyond in QFT encounter ultraviolet (UV) and infrared (soft and collinear) divergences in the intermediate steps. While the UV divergences are removed through the renormalisation procedure, handling infrared (IR) divergences is more intricate. First, IR divergences can only be cancelled for so-called IR-safe observables, thanks to the Kinoshita-Lee-Nauenberg (KLN) theorem [48], [49] and factorisation theorems/conjectures. Second, in a generic situation, phase-space integration must be carried out numerically using Monte Carlo importance sampling methods. However, this is hindered by the IR singularities present in real radiative corrections. Both IR subtraction and phase-space slicing approaches are employed to overcome such complications in the real-emission contributions, with IR subtraction methods known to outperform slicing approaches. Therefore, the former serves as the backbone of contemporary NLO automation codes. The two widely adopted NLO subtraction methods were originally proposed by Frixione, Kunszt and Signer (FKS) [50], [51], and Catani and Seymour [52], [53]. They are usually referred to as the FKS and dipole subtraction schemes, respectively. Both methods work well for both simple and complicated processes involving elementary particles. Initially devised for massless coloured particles, they have been generalised to include massive coloured particles [54]–[56]. Both schemes have been implemented in various public computer programmes [56]–[61].

On the other hand, regarding the problem of heavy quarkonium production in NRQCD, NLO calculations in the literature are almost exclusively carried out using the slicing methods [62], with only a few exceptions. The earliest exception ³ pertains to inclusive colour-singlet S-wave quarkonium production at hadron colliders [64], where the dipole counterterms for massless quarks and gluons suffice, as the colour-singlet S-wave quarkonium does not exhibit any IR singularities. The remaining exceptional NLO calculations [65], [66] employ the dipole formalism ⁴ developed for processes featuring an S- or P-wave quarkonium alongside massless quarks and gluons, as outlined in ref. [67], [68]. If one considers a process involving both a P-wave quarkonium and massive partons, such as the associated production processes of quarkonium and heavy quarks, new (yet unknown) dipole terms may be required. The aim of this paper is to incorporate heavy quarkonium into the FKS subtraction scheme. As we demonstrate later, the formalism is general enough to be applied to arbitrary processes involving a quarkonium and massless/massive partons. Therefore, the scope of the phenomenological applications using our formalism is anticipated to be broader than that of the dipole formalism derived in refs. [67], [68].

Since our ultimate goal is to automate NLO computations for quarkonium production processes within the MadGraph5_aMC@NLOframework, we will closely adhere to the notations and conventions of the original FKS formulation [50] and the MadFKSpaper [56]. For the sake of completeness and self-consistency, we will reproduce some known equations from the literature. We hope this will aid in improving the readability of the article, especially for readers who may not be familiar with the two aforementioned papers.

The remaining context of this paper is organised as follows. In section 2, we elucidate how to obtain short-distance cross sections for quarkonium production within NRQCD factorisation. This enables us to establish a few notations used throughout the paper. We derive the soft limit of the (squared) amplitudes in the real radiative corrections for a single quarkonium production in section 3. The local and integrated FKS subtraction counterterms are given in section 4. We perform a few cross-checks to ensure the validity of our formalism in section 5, and finally draw our conclusions in section 6. Appendix 7 presents the analytic expressions for eikonal tensor integrals that appear in the integrated counterterms. The universal IR poles of one-loop matrix elements can be found in appendix 8.

2 Quarkonium production in NRQCD factorisation↩︎

In NRQCD factorisation [7], the inclusive production of a heavy quarkonium factorises into the perturbative short-distance cross section and the non-perturbative long-distance matrix elements (LDMEs): \[d\sigma(AB\rightarrow H+X) = \sum_n \Big( \sum_{a,b,X}\int dx_a dx_b f_{a/A}(x_a) f_{b/B}(x_b) d\hat{\sigma}(ab \rightarrow Q\bar{Q}^\prime[n] + X) \Big) \braket{{\mathcal{O}}^H_n},\] where \(f_{a/A}\) and \(f_{b/B}\) are parton distribution functions (PDFs) of partons \(a\) and \(b\) in the initial hadrons \(A\) and \(B\). \(d\hat{\sigma}(ab \rightarrow Q\bar{Q}^\prime[n] + X)\) describes the short distance production of a \(Q\bar{Q}^\prime\) pair ⁵ in a specific colour irreducible representation \(C\), with spin \(S\) and orbital angular momentum state \(L\) denoted as \(n=\bigl.^{2S+1}L^{[C]}_J\) following the usual spectroscopic notation, and the LDME, \(\braket{{\mathcal{O}}^H_n}\), represents the hadronisation of the heavy quark pair into the physical quarkonium state \(H\). An important consequence of NRQCD factorisation is the prediction that the LDMEs do not depend on the details of the hard process, and their values can be extracted from experiments, lattice QCD calculations [69], [70] or potential models [71].

In principle, for a specific quarkonium, there is an infinite number of Fock states \(n\) and an infinite number of LDMEs \(\braket{{\mathcal{O}}^H_n}\) to be determined, which limits the prediction power. Thanks to the power counting rules in NRQCD, only a limited number of Fock states should be involved in the calculations up to a specific order of \(v\), where \(v\) (\(v\ll 1\)) is the relative velocity of the heavy quark pair \(Q\bar{Q}^\prime\). We express a given Fock state using the spectroscopic notation \(\bigl.^{2S+1}L^{[C]}_J\).

Our focus is to evaluate the perturbative short-distance coefficients, which can be determined from the amplitudes of \(Q\) and \(\bar{Q}^\prime\) production with necessary operations in order to constrain the heavy quark pair \(Q\bar{Q}^\prime\) into a specific quantum state \(n\). A convenient way to do so is by performing projections. Let us consider a general \(2\rightarrow n\) process involving only open \(Q\) and \(\bar{Q}^\prime\) quarks ⁶, denoted as \({\cal I}_1{\cal I}_2\rightarrow{\cal I}_3{\cal I}_4\cdots {\cal I}_{n+2}\), where the identity of the \(k\)-th particle is denoted by \({\cal I}_k\), and \({\cal I}_1=a, {\cal I}_2=b, {\cal I}_3=Q\), and \({\cal I}_4=\bar{Q}^\prime\). Following the same notation as ref. [56], we can write the process as \(r=\left({\cal I}_1,\ldots,{\cal I}_{n+2}\right)\). We denote the corresponding (tree-level) amplitude as \({\cal A}^{(n,0)}(r)\). Let us also define the amputated amplitude \(\Gamma^{(n,0)}(r)\) by removing the external wavefunctions of \(Q\) and \(\bar{Q}^\prime\) that will form a bound state \(Q\bar{Q}^\prime\), i.e., \[\begin{align} {\cal A}^{(n,0)}(r)&=&\bar{u}_{\lambda_Q}(k_Q)\Gamma^{(n,0)}(r)v_{\lambda_{\bar{Q}^\prime}}(k_{\bar{Q}^\prime}), \end{align}\] where \(u_\lambda\) and \(v_\lambda\) are Dirac spinors and \(\lambda_{Q/\bar{Q}^\prime}\) are helicities of \(Q\) and \(\bar{Q}^\prime\), respectively.

Since \(\mathbf{3}^{\count 0=0 \loop \ifnum\count 0>0 \advance\count 0 by -1 \prime\repeat}\otimes \settoheight{\irrepbarheight}{\mathbf{3}} \settowidth{\irrepwidth}{\mathbf{3}} \makebox[0pt][l]{\mathbf{3}} \rule[1.2\irrepbarheight]{\irrepwidth}{\irrepbarthickness}^{\count 0=0 \loop \ifnum\count 0>0 \advance\count 0 by -1 \prime\repeat}=\mathbf{1}^{\count 0=0 \loop \ifnum\count 0>0 \advance\count 0 by -1 \prime\repeat}\oplus \mathbf{8}^{\count 0=0 \loop \ifnum\count 0>0 \advance\count 0 by -1 \prime\repeat}\), we only have colour singlet \(C=1\) and colour octet \(C=8\) in the decomposition of \({\cal A}^{(n,0)}(r)\). The colour projectors are : \[\begin{align} \mathbb{P}_{C=1}&=&\frac{\delta_{c_4c_3}}{\sqrt{N_{c}}}, \nonumber \\ \mathbb{P}_{C=8} &=&\sqrt{2} t_{c_4c_3}^{c_{34}}, \end{align}\] where \(c_3,c_4\) are the colour indices of \(Q\) and \(\bar{Q}^\prime\), and \(t^{c_{34}}\) is the Gell-Mann matrix. In other words, we can define the following two amplitudes from \({\cal A}^{(n,0)}(r)\): \[\begin{align} {\cal A}^{(n,0)}_{\left\{[C=1]\right\}}(r)&=&\sum_{c_3,c_4}{\mathbb{P}_{C=1}{\cal A}^{(n,0)}(r)},\nonumber\\ {\cal A}^{(n,0)}_{\left\{[C=8]\right\}}(r)&=&\sum_{c_3,c_4}{\mathbb{P}_{C=8}{\cal A}^{(n,0)}(r)}, \end{align}\] where we have explicitly summed over the colour indices \(c_3, c_4\).

Similarly, we only have spin singlet \(S=0\) and spin triplet \(S=1\) for the \(Q\bar{Q}^\prime\) pair. The spin projectors for the heavy quark momenta \(k_{Q}^\mu=\frac{m_Q}{m_{Q}+m_{\bar{Q}^\prime}}K^\mu + q^\mu\) and \(k_{\bar{Q}^\prime}^\mu = \frac{m_{\bar{Q}^\prime}}{m_{Q}+m_{\bar{Q}^\prime}}K^\mu - q^\mu\) are given by \[\begin{align} \mathbb{P}_{S=0} &= & \frac{1}{2\sqrt{2m_{Q}m_{\bar{Q}^\prime}}} \bar{v}_{\lambda_{\bar{Q}^\prime}}(k_{\bar{Q}^\prime})\gamma_5 u_{\lambda_{Q}}(k_{Q}), \nonumber \\ \mathbb{P}_{S=1} &= & \frac{1}{2\sqrt{2 m_{Q}m_{\bar{Q}^\prime}}} \bar{v}_{\lambda_{\bar{Q}^\prime}}(k_{\bar{Q}^\prime})\cancel{\varepsilon}_{\lambda_s}^*(K)u_{\lambda_{Q}}(k_{Q}), \end{align}\] where \(\varepsilon_{\lambda_s}^*(K)\) is the polarisation vector for the spin-\(1\) \(Q\bar{Q}^\prime\) with its spin quantum number as \(\lambda_s=\pm1,0\). Here, \(K\) is the four-momentum of the \(Q\bar{Q}^\prime\) pair, \(q\) is the relative momentum between the two constituent heavy quarks, and \(m_{Q}, m_{\bar{Q}^\prime}\) are the heavy quark masses of \(Q\) and \(\bar{Q}^\prime\), respectively. For simplicitly, we just denote \(\tilde{\gamma}_0=\gamma_5\) and \(\tilde{\gamma}_1=\cancel{\varepsilon}_{\lambda_s}^*(K)\). Thus, the amplitudes can be further decomposed into two spin configurations \[\begin{align} {\cal A}^{(n,0)}_{\left\{S\right\}}(r)&=&\sum_{\lambda_{Q},\lambda_{\bar{Q}^\prime}}{\mathbb{P}_{S}{\cal A}^{(n,0)}(r)}\nonumber\\ &=&\sum_{\lambda_{Q},\lambda_{\bar{Q}^\prime}}{\frac{1}{2\sqrt{2m_{Q}m_{\bar{Q}^\prime}}} \bar{v}_{\lambda_{\bar{Q}^\prime}}(k_{\bar{Q}^\prime})\tilde{\gamma}_Su_{\lambda_{Q}}(k_{Q})\bar{u}_{\lambda_Q}(k_Q)\Gamma^{(n,0)}(r)v_{\lambda_{\bar{Q}^\prime}}(k_{\bar{Q}^\prime})}\nonumber\\ &=&\frac{1}{2\sqrt{2m_{Q}m_{\bar{Q}^\prime}}}{\rm Tr}_{\gamma}\left[\left(\cancel{k}_{\bar{Q}^\prime}-m_{\bar{Q}^\prime}\right)\tilde{\gamma}_S \left(\cancel{k}_{Q}+m_{Q}\right)\Gamma^{(n,0)}(r)\right]. \end{align}\] Together with both spin and colour configurations, the amplitude takes the form: \[\begin{align} {\cal A}^{(n,0)}_{\left\{[C],S\right\}}(r)&=&\sum_{\lambda_{Q},\lambda_{\bar{Q}^\prime}}{\mathbb{P}_{S}{\cal A}^{(n,0)}_{\left\{[C]\right\}}(r)}\nonumber\\ &=&\sum_{\lambda_{Q},\lambda_{\bar{Q}^\prime}}{\sum_{c_3,c_4}{\mathbb{P}_{S}\mathbb{P}_{[C]}{\cal A}^{(n,0)}(r)}}\nonumber\\ &=&\sum_{c_3,c_4}{\mathbb{P}_{[C]}\frac{1}{2\sqrt{2m_{Q}m_{\bar{Q}^\prime}}}{\rm Tr}_{\gamma}\left[\left(\cancel{k}_{\bar{Q}^\prime}-m_{\bar{Q}^\prime}\right)\tilde{\gamma}_S \left(\cancel{k}_{Q}+m_{Q}\right)\Gamma^{(n,0)}(r)\right]}. \end{align}\] Here, the spin projection operator commutes with the colour projection operator. Note that the trace \({\rm Tr}_{\gamma}\) in the Dirac spinor space does not necessarily represent a real trace that we need to compute. Its evaluation depends on how the other fermion lines are organised in the amputated amplitude \(\Gamma^{(n,0)}(r)\).

The non-relativistic nature, i.e., in the rest frame of \(Q\bar{Q}^\prime\), \(q\ll \sqrt{m_{Q}m_{\bar{Q}^\prime}}\), allows us to expand the amplitudes into the series of \(v\sim q/\sqrt{m_{Q}m_{\bar{Q}^\prime}}\ll 1\). This gives us the eigenfunctions of the orbital angular momentum operator. The projection on a state with orbital angular momentum \(L\) is obtained by differentiating \(L=0,1,\ldots\) (a la \(S,P,\ldots\) waves) times the spin-colour projected amplitude with respect to the relative momentum \(q\) of the heavy quarks in the \(Q\bar{Q}^\prime\) rest frame, and then setting \(q\rightarrow 0\). Considering only \(L=0,1\) states, which we are only interested in at this stage, the amplitude takes the form: \[\begin{align} {\cal A}^{(n,0)}_{\left\{[C], S, L\right\}}(r) & = & \left[\left(\varepsilon^{\mu,*}_{\lambda_l}(K)\frac{d}{dq^\mu}\right)^L {\cal A}^{(n,0)}_{[C],S}(r)\right]_{q=0}, \label{AmpQ} \end{align}\tag{1}\] where \(\varepsilon^{\mu,*}_{\lambda_l}(K)\) is the polarisation vector for \(L=1\) orbital angular momentum with \(\lambda_l=\pm1,0\). Since the spin projectors depend on the relative momentum \(q\), the orbital angular momentum expansion must be carried out after projecting onto the given spin configuration.

Finally, the total angular momentum \(J\) is uniquely determined by \(L\) or \(S\) unless \(L\neq 0\) and \(S\neq 0\). In our specific case of interest, this can only occur when \(L=1\) and \(S=1\). In the latter case, we know how to determine \(J=0,1,2\) and \(\lambda_j=-J,-J+1,\ldots, J-1,J\) from quantum mechanics, i.e., \[\begin{align} \varepsilon^{\mu \nu,*}_{J,\lambda_j}(K)&=&\sum_{\lambda_s,\lambda_l}{\langle J, \lambda_j| 1,\lambda_l; 1,\lambda_s\rangle \varepsilon^{\mu*}_{\lambda_l}(K)\varepsilon^{\nu*}_{\lambda_s}(K)},\label{eq:projLS2J} \end{align}\tag{2}\] where \(\langle J, \lambda_j| 1,\lambda_l; 1,\lambda_s\rangle\) is the Clebsch-Gordan coefficient. For \(J=0, 1\), the expressions are \[\begin{align} \varepsilon^{\mu \nu,*}_{0,0}(K)&=&\frac{1}{\sqrt{3}}\left(g^{\mu \nu}-\frac{K^\mu K^\nu}{K^2}\right),\nonumber\\ \varepsilon^{\mu \nu,*}_{1,\lambda_j}(K)&=&-\frac{i}{\sqrt{2}}\epsilon^{\mu \nu \alpha \beta }\frac{K_{\alpha}}{\sqrt{K^2}}\varepsilon^{*}_{\lambda_{j},\beta}(K), \end{align}\] where \(\epsilon^{\mu \nu \alpha \beta }\) is the Levi-Civita tensor. Thus, we obtain the amplitude \({\cal A}^{(n,0)}_{\left\{[C], S, L,J\right\}}(r)\) for a given quantum number \(\bigl.^{2S+1}L^{[C]}_J\). When the product \(L S=0\), we find \({\cal A}^{(n,0)}_{\left\{[C], S, L,J\right\}}(r)={\cal A}^{(n,0)}_{\left\{[C], S, L\right\}}(r)\). On the other hand, if \(L=S=1\), we have \[\begin{align} {\cal A}^{(n,0)}_{\left\{[C],1,1,J\right\}}(r)&=&\sum_{\lambda_s,\lambda_l}{\langle J, \lambda_j| 1,\lambda_l; 1,\lambda_s\rangle {\cal A}^{(n,0)}_{\left\{[C],1,1\right\}}(r)}\nonumber\\ &=&\frac{1}{2\sqrt{2m_Qm_{\bar{Q}^\prime}}}\sum_{c_3,c_4}\mathbb{P}_{[C]}\varepsilon_{J,\lambda_j}^{\mu \nu,*}(K)\nonumber\\ &&\times\left.\frac{d}{dq^\mu}{\rm Tr}_\gamma\left[\left(\cancel{k}_{\bar{Q}^\prime}-m_{\bar{Q}^\prime}\right)\gamma_\nu \left(\cancel{k}_{Q}+m_{Q}\right)\Gamma^{(n,0)}(r)\right]\right|_{q=0}\,. \end{align}\]

After all of the above preparations, we can now glue \({\cal I}_3=Q\) and \({\cal I}_4=\bar{Q}^\prime\) as a new single particle \({\cal I}_{3\oplus 4}=Q\bar{Q}^\prime[\bigl.^{2S+1}L^{[C]}_J]\) with four-momentum \(K^\mu\) and invariant mass \(\sqrt{K^2}=m_{Q\bar{Q}^\prime}=m_{Q}+m_{\bar{Q}^\prime}\) in the non-relativistic limit. The new process is denoted as \[\begin{align} \dot{r}&=&r^{3\oplus 4,4\backslash}= \left({\cal I}_1,{\cal I}_2,{\cal I}_{3\oplus 4},{\cal I}\backslash_{4},\ldots{\cal I}_{n+2}\right)\,. \label{QQproc} \end{align}\tag{3}\] The amplitude for the process \(\dot{r}=r^{3\oplus 4,4\backslash}\) is \[\begin{align} {\mathbb{A}}^{(n-1,0)}(\dot{r})&=&{\cal A}^{(n,0)}_{\left\{[C], S, L,J\right\}}(r). \end{align}\] Note that the final state symmetry must be applied at the level of \(\dot{r}\) not the initial \(r\), while the phase space integration should also be carried out at the level of \(\dot{r}\). The partonic cross section can be written as \[\begin{align} d\hat{\sigma}(\dot{r})&=&\frac{1}{{\cal N}(\dot{r})}\underbrace{\frac{1}{(2J+1)N_{[C]}}\frac{m_{Q}+m_{\bar{Q}^\prime}}{2m_{Q}m_{\bar{Q}^\prime}}}_{={\cal G}(\dot{r})}\left({\mathbb{M}}^{(n-1,0)}(\dot{r})J^{n_{L}^{(B)}}\right)d\phi_{n-1}(\dot{r}), \end{align}\] where the amplitude square is given by \[\begin{align} {\mathbb{M}}^{(n-1,0)}(\dot{r})&=&\frac{1}{2s}\frac{1}{\omega({\cal I}_1)\omega({\cal I}_2)} \mathop{\sum_{\rm colour}}_{\rm spin}\left|{\mathbb{A}}^{(n-1,0)}(\dot{r})\right|^2\nonumber\\ ={\cal M}^{(n,0)}_{\left\{[C],S, L,J\right\}}(r)&=&\frac{1}{2s}\frac{1}{\omega({\cal I}_1)\omega({\cal I}_2)} \mathop{\sum_{\rm colour}}_{\rm spin}\left|{\cal A}^{(n,0)}_{\left\{[C],S, L,J\right\}}(r)\right|^2, \end{align}\] \({\cal N}(\dot{r})\) is the final state symmetry factor, \(N_{[C=1]}=2N_{c}, N_{[C=8]}=N_{c}^2-1\) with \(N_c=3\) in QCD, and \(d\phi_{n-1}(\dot{r})\) is the \((n-1)\)-body phase space measure. The Mandelstam variable \(s=(k_1+k_2)^2=2k_1\cdot k_2\), and \(\omega({\cal I})\) is the product of spin and colour degrees of freedom for the particle \({\cal I}\). The condition with the measurement \(n_{L}^{(B)}\)-jet function \(J^{n_{L}^{(B)}}\), where \(n_{L}^{(B)}\) is the number of the light partons in the underlying Born, is sufficient to prevent the appearance of phase-space singularities in the Born-like quantities. Without losing generality, we can always assume these cuts to be equivalent to the request of having either \(n_{L}^{(B)}\) and \(n_{L}^{(B)}+1\) jets in the final state for an NLO computation. The same procedure can be iterated if we have more-than-one quarkonia.

3 Soft limit↩︎

At NLO, we have contributions coming from one-loop virtual corrections and real emissions besides Born. Due to the complexities introduced by bound states, it is necessary to derive new local and integrated FKS counterterms to handle the IR divergences in real contributions at NLO. We observe that, since the constituent quarks \(Q\) and \(\bar{Q}^\prime\) are massive, we can recycle the counterterms for collinear and soft-collinear origins. What we need to deal with beyond elementary particle production is the soft but non-collinear part. This includes two new components. The first involves the usual soft counterterms that locally cancel singularities of real emissions and their one-body phase space integrated counterparts. The second consists of additional integrated counterterms resulting from the renormalisation of LDMEs, analogous to the usual initial/final collinear counterterms but originating from a soft origin. Similar to the latter, which are necessary to cancel the remaining IR divergences protected by the collinear factorisation in perturbative QCD, the former are a consequence of the NRQCD factorisation formalism. In this section, we begin by considering the soft limit ⁷ of the quarkonium real emission amplitudes and their squares, following the procedure outlined before.

3.1 Soft limit at the amplitude level↩︎

Let us consider a \(2\rightarrow n+1\) real process denoted as \(r=({\cal I}_1,{\cal I}_2,{\cal I}_3,{\cal I}_4,\ldots,{\cal I}_{n+3})\). Following the approach in ref. [56], and without loss of generality, we can always reorder the final-state particles such that the final massless (anti-)quarks and gluons are \(3\leq i\leq n_{L}^{(R)}+2\) and the strongly-interacting particles are \(n_{I}\leq j \leq n_{L}^{(R)}+n_{H}+2\), where \(n_{I}=1, 2, 3\) for hadron-hadron, lepton-hadron, and lepton-lepton collisions. Here, \(n_{H}\) represents the number of massive coloured partons in the final state. Simultaneously, we have \({\cal I}_{n_{L}^{(R)}+n_{H}+1}=Q\) and \({\cal I}_{n_{L}^{(R)}+n_{H}+2}=\bar{Q}^\prime\). Consequently, we always have \(n_{H}\geq 2\). For simplicity, in the following discussion, we denote \(j_Q=n_L^{(R)}+n_H+1\) and \(j_{\bar{Q}^\prime}=n_L^{(R)}+n_H+2\). In the soft limit of \({\cal I}_i=g\) with \(3\leq i\leq n_{L}^{(R)}+2\), the real amplitude can be expressed in terms of the reduced Born amplitude multiplied with an additional eikonal factor. If the soft gluon with four-momentum \(k_i\) is emitted from an external leg \({\cal I}_j\) (\(j\neq i\)) with four-momentum \(k_j\), the real emission amplitude takes the form in the soft or eikonal approximation: \[\begin{align} \lim_{k_i \rightarrow 0}{{\cal A}^{(n+1,0)}(r)} = g_s\frac{k_j \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_j\cdot k_i} \vec{Q}({\cal I}_j) {\cal A}^{(n,0)}(r^{i\backslash}), \end{align}\] where the reduced process is \(r^{i\backslash}=\left({\cal I}_1,\ldots{\cal I}\backslash_i,\ldots{\cal I}_j, \ldots{\cal I}_{n+3}\right)\), and \(\vec{Q}({\cal I})\) represents the colour generator associated with the particle \({\cal I}\), and \(g_s=\sqrt{4\pi\alpha_s}\) is the strong coupling. Depending on the particle species \({\cal I}\), we have \[\begin{align} \label{SUN} \vec{Q}({\cal I}) = \{t^a\}_{a=1}^{8} ,\quad \{-t^{aT}\}_{a=1}^8, \quad \{T^a\}_{a=1}^8 \quad {\cal I}\in {\mathbf{3}^{\count 0=0 \loop \ifnum\count 0>0 \advance\count 0 by -1 \prime\repeat}, \settoheight{\irrepbarheight}{\mathbf{3}} \settowidth{\irrepwidth}{\mathbf{3}} \makebox[0pt][l]{\mathbf{3}} \rule[1.2\irrepbarheight]{\irrepwidth}{\irrepbarthickness}^{\count 0=0 \loop \ifnum\count 0>0 \advance\count 0 by -1 \prime\repeat},\mathbf{8}^{\count 0=0 \loop \ifnum\count 0>0 \advance\count 0 by -1 \prime\repeat}} \end{align}\tag{4}\] with \(t^a\) and \(T^a\) being the SU(3) generators in the fundamental and adjoint representations, respectively. We take the final quark or initial antiquark as \(\mathbf{3}^{\count 0=0 \loop \ifnum\count 0>0 \advance\count 0 by -1 \prime\repeat}\), while the initial quark and final antiquark are taken as $

^{=0 >0 by -1 }$. The matrix element of the adjoint representation is \(T^a_{bc}=-if_{abc}\) with \(f_{abc}\) being the anti-symmetric structure constants.

After performing the quantum number projection, as described in sect. 2, we obtain a similar eikonal decomposition as long as \(j\neq j_{Q}, j_{\bar{Q}^\prime}\), i.e., \[\begin{align} \lim_{k_i \rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[C],S, L,J\right\}}(r)} = g_s\frac{k_j \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_j\cdot k_i} \vec{Q}({\cal I}_j) {\cal A}^{(n,0)}_{\left\{[C],S, L,J\right\}}(r^{i\backslash}). \end{align}\] However, when \(j=j_{Q}\) or \(j=j_{\bar{Q}^\prime}\), we should pay special attention to it. It will be always convenient to sum the contributions of \(j=j_{Q}\) and \(j=j_{\bar{Q}^\prime}\) together, which we will adopt in the following. Now, let us consider the case of \(j=j_{Q},j_{\bar{Q}^\prime}\).

3.1.1 Colour projection↩︎

With the procedure outlined in sect. 2, the colour projected amplitudes are \[\begin{align} &&\lim_{k_i \rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[C=1]\right\}}(r)} \nonumber\\ &=& g_s\sum_{c_{j_{Q}},c_{j_{\bar{Q}^\prime}}}{\mathbb{P}_{[C=1]}\left[\frac{k_{j_{Q}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{Q}}\cdot k_i} \vec{Q}({\cal I}_{j_{Q}})+\frac{k_{j_{\bar{Q}^\prime}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{\bar{Q}^\prime}}\cdot k_i} \vec{Q}({\cal I}_{j_{\bar{Q}^\prime}}) \right] {\cal A}^{(n,0)}(r^{i\backslash})}\nonumber\\ &=&g_s\sum_{c_{j_{Q}},c_{j_{\bar{Q}^\prime}}}{\frac{\delta_{c_{j_{Q}}c_{j_{\bar{Q}^\prime}}}}{\sqrt{N_{c}}}\left[\frac{k_{j_{Q}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{Q}}\cdot k_i}t^a_{c_{j_{Q}}c_{j_{Q}}^\prime}{\cal A}^{(n,0)}_{c_{j_{Q}}^\prime c_{j_{\bar{Q}^\prime}}}(r^{i\backslash})-\frac{k_{j_{\bar{Q}^\prime}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{\bar{Q}^\prime}}\cdot k_i}t^a_{c_{j_{\bar{Q}^\prime}}^\prime c_{j_{\bar{Q}^\prime}}}{\cal A}^{(n,0)}_{c_{j_{Q}}c_{j_{\bar{Q}^\prime}}^\prime}(r^{i\backslash})\right]}\nonumber\\ &=&g_s\left[\frac{k_{j_{Q}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{Q}}\cdot k_i}-\frac{k_{j_{\bar{Q}^\prime}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{\bar{Q}^\prime}}\cdot k_i}\right]\frac{\delta_{a c_{j_{Q}j_{\bar{Q}^\prime}}^\prime}}{\sqrt{2N_{c}}}{\cal A}^{(n,0)}_{\left\{[C=8]\right\}}(r^{i\backslash}), \end{align}\] and \[\begin{align} &&\lim_{k_i \rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[C=8]\right\}}(r)}\nonumber\\ &=& g_s\sum_{c_{j_{Q}},c_{j_{\bar{Q}^\prime}}}{\mathbb{P}_{[C=8]}\left[\frac{k_{j_{Q}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{Q}}\cdot k_i} \vec{Q}({\cal I}_{j_{Q}})+\frac{k_{j_{\bar{Q}^\prime}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{\bar{Q}^\prime}}\cdot k_i} \vec{Q}({\cal I}_{j_{\bar{Q}^\prime}}) \right] {\cal A}^{(n,0)}(r^{i\backslash})}\nonumber\\ &=&g_s\!\!\sum_{c_{j_{Q}},c_{j_{\bar{Q}^\prime}}}{\sqrt{2}t^{c_{j_{Q}j_{\bar{Q}^\prime}}}_{c_{j_{\bar{Q}^\prime}}c_{j_{Q}}}\!\!\!\left[\frac{k_{j_{Q}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{Q}}\cdot k_i}t^a_{c_{j_{Q}}c_{j_{\bar{Q}^\prime}}^\prime}{\cal A}^{(n,0)}_{c_{j_{Q}}^\prime c_{j_{\bar{Q}^\prime}}}\!(r^{i\backslash})-\frac{k_{j_{\bar{Q}^\prime}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{\bar{Q}^\prime}}\cdot k_i}t^a_{c_{j_{\bar{Q}^\prime}}^\prime c_{j_{\bar{Q}^\prime}}}{\cal A}^{(n,0)}_{c_{j_{Q}}c_{j_{\bar{Q}^\prime}}^\prime}\!(r^{i\backslash})\right]}\nonumber\\ &=&g_s\left[\frac{k_{j_{Q}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{Q}}\cdot k_i}-\frac{k_{j_{\bar{Q}^\prime}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{\bar{Q}^\prime}}\cdot k_i}\right]\left[\frac{\delta_{a c_{j_{Q}j_{\bar{Q}^\prime}}}}{\sqrt{2N_{c}}}{\cal A}^{(n,0)}_{\left\{[C=1]\right\}}(r^{i\backslash})+\frac{1}{2}d_{ac_{j_{Q}j_{\bar{Q}^\prime}}c_{j_{Q}j_{\bar{Q}^\prime}}^\prime}{\cal A}^{(n,0)}_{\left\{[C=8]\right\}}(r^{i\backslash})\right]\nonumber\\ &&+g_s\left[\frac{k_{j_{Q}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{Q}}\cdot k_i}+\frac{k_{j_{\bar{Q}^\prime}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{\bar{Q}^\prime}}\cdot k_i}\right]\left(-\frac{i}{2}f_{a c_{j_{Q}j_{\bar{Q}^\prime}}c_{j_{Q}j_{\bar{Q}^\prime}}^\prime}\right){\cal A}^{(n,0)}_{\left\{[C=8]\right\}}(r^{i\backslash}), \end{align}\] where \(d_{abc}\)’s are the symmetric structure constants. Moreover, we have used the following relation of Gell-Mann matrices \[\begin{align} (t^at^b)_{kl} = \displaystyle\frac{\delta_{ab}}{2N_{c}}\delta_{kl} + \displaystyle\frac{1}{2}(d_{abc} + if_{abc})t^c_{kl} \end{align}\] and have assumed the colour index for \(Q\bar{Q}^\prime\) in the real process \(r\) (the reduced Born process \(r^{i\backslash}\)) to be \(c_{j_{Q}j_{\bar{Q}^\prime}}\) (\(c_{j_{Q}j_{\bar{Q}^\prime}}^\prime\)). We can put the two equations into the compact matrix form by using the colour nonet \(\mathbf{9}^{\count 0=0 \loop \ifnum\count 0>0 \advance\count 0 by -1 \prime\repeat}\) index \(\boldsymbol{b}=0,1,2,\ldots,8\), where \(\boldsymbol{b}=0\) corresponds to \({\cal A}^{(n+1,0)}_{\left\{[C=1]\right\}}(r)\) and \({\cal A}^{(n,0)}_{\left\{[C=1]\right\}}(r^{i\backslash})\), while \(\boldsymbol{b}=1,\ldots,8\) are \({\cal A}^{(n+1,0)}_{\left\{[C=8]\right\}}(r)\) and \({\cal A}^{(n,0)}_{\left\{[C=8]\right\}}(r^{i\backslash})\) with the colour index of the \(Q\bar{Q}^\prime\) as \(\boldsymbol{b}\). In other words, with the emitters being \(Q\) and \(\bar{Q}^\prime\), we have \[\begin{align} \lim_{k_i \rightarrow 0}{\left(\begin{array}{c}{\cal A}^{(n+1,0)}_{\left\{[C=1]\right\}}(r) \\ {\cal A}^{(n+1,0)}_{\left\{[C=8]\right\},\boldsymbol{b}=1}(r)\\ \vdots\\ {\cal A}^{(n+1,0)}_{\left\{[C=8]\right\},\boldsymbol{b}=8}(r) \end{array}\right)}&=&g_s\left\{\left[\frac{k_{j_{Q}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{Q}}\cdot k_i}-\frac{k_{j_{\bar{Q}^\prime}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{\bar{Q}^\prime}}\cdot k_i}\right]\vec{Q}_1(Q\bar{Q}^\prime)\right.\nonumber\\ &&\left.+\left[\frac{k_{j_{Q}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{Q}}\cdot k_i}+\frac{k_{j_{\bar{Q}^\prime}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{\bar{Q}^\prime}}\cdot k_i}\right]\vec{Q}_2(Q\bar{Q}^\prime)\right\}\nonumber\\ &&\times\left(\begin{array}{c}{\cal A}^{(n,0)}_{\left\{[C=1]\right\}}(r^{i\backslash})\\ {\cal A}^{(n,0)}_{\left\{[C=8]\right\},\boldsymbol{b}=1}(r^{i\backslash})\\ \vdots\\ {\cal A}^{(n,0)}_{\left\{[C=8]\right\},\boldsymbol{b}=8}(r^{i\backslash}) \end{array}\right), \end{align}\] where the colour generators are \[\begin{align} \vec{Q}_1(Q\bar{Q}^\prime)&=&\left(\begin{array}{cccc} 0 & \frac{\delta_{a1}}{\sqrt{2N_{c}}} & \cdots & \frac{\delta_{a8}}{\sqrt{2N_{c}}} \\ \frac{\delta_{a1}}{\sqrt{2N_{c}}} & \ddots & & \reflectbox{\ddots} \\ \vdots & & \frac{D^a}{2} & \\ \frac{\delta_{a8}}{\sqrt{2N_{c}}} & \reflectbox{\ddots} & & \ddots\\ \end{array}\right),\\ \vec{Q}_2(Q\bar{Q}^\prime)&=&\left(\begin{array}{cccc} 0 & 0 & \cdots & 0 \\ 0 & \ddots & & \reflectbox{\ddots} \\ \vdots & & \frac{T^a}{2} & \\ 0 & \reflectbox{\ddots} & & \ddots \\ \end{array}\right), \end{align}\] with the elements of the matrix \(D^a\) as \(D^a_{bc}=d_{abc}\).

3.1.2 Spin and orbital angular momentum projections↩︎

The next step is to perform the projection for the spin and the orbital angular momentum via \[\begin{align} &&\lim_{k_i \rightarrow 0}{\left(\begin{array}{c}{\cal A}^{(n+1,0)}_{\left\{[C=1],S,L\right\}}(r) \\ {\cal A}^{(n+1,0)}_{\left\{[C=8],S,L\right\},\boldsymbol{b}=1}(r)\\ \vdots\\ {\cal A}^{(n+1,0)}_{\left\{[C=8],S,L\right\},\boldsymbol{b}=8}(r) \end{array}\right)}\nonumber\\ &=&\frac{g_s}{2\sqrt{2m_{Q}m_{\bar{Q}^\prime}}}\left\{\left(\varepsilon_{\lambda_l}^{\mu,*}(K)\frac{d}{dq^\mu}\right)^L\left[\left(\frac{k_{j_{Q}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{Q}}\cdot k_i}-\frac{k_{j_{\bar{Q}^\prime}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{\bar{Q}^\prime}}\cdot k_i}\right)\vec{Q}_1(Q\bar{Q}^\prime)\right.\right.\nonumber\\ &&\left.+\left(\frac{k_{j_{Q}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{Q}}\cdot k_i}+\frac{k_{j_{\bar{Q}^\prime}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{\bar{Q}^\prime}}\cdot k_i}\right)\vec{Q}_2(Q\bar{Q}^\prime)\right]\nonumber\\ &&\times \left.\sum_{c_{j_{Q}},c_{j_{\bar{Q}^\prime}}}{\left(\begin{array}{c} \frac{\delta_{c_{j_{Q}}c_{j_{\bar{Q}^\prime}}}}{\sqrt{N_{c}}}\\ \sqrt{2}t^{1}_{c_{j_{\bar{Q}^\prime}}c_{j_{Q}}}\\ \vdots \\ \sqrt{2}t^{8}_{c_{j_{\bar{Q}^\prime}}c_{j_{Q}}}\\\end{array}\right){\rm Tr}_{\gamma}\left(\left(\cancel{k}_{\bar{Q}^\prime}-m_{\bar{Q}^\prime}\right)\tilde{\gamma}_S \left(\cancel{k}_{Q}+m_{Q}\right)\Gamma^{(n,0)}(r^{i\backslash})\right)}\right\}_{q=0}. \end{align}\] When \(L=0\), we just set the relative momentum \(q\) to be zero. The coefficient of \(\vec{Q}_1(Q\bar{Q}^\prime)\) vanishes. This implies the following two consequences:

• For the colour singlet \(C=1\) with \(L=0\) (S-wave), regardless of the value of the spin \(S\), there are no soft divergences. Thus, we can treat a colour-singlet S-wave state as any other elementary colour-singlet particle, such as \(Z\) and \(H\) bosons, from the IR perspective.

• For the colour octet \(C=8\) S-wave states, we have the following soft limit relation \[\begin{align} \lim_{k_i \rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[8],S,0\right\}}(r)}&=&g_s\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}\vec{Q}(Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[8]}_J]){\cal A}^{(n,0)}_{\left\{[8],S,0\right\}}(r^{i\backslash}),\label{eq:softOctetS} \end{align}\tag{5}\] where now \(Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[8]}_J] \in \mathbf{8}^{\count 0=0 \loop \ifnum\count 0>0 \advance\count 0 by -1 \prime\repeat}\) and \(\vec{Q}(Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[8]}_J])=\left\{T^a\right\}_{a=1}^{8}\). It means that, in the soft limit, a colour-octet S-wave state behaves as if \(Q\bar{Q}^\prime\) were an elementary colour-octet particle, akin to a sgluon or a Kaluza-Klein massive gluon.

For P-wave (\(L=1\)) Fock states, we have \[\begin{align} \lim_{k_i \rightarrow 0}{\left(\begin{array}{c}{\cal A}^{(n+1,0)}_{\left\{[C=1],S,1\right\}}(r) \\ {\cal A}^{(n+1,0)}_{\left\{[C=8],S,1\right\},\boldsymbol{b}=1}(r)\\ \vdots\\ {\cal A}^{(n+1,0)}_{\left\{[C=8],S,1\right\},\boldsymbol{b}=8}(r) \end{array}\right)}&=&g_s\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}2\vec{Q}_2(Q\bar{Q}^\prime)\left(\begin{array}{c}{\cal A}^{(n,0)}_{\left\{[C=1],S,1\right\}}(r^{i\backslash}) \\ {\cal A}^{(n,0)}_{\left\{[C=8],S,1\right\},\boldsymbol{b}=1}(r^{i\backslash})\\ \vdots\\ {\cal A}^{(n,0)}_{\left\{[C=8],S,1\right\},\boldsymbol{b}=8}(r^{i\backslash}) \end{array}\right)\nonumber\\ &&+g_s\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]\nonumber\\ &&\times\!\left[\left(2+\frac{m_{Q}}{m_{\bar{Q}^\prime}}+\frac{m_{\bar{Q}^\prime}}{m_{Q}}\right)\vec{Q}_{1}(Q\bar{Q}^\prime)\!+\!\left(\frac{m_{\bar{Q}^\prime}}{m_{Q}}\!-\!\frac{m_{Q}}{m_{\bar{Q}^\prime}}\right)\vec{Q}_{2}(Q\bar{Q}^\prime)\right]\nonumber\\ &&\times\left(\begin{array}{c}{\cal A}^{(n,0)}_{\left\{[C=1],S,0\right\}}(r^{i\backslash}) \\ {\cal A}^{(n,0)}_{\left\{[C=8],S,0\right\},\boldsymbol{b}=1}(r^{i\backslash})\\ \vdots\\ {\cal A}^{(n,0)}_{\left\{[C=8],S,0\right\},\boldsymbol{b}=8}(r^{i\backslash}) \end{array}\right), \end{align}\] because of \[\begin{align} \varepsilon_{\lambda_l}^{\mu,*}(K)\frac{d}{dq^\mu}\frac{k_{j_{Q}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{Q}}\cdot k_i}&=&\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{Q}}\cdot k_i}-\frac{k_{j_{Q}} \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(k_{j_{Q}}\cdot k_i\right)^2}\nonumber\\ &\overset{q=0}{=}&\frac{m_{Q}+m_{\bar{Q}^\prime}}{m_{Q}}\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right],\nonumber\\ \varepsilon_{\lambda_l}^{\mu,*}(K)\frac{d}{dq^\mu}\frac{k_{j_{\bar{Q}^\prime}} \cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{\bar{Q}^\prime}}\cdot k_i}&=&-\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{k_{j_{\bar{Q}^\prime}}\cdot k_i}+\frac{k_{j_{\bar{Q}^\prime}} \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(k_{j_{\bar{Q}^\prime}}\cdot k_i\right)^2}\nonumber\\ &\overset{q=0}{=}&-\frac{m_{Q}+m_{\bar{Q}^\prime}}{m_{\bar{Q}^\prime}}\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right].\nonumber\\ \end{align}\] This essentially means that

• For the colour-singlet (\(C=1\)) P-wave states, we get \[\begin{align} \lim_{k_i \rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[C=1],S,1\right\}}(r)}&=&g_s\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]\nonumber\\ &&\times\left(2+\frac{m_{Q}}{m_{\bar{Q}^\prime}}+\frac{m_{\bar{Q}^\prime}}{m_{Q}}\right)\frac{\delta_{a c_{j_{Q}j_{\bar{Q}^\prime}}^\prime}}{\sqrt{2N_{c}}}{\cal A}^{(n,0)}_{\left\{[C=8],S,0\right\}}(r^{i\backslash}).\label{eq:softSingletP1} \end{align}\tag{6}\] The object \({\cal I}_{j_{Q}\oplus j_{\bar{Q}^\prime}\oplus i}\) may be interpreted as a new colour-octet S-wave particle (with colour index \(c_{j_{Q}j_{\bar{Q}^\prime}}^\prime=a\)) with a non-standard eikonal factor.

• For the colour-octet (\(C=8\)) P-wave states, the soft limit yields a more complicated expression \[\begin{align} \lim_{k_i \rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[C=8],S,1\right\}}(r)}&=&g_s\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}\vec{Q}(Q\bar{Q}^\prime[\bigl.^{2S+1}P^{[8]}_J]){\cal A}^{(n,0)}_{\left\{[C=8],S,1\right\}}(r^{i\backslash})\nonumber\\ &&+g_s\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]\nonumber\\ &&\times\left[\left(2+\frac{m_{Q}}{m_{\bar{Q}^\prime}}+\frac{m_{\bar{Q}^\prime}}{m_{Q}}\right)\left(\frac{\delta_{a c_{j_{Q}j_{\bar{Q}^\prime}}}}{\sqrt{2N_{c}}}{\cal A}^{(n,0)}_{\left\{[C=1],S,0\right\}}(r^{i\backslash})\right.\right.\nonumber\\ &&\left.+\frac{d_{ac_{j_{Q}j_{\bar{Q}^\prime}}c_{j_{Q}j_{\bar{Q}^\prime}}^\prime}}{2}{\cal A}^{(n,0)}_{\left\{[C=8],S,0\right\}}(r^{i\backslash})\right)\nonumber\\ &&\left.+\frac{1}{2}\left(\frac{m_{\bar{Q}^\prime}}{m_{Q}}-\frac{m_{Q}}{m_{\bar{Q}^\prime}}\right)\vec{Q}(Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[8]}_J]){\cal A}^{(n,0)}_{\left\{[C=8],S,0\right\}}(r^{i\backslash})\right].\nonumber\\ \label{eq:softOctetP1} \end{align}\tag{7}\]

If we define the following effective colour generators \[\begin{align} \vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[18]})&=&\left\{\left(2+\frac{m_{Q}}{m_{\bar{Q}^\prime}}+\frac{m_{\bar{Q}^\prime}}{m_{Q}}\right)\frac{\Delta^{a}}{\sqrt{2N_{c}}}\right\}_{a=1}^{8},\nonumber\\ \vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})&=&\left\{\left(2+\frac{m_{Q}}{m_{\bar{Q}^\prime}}+\frac{m_{\bar{Q}^\prime}}{m_{Q}}\right)\frac{\Delta^{aT}}{\sqrt{2N_{c}}}\right\}_{a=1}^{8},\nonumber\\ \vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})&=&\left\{\left(2+\frac{m_{Q}}{m_{\bar{Q}^\prime}}+\frac{m_{\bar{Q}^\prime}}{m_{Q}}\right)\frac{D^a}{2}+\left(\frac{m_{\bar{Q}^\prime}}{m_{Q}}-\frac{m_{Q}}{m_{\bar{Q}^\prime}}\right)\frac{T^a}{2}\right\}_{a=1}^{8}, \end{align}\] where \(\Delta^a_{\circ b}=\delta_{ab}\) and \(\circ\) represents no colour index for a colour-singlet state, the soft-limit expressions can be written in a more compact form \[\begin{align} \lim_{k_i \rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[C=1],S,1\right\}}(r)}&=&g_s\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]\nonumber\\ &&\times\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[18]}){\cal A}^{(n,0)}_{\left\{[C=8],S,0\right\}}(r^{i\backslash}),\tag{8}\\ \lim_{k_i \rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[C=8],S,1\right\}}(r)}&=&g_s\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}\vec{Q}(Q\bar{Q}^\prime[\bigl.^{2S+1}P^{[8]}_J]){\cal A}^{(n,0)}_{\left\{[C=8],S,1\right\}}(r^{i\backslash})\nonumber\\ &&+g_s\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]\nonumber\\ &&\times\left[\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]}){\cal A}^{(n,0)}_{\left\{[C=1],S,0\right\}}(r^{i\backslash})+\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]}){\cal A}^{(n,0)}_{\left\{[C=8],S,0\right\}}(r^{i\backslash})\right].\nonumber\\ \tag{9} \end{align}\] Here, we have defined the effective Casimir constants through the following relations \[\begin{align} C_{{\rm eff}} (Q\bar{Q}^\prime_{[18]})&=&\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[18]})\!\cdot\!\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[18]})=\left(2+\frac{m_{Q}}{m_{\bar{Q}^\prime}}+\frac{m_{\bar{Q}^\prime}}{m_{Q}}\right)^2\frac{1}{2N_{c}}\,,\nonumber\\ C_{{\rm eff}} (Q\bar{Q}^\prime_{[81]})&=&\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})\!\cdot\!\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})=\left(2+\frac{m_{Q}}{m_{\bar{Q}^\prime}}+\frac{m_{\bar{Q}^\prime}}{m_{Q}}\right)^2C_{F}\,, \nonumber\\ C_{{\rm eff}} (Q\bar{Q}^\prime_{[88]})&=&\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})\!\cdot\!\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})=\left(2+\frac{m_{Q}}{m_{\bar{Q}^\prime}}+\frac{m_{\bar{Q}^\prime}}{m_{Q}}\right)\nonumber\\ &&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\times\left(-\frac{2}{N_{c}}+\frac{N_{c}^2-2}{2N_{c}}\left(\frac{m_{\bar{Q}^\prime}}{m_{Q}}+\frac{m_{Q}}{m_{\bar{Q}^\prime}}\right)\right).\nonumber\\ \end{align}\] In addition, we also have \[\begin{align} \vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})\!\cdot\!\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})&=&\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})\!\cdot\!\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})=0,\nonumber\\ \vec{Q}(Q\bar{Q}^\prime[\bigl.^{2S+1}P^{[8]}_J])\!\cdot\!\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})&=&\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})\!\cdot\!\vec{Q}(Q\bar{Q}^\prime[\bigl.^{2S+1}P^{[8]}_J])=0\,,\nonumber\\ \vec{Q}(Q\bar{Q}^\prime[\bigl.^{2S+1}P^{[8]}_J])\!\cdot\!\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})&=&\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})\!\cdot\!\vec{Q}(Q\bar{Q}^\prime[\bigl.^{2S+1}P^{[8]}_J])=\frac{1}{2}\left(\frac{m_{\bar{Q}^\prime}}{m_{Q}}-\frac{m_{Q}}{m_{\bar{Q}^\prime}}\right)C_{A}\,.\nonumber\\ \end{align}\]

Finally, in order to get a given total angular momentum \(J\) when \(L=S=1\), we need to use eq.@eq:eq:projLS2J at the amplitude level: \[\begin{align} \lim_{k_i \rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[C=1],1,1,J\right\}}(r)}&=&\sum_{\lambda_l,\lambda_s}{\langle J, \lambda_j|1,\lambda_l; 1, \lambda_s\rangle g_s\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]}\nonumber\\ &&\times\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[18]}){\cal A}^{(n,0)}_{\left\{[C=8],1,0\right\}}(r^{i\backslash}),\label{eq:softSingletP3} \end{align}\tag{10}\] and \[\begin{align} \lim_{k_i \rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[C=8],1,1,J\right\}}(r)}&=&g_s\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}\vec{Q}(Q\bar{Q}^\prime[\bigl.^{3}P^{[8]}_J]){\cal A}^{(n,0)}_{\left\{[C=8],1,1,J\right\}}(r^{i\backslash})\nonumber\\ &&+\sum_{\lambda_l,\lambda_s}{\langle J,\lambda_j|1,\lambda_l; 1,\lambda_s\rangle g_s\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]}\nonumber\\ &&\times\left[\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]}){\cal A}^{(n,0)}_{\left\{[C=1],S,0\right\}}(r^{i\backslash})+\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]}){\cal A}^{(n,0)}_{\left\{[C=8],S,0\right\}}(r^{i\backslash})\right].\nonumber\\ &&\label{eq:softOctetP3} \end{align}\tag{11}\] To rewrite this in a unified compact form, we can further define the following effective identity (called “flavour") operators \[\begin{align} \vec{F}_{{\rm eff}}(Q\bar{Q}^\prime_{[18]})&:&Q\bar{Q}^\prime[\bigl.^{2S+1}P^{[1]}_J]\rightarrow Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[8]}_{S}],\nonumber\\ \vec{F}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})&:&Q\bar{Q}^\prime[\bigl.^{2S+1}P^{[8]}_J]\rightarrow Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[1]}_{S}],\nonumber\\ \vec{F}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})&:&Q\bar{Q}^\prime[\bigl.^{2S+1}P^{[8]}_J]\rightarrow Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[8]}_{S}]. \end{align}\] The eikonal current operator for an elementary coloured particle \({\cal I}_j\) is \[\begin{align} \vec{J}({\cal I}_j)&=&\frac{k_j\cdot \varepsilon^*_{\lambda_i}(k_i)}{k_j\cdot k_i}\vec{Q}({\cal I}_j),\label{eq:eikonalcurrent1} \end{align}\tag{12}\] while for bound states it is \[\begin{align} \vec{J}(Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[1]}_{S}])&=&0,\nonumber\\ \vec{J}(Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[8]}_{S}])&=&\frac{K\cdot \varepsilon^*_{\lambda_i}(k_i)}{K\cdot k_i}\vec{Q}(Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[8]}_{S}]),\nonumber\\ \vec{J}(Q\bar{Q}^\prime[\bigl.^{1}P^{[1]}_{1}])&=&\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[18]})\vec{F}_{{\rm eff}}(Q\bar{Q}^\prime_{[18]}),\nonumber\\ \vec{J}(Q\bar{Q}^\prime[\bigl.^{1}P^{[8]}_{1}])&=&\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}\vec{Q}(Q\bar{Q}^\prime[\bigl.^{2S+1}P^{[8]}_J])\nonumber\\ &&+\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]\nonumber\\ &&\times\left[\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})\vec{F}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})+\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})\vec{F}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})\right],\nonumber\\ \vec{J}(Q\bar{Q}^\prime[\bigl.^{3}P^{[1]}_{J}])&=&\sum_{\lambda_l,\lambda_s}{\langle J,\lambda_j|1,\lambda_1;1,\lambda_s\rangle\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]}\nonumber\\ &&\times\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[18]})\vec{F}_{{\rm eff}}(Q\bar{Q}^\prime_{[18]}),\nonumber\\ \vec{J}(Q\bar{Q}^\prime[\bigl.^{3}P^{[8]}_{J}])&=&\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}\vec{Q}(Q\bar{Q}^\prime[\bigl.^{2S+1}P^{[8]}_J])\nonumber\\ &&+\sum_{\lambda_l,\lambda_s}{\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]}\nonumber\\ &&\times\langle J,\lambda_j|1,\lambda_l;1,\lambda_s\rangle\!\!\left[\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})\vec{F}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})\!+\!\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})\vec{F}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})\right].\nonumber\\ \label{eq:eikonalcurrent2} \end{align}\tag{13}\] In such a case, the amplitude in the soft limit of \(k_i\rightarrow 0\) can be rewritten as \[\begin{align} \lim_{k_i\rightarrow 0}{{\mathbb{A}}^{(n,0)}(\dot{r})}&=&g_s\mathop{\sum_{j=n_{I}}}_{j\neq i}^{n_{L}^{(R)}+n_{H}+1}{\vec{J}({\cal I}_j){\mathbb{A}}^{(n-1,0)}(\dot{r}^{i\backslash})},\label{eq:softamp} \end{align}\tag{14}\] where we have glued the \(Q\bar{Q}^\prime\) pair into a single particle in processes with a dot, i.e., \[\begin{align} \dot{r}&=&r^{j_{Q}\oplus j_{\bar{Q}^\prime}, j_{\bar{Q}^\prime}\backslash}=\left({\cal I}_1,\ldots,{\cal I}_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{2S+1}L^{[C]}_J],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{2S+1}L^{[C]}_J],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right). \end{align}\] Note that \(\dot{r}\) is essentially equivalent to eq.@eq:QQproc except that we have reordered the final particles, and the total number of external parton legs is \(n+3\) instead of \(n+2\).

3.2 Soft limit of real matrix elements↩︎

We can now examine the soft limit of the real amplitude square by considering the \(\ell\)-loop amplitude \({\cal A}^{(n,\ell)}(r)\) for a generic \(2\rightarrow n\) process \(r\). To facilitate our later discussion, we introduce the following amplitude squares: \[\begin{align} {\mathbb{M}}^{(n,0)}(\dot{r})&=&{\cal M}^{(n+1,0)}_{\left\{[C],S,L,J\right\}}(r)=\frac{1}{2s}\frac{1}{\omega({\cal I}_1)\omega({\cal I}_2)} \mathop{\sum_{\rm colour}}_{\rm spin}\left|{\cal A}^{(n+1,0)}_{\left\{[C],S,L,J\right\}}(r)\right|^2, \tag{15} \\ {\mathbb{M}}^{(n,0)}_J(\dot{r})^{\mu\nu}&=&\frac{1}{2s}\frac{1}{\omega({\cal I}_1)\omega({\cal I}_2)}\nonumber\\ &&\times\mathop{\sum_{\rm colour}}_{\rm spin}\left[\left(\sum_{\lambda_l,\lambda_s}{\langle J,\lambda_j| 1,\lambda_l;1,\lambda_s\rangle\varepsilon^{\mu,*}_{\lambda_l}(K){\mathbb{A}}^{(n-1,0)}(\dot{r})}\right)\right.\nonumber\\ &&\left.\times\left(\sum_{\lambda_l^\prime,\lambda_s^\prime}{\langle J,\lambda_j| 1,\lambda_l^\prime;1,\lambda_s^\prime\rangle\varepsilon^{\nu,*}_{\lambda_l^\prime}(K){\mathbb{A}}^{(n-1,0)}(\dot{r})}\right)^\star\right], \tag{16} \\ {\mathbb{M}}^{(n-1,0)}(\dot{r})&=&{\cal M}^{(n,0)}_{\left\{[C],S,L,J\right\}}(r)=\frac{1}{2s}\frac{1}{\omega({\cal I}_1)\omega({\cal I}_2)} \mathop{\sum_{\rm colour}}_{\rm spin}\left|{\cal A}^{(n,0)}_{\left\{[C],S,L,J\right\}}(r)\right|^2, \tag{17} \\ {\cal M}^{(n,0)}_{\left\{[C],S,L,J\right\},kl}(r)&=&-\frac{1}{2s} \frac{2-\delta_{kl}}{\omega({\cal I}_1)\omega({\cal I}_2)}\nonumber\\ &&\times\mathop{\sum_{\rm colour}}_{\rm spin} {\cal A}^{(n,0)}_{\left\{[C],S,L,J\right\}}(r) \vec{Q}({\cal I}_k)\!\cdot\!\vec{Q}({\cal I}_l) {{\cal A}^{(n,0)}_{\left\{[C],S,L,J\right\}}(r)}^{\star},\phantom{aa} \tag{18} \\ {\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r})&=&-\frac{1}{2s} \frac{2-\delta_{kl}}{\omega({\cal I}_1)\omega({\cal I}_2)}\nonumber\\ &&\times\mathop{\sum_{\rm colour}}_{\rm spin} {\mathbb{A}}^{(n-1,0)}(\dot{r}) \vec{Q}({\cal I}_k)\!\cdot\!\vec{Q}({\cal I}_l) {{\mathbb{A}}^{(n-1,0)}(\dot{r})}^{\star},\phantom{aa} \tag{19} \\ {\mathbb{M}}^{(n-1,0)}_{k[C_1C_2]}(\dot{r}_1,\dot{r}_2)^{\mu}&=&-\frac{1}{2s} \frac{1}{\omega({\cal I}_1)\omega({\cal I}_2)}\nonumber\\ &&\mathop{\sum_{\rm colour}}_{\rm spin} 2\Re{\left\{\varepsilon^{*,\mu}_{\lambda_l}(K){\mathbb{A}}^{(n-1,0)}(\dot{r}_1)^{\star}\vec{Q}({\cal I}_k)\!\cdot\!\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[C_1C_2]}) {{\mathbb{A}}^{(n-1,0)}(\dot{r}_2)}\right\}},\nonumber\\[-10pt] \tag{20} \\ {\mathbb{M}}^{(n-1,0)}_{J,k[C_1C_2]}(\dot{r}_1,\dot{r}_2)^{\mu}&=&-\frac{1}{2s} \frac{1}{\omega({\cal I}_1)\omega({\cal I}_2)}\nonumber\\ &&\mathop{\sum_{\rm colour}}_{\rm spin} 2\Re{\left\{{\mathbb{A}}^{(n-1,0)}(\dot{r}_1)^{\star}\vec{Q}({\cal I}_k)\!\cdot\!\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[C_1C_2]})\right.}\nonumber\\ &&\left.\times\sum_{\lambda_l,\lambda_s}{\langle J,\lambda_j|1,\lambda_l;1,\lambda_s\rangle\varepsilon^{*,\mu}_{\lambda_l}(K){\mathbb{A}}^{(n-1,0)}(\dot{r}_2)}\right\},\phantom{aa} \tag{21} \\ {\cal M}^{(n,1)}(r)&=&\frac{1}{2s}\frac{1}{\omega({\cal I}_1)\omega({\cal I}_2)} \mathop{\sum_{\rm colour}}_{\rm spin} 2\Re\left\{{\cal A}^{(n,0)}(r){{\cal A}^{(n,1)}(r)}^{\star}\right\}. \tag{22} \end{align}\] Once again, in these equations, \(s=(k_1+k_2)^2=2k_1\cdot k_2\), and \(\omega({\cal I})\) represents the product of spin and colour degrees of freedom for particle \({\cal I}\). In \(d=4-2\epsilon\) dimensions, the average factors are given by \(\omega(q)=\omega(\bar{q})=2N_{c}=6\) and \(\omega(g)=2(1-\epsilon)(N_{c}^2-1)=16(1-\epsilon)\). These average factors are utilised to fully specify the divergent part of the one-loop contribution within the conventional dimensional regularisation (CDR) scheme. It is important to note that while these factors include \(\epsilon\) for completeness, in numerical calculations, \(d=4\) tree-level amplitudes are typically evaluated, and the \(\epsilon\) dependence is dropped.

In general, we can write down the soft limit of the amplitude square in the form \[\begin{align} \lim_{k_i\rightarrow 0}{{\mathbb{M}}^{(n,0)}(\dot{r})}&=&g_s^2\mathop{\sum_{k,l=n_{I}}}_{k,l\neq i, k\leq l}^{n_{L}^{(R)}+n_{H}+1}{{\mathsf M}^{(n-1,0)}_{{\rm soft}, kl}(\dot{r}^{i\backslash})}.\label{eq:softampsq0} \end{align}\tag{23}\] It is well known that if both \({\cal I}_k\) and \({\cal I}_l\) are elementary particles, we have \[\begin{align} {\mathsf M}^{(n-1,0)}_{{\rm soft}, kl}(\dot{r}^{i\backslash})&=&\frac{k_k\!\cdot\!k_l}{k_k\!\cdot\!k_ik_l\!\cdot\!k_i}{\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}).\label{eq:softampsq} \end{align}\tag{24}\] In the following, we will try to derive \({\mathsf M}^{(n-1,0)}_{{\rm soft}, kl}(\dot{r}^{i\backslash})\) when at least \({\cal I}_k\) or \({\cal I}_l\) is a quarkonium state.

3.2.1 Colour-singlet S-wave state↩︎

We have derived the soft limit of colour-singlet S-wave states at the amplitude level in sect. 3.1 as \[\begin{align} \lim_{k_i\rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[1],S,0,J=S\right\}}(r)}&=&\mathop{\sum_{j=n_{I}}}_{j\neq i}^{n_{L}^{(R)}+n_{H}}{g_s\frac{k_j\cdot \varepsilon_{\lambda_i}^*(k_i)}{k_j\cdot k_i}\vec{Q}({\cal I}_j){\cal A}^{(n,0)}_{\left\{[1],S,0,J=S\right\}}(r^{i\backslash})}. \end{align}\] It follows that the amplitude square in the soft limit is given by \[\begin{align} \lim_{k_i\rightarrow 0}{{\mathbb{M}}^{(n,0)}(\dot{r})}=\lim_{k_i\rightarrow 0}{{\cal M}^{(n+1,0)}_{\left\{[1],S,0,J=S\right\}}(r)}&=&g_s^2\mathop{\sum_{k,l=n_{I}}}_{k,l\neq i, k\leq l}^{n_{L}^{(R)}+n_{H}}{\frac{k_k\cdot k_l}{k_k\cdot k_i k_l\cdot k_i}{\cal M}^{(n,0)}_{\left\{[1],S,0,J=S\right\},kl}(r^{i\backslash})}\nonumber\\ &=&g_s^2\mathop{\sum_{k,l=n_{I}}}_{k,l\neq i, k\leq l}^{n_{L}^{(R)}+n_{H}}{\frac{k_k\cdot k_l}{k_k\cdot k_i k_l\cdot k_i}{\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash})},\label{eq:softCSSwave} \end{align}\tag{25}\] where \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[1]}_J],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right). \end{align}\] In eq.@eq:eq:softCSSwave , we have used the light-cone axial gauge \[\begin{align} \sum_{\lambda_i}{\varepsilon^{\mu *}_{\lambda_i}(k_i)\varepsilon^{\nu}_{\lambda_i}(k_i)}&=&-g^{\mu\nu}+\frac{k_i^\mu n^\nu+k_i^\nu n^\mu}{k_i\cdot n},\label{eq:axialgaugesum} \end{align}\tag{26}\] where \(n^\mu\) is a light-like auxiliary vector (\(n^2=0\)). The terms involving \(n^{\mu}\) cancel out in eq.@eq:eq:softCSSwave as a result of colour conservation. Therefore, we have, for \(n_{I}\leq k \leq n_{L}^{(R)}+n_{H}+1\), \[\begin{align} {\mathsf M}^{(n-1,0)}_{{\rm soft}, kj_{Q}}(\dot{r}^{i\backslash})&=&0. \end{align}\]

3.2.2 Colour-octet S-wave state↩︎

For colour-octet S-wave states, we have derived the soft limit at the amplitude level in sect. 3.1 (cf. eq.@eq:eq:softOctetS ) as \[\begin{align} \lim_{k_i\rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[8],S,0,J=S\right\}}(r)}&=&\mathop{\sum_{j=n_{I}}}_{j\neq i}^{n_{L}^{(R)}+n_{H}}{g_s\frac{k_j\cdot \varepsilon_{\lambda_i}^*(k_i)}{k_j\cdot k_i}\vec{Q}({\cal I}_j){\cal A}^{(n,0)}_{\left\{[8],S,0,J=S\right\}}(r^{i\backslash})}\nonumber\\ &&+g_s\frac{K\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}\vec{Q}(Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[8]}_J]){\cal A}^{(n,0)}_{\left\{[8],S,0,J=S\right\}}(r^{i\backslash})\,. \end{align}\] The amplitude square in the soft limit is \[\begin{align} \lim_{k_i\rightarrow 0}{{\mathbb{M}}^{(n,0)}(\dot{r})}=\lim_{k_i\rightarrow 0}{{\cal M}^{(n+1,0)}_{\left\{[8],S,0,J=S\right\}}(r)}&=&g_s^2\!\mathop{\sum_{k,l=n_{I}}}_{k,l\neq i,k\leq l}^{n_{L}^{(R)}+n_{H}+1}\!\!{\frac{\tilde{k}_k\cdot \tilde{k}_l}{\tilde{k}_k\cdot k_i \tilde{k}_l\cdot k_i}{\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash})}\,,~~~\label{eq:softCOSwave} \end{align}\tag{27}\] where \(\tilde{k}_{n_{L}^{(R)}+n_{H}+1}=K\) and \(\tilde{k}_j=k_j\) when \(j\neq j_{Q}, j_{\bar{Q}^\prime}\) and \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[8]}_J],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right)\,. \end{align}\] Thanks to colour conservation, the \(n^\mu\) dependent terms stemming from eq.@eq:eq:axialgaugesum vanish on the right-hand side (r.h.s) of eq.@eq:eq:softCOSwave . Thus, we have, for \(n_{I}\leq k \leq n_{L}^{(R)}+n_{H}\), \[\begin{align} {\mathsf M}^{(n-1,0)}_{{\rm soft}, kj_{Q}}(\dot{r}^{i\backslash})&=&\frac{k_k\cdot K}{k_k\cdot k_i K\cdot k_i}{\mathbb{M}}^{(n-1,0)}_{kj_{Q}}(\dot{r}^{i\backslash})\,,\nonumber\\ {\mathsf M}^{(n-1,0)}_{{\rm soft}, j_{Q}j_{Q}}(\dot{r}^{i\backslash})&=&-\frac{K^2}{\left(K\cdot k_i\right)^2}C_{A}{\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash})\,. \end{align}\]

3.2.3 Colour-singlet spin-singlet P-wave state↩︎

With eq.@eq:eq:softSingletP2 in sect. 3.1, the soft limit of the amplitude for a colour-singlet spin-singlet P-wave quarkonium is \[\begin{align} \lim_{k_i\rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[1],0,1,1\right\}}(r)}&=&\mathop{\sum_{j=n_{I}}}_{j\neq i}^{n_{L}^{(R)}+n_{H}}{g_s\frac{k_j\cdot \varepsilon_{\lambda_i}^*(k_i)}{k_j\cdot k_i}\vec{Q}({\cal I}_j){\cal A}^{(n,0)}_{\left\{[1],0,1,1\right\}}(r^{i\backslash})}\nonumber\\ &&+g_s\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]\nonumber\\ &&\times\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[18]}){\cal A}^{(n,0)}_{\left\{[8],0,0,0\right\}}(r^{i\backslash}). \end{align}\] Then, the amplitude square in the soft limit becomes \[\begin{align} &&\lim_{k_i\rightarrow 0}{{\mathbb{M}}^{(n,0)}(\dot{r})}=\lim_{k_i\rightarrow 0}{{\cal M}^{(n+1,0)}_{\left\{[1],0,1,1\right\}}(r)}\nonumber\\ &=&g_s^2\mathop{\sum_{k,l=n_{I}}}_{k,l\neq i, k\leq l}^{n_{L}^{(R)}+n_{H}}{\frac{k_k\cdot k_l}{k_k\cdot k_i k_l\cdot k_i}{\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash})}\nonumber\\ &&+g_s^2\mathop{\sum_{k=n_{I}}}_{k\neq i}^{n_{L}^{(R)}+n_{H}}{\left[\frac{k_{k,\mu}}{k_k\cdot k_i K\cdot k_i}-\frac{K\cdot k_k k_{i,\mu}}{k_k\cdot k_i \left(K\cdot k_i\right)^2}\right]{\mathbb{M}}^{(n-1,0)}_{k[18]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu}\nonumber\\ &&-g_s^2\frac{2\epsilon-2}{\left(K\cdot k_i\right)^2}C_{{\rm eff}}(Q\bar{Q}^\prime_{[18]}){\mathbb{M}}^{(n-1,0)}(\dot{r}_1^{i\backslash}),\label{eq:softCSPwavesinglet} \end{align}\tag{28}\] where \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{1}P^{[1]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_1^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{1}S^{[8]}_0],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right), \end{align}\] and we have used the relations \[\begin{align} K\cdot \varepsilon^{(*)}_{\lambda_l}(K)&=&0,\tag{29}\\ \sum_{\lambda_l}{\varepsilon^{\mu,*}_{\lambda_l}(K)\varepsilon^{\nu}_{\lambda_l}(K)}&=&-g^{\mu\nu}+\frac{K^\mu K^\nu}{K^2}.\tag{30} \end{align}\] In numerical calculations, setting the dimensional regulator \(\epsilon\) to zero in the prefactor of the soft matrix element ensures consistency with the use of \(4\)-dimensional external wave functions. The terms dependent on \(n^\mu\) from eq.@eq:eq:axialgaugesum are absent on the r.h.s of eq.@eq:eq:softCSPwavesinglet due to colour conservation and \[\begin{align} \left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]_{\varepsilon_{\lambda_i}^*(k_i)\rightarrow k_i}&=&0. \end{align}\] The same observation holds for the other P-wave states discussed later, and there is no need to reiterate this point.

Thus, we have, for \(n_{I}\leq k \leq n_{L}^{(R)}+n_{H}\), \[\begin{align} {\mathsf M}^{(n-1,0)}_{{\rm soft}, kj_{Q}}(\dot{r}^{i\backslash})&=&\left[\frac{k_{k,\mu}}{k_k\cdot k_i K\cdot k_i}-\frac{K\cdot k_k k_{i,\mu}}{k_k\cdot k_i \left(K\cdot k_i\right)^2}\right]{\mathbb{M}}^{(n-1,0)}_{k[18]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu,\nonumber\\ {\mathsf M}^{(n-1,0)}_{{\rm soft}, j_{Q}j_{Q}}(\dot{r}^{i\backslash})&=&-\frac{2\epsilon-2}{\left(K\cdot k_i\right)^2}C_{{\rm eff}}(Q\bar{Q}^\prime_{[18]}){\mathbb{M}}^{(n-1,0)}(\dot{r}_1^{i\backslash}). \end{align}\]

3.2.4 Colour-octet spin-singlet P-wave state↩︎

For the production of a colour-octet spin-singlet P-wave quarkonium, we have established the soft limit at the amplitude level (cf. eq.@eq:eq:softOctetP2 in sect. 3.1) as follows: \[\begin{align} \lim_{k_i\rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[8],0,1,1\right\}}(r)}&=&\mathop{\sum_{j=n_{I}}}_{j\neq i}^{n_{L}^{(R)}+n_{H}}{g_s\frac{k_j\cdot \varepsilon_{\lambda_i}^*(k_i)}{k_j\cdot k_i}\vec{Q}({\cal I}_j){\cal A}^{(n,0)}_{\left\{[8],0,1,1\right\}}(r^{i\backslash})}\nonumber\\ &&+g_s\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}\vec{Q}(Q\bar{Q}^\prime[\bigl.^{1}P^{[8]}_1]){\cal A}^{(n,0)}_{\left\{[8],0,1,1\right\}}(r^{i\backslash})\nonumber\\ &&+g_s\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]\nonumber\\ &&\times\left[\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]}){\cal A}^{(n,0)}_{\left\{[1],0,0,0\right\}}(r^{i\backslash})+\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]}){\cal A}^{(n,0)}_{\left\{[8],0,0,0\right\}}(r^{i\backslash})\right].\nonumber\\ \end{align}\] The amplitude square in the soft limit is \[\begin{align} &&\lim_{k_i\rightarrow 0}{{\mathbb{M}}^{(n,0)}(\dot{r})}=\lim_{k_i\rightarrow 0}{{\cal M}^{(n+1,0)}_{\left\{[8],0,1,1\right\}}(r)}\nonumber\\ &=&g_s^2\mathop{\sum_{k,l=n_{I}}}_{k,l\neq i, k\leq l}^{n_{L}^{(R)}+n_{H}+1}{\frac{\tilde{k}_k\cdot \tilde{k}_l}{\tilde{k}_k\cdot k_i \tilde{k}_l\cdot k_i}{\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash})}\nonumber\\ &&+g_s^2\mathop{\sum_{k=n_{I}}}_{k\neq i}^{n_{L}^{(R)}+n_{H}+1}{\left[\frac{\tilde{k}_{k,\mu}}{\tilde{k}_k\cdot k_i K\cdot k_i}-\frac{K\cdot \tilde{k}_k k_{i,\mu}}{\tilde{k}_k\cdot k_i \left(K\cdot k_i\right)^2}\right]\left[{\mathbb{M}}^{(n-1,0)}_{k[88]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu+ {\mathbb{M}}^{(n-1,0)}_{k[81]}(\dot{r}^{i\backslash},\dot{r}_2^{i\backslash})^\mu\right]}\nonumber\\ &&-g_s^2\frac{2\epsilon-2}{\left(K\cdot k_i\right)^2}\left[C_{{\rm eff}} (Q\bar{Q}^\prime_{[88]}){\mathbb{M}}^{(n-1,0)}(\dot{r}_1^{i\backslash})+C_{{\rm eff}} (Q\bar{Q}^\prime_{[81]}){\mathbb{M}}^{(n-1,0)}(\dot{r}_2^{i\backslash})\right], \end{align}\] where \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{1}P^{[8]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_1^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{1}S^{[8]}_0],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_2^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{1}S^{[1]}_0],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right), \end{align}\] and \(\tilde{k}_{n_{L}^{(R)}+n_{H}+1}=K\) and \(\tilde{k}_j=k_j\) when \(j\neq j_{Q}, j_{\bar{Q}^\prime}\). Additionally, we have used the relations eq.@eq:eq:polrelation1 .

Thus, we have, for \(n_{I}\leq k \leq n_{L}^{(R)}+n_{H}\), \[\begin{align} {\mathsf M}^{(n-1,0)}_{{\rm soft}, kj_{Q}}(\dot{r}^{i\backslash})&=&\frac{k_k\cdot K}{k_k\cdot k_i K\cdot k_i}{\mathbb{M}}^{(n-1,0)}_{kj_{Q}}(\dot{r}^{i\backslash})\nonumber\\ &&+\left[\frac{k_{k,\mu}}{k_k\cdot k_i K\cdot k_i}-\frac{K\cdot k_k k_{i,\mu}}{k_k\cdot k_i \left(K\cdot k_i\right)^2}\right]\left[{\mathbb{M}}^{(n-1,0)}_{k[88]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu+ {\mathbb{M}}^{(n-1,0)}_{k[81]}(\dot{r}^{i\backslash},\dot{r}_2^{i\backslash})^\mu\right],\nonumber\\ {\mathsf M}^{(n-1,0)}_{{\rm soft}, j_{Q}j_{Q}}(\dot{r}^{i\backslash})&=&-\frac{K^2}{\left(K\cdot k_i\right)^2}C_{A}{\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash})\nonumber\\ &&+\left[\frac{K_{\mu}}{\left(K\cdot k_i\right)^2}-\frac{K^2 k_{i,\mu}}{\left(K\cdot k_i\right)^3}\right]\left[{\mathbb{M}}^{(n-1,0)}_{j_{Q}[88]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu+{\mathbb{M}}^{(n-1,0)}_{j_{Q}[81]}(\dot{r}^{i\backslash},\dot{r}_2^{i\backslash})^\mu\right]\nonumber\\ &&-\frac{2\epsilon-2}{\left(K\cdot k_i\right)^2}\left[C_{{\rm eff}} (Q\bar{Q}^\prime_{[88]}){\mathbb{M}}^{(n-1,0)}(\dot{r}_1^{i\backslash})+C_{{\rm eff}} (Q\bar{Q}^\prime_{[81]}){\mathbb{M}}^{(n-1,0)}(\dot{r}_2^{i\backslash})\right]. \end{align}\]

3.2.5 Colour-singlet spin-triplet P-wave state↩︎

We have derived the soft limit at the amplitude level for a colour-singlet spin-triplet P-wave state in sect. 3.1 (cf. eq.@eq:eq:softSingletP3 ) as \[\begin{align} \lim_{k_i\rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[1],1,1,J\right\}}(r)}&=&\mathop{\sum_{j=n_{I}}}_{j\neq i}^{n_{L}^{(R)}+n_{H}}{g_s\frac{k_j\cdot \varepsilon_{\lambda_i}^*(k_i)}{k_j\cdot k_i}\vec{Q}({\cal I}_j){\cal A}^{(n,0)}_{\left\{[1],1,1,J\right\}}(r^{i\backslash})}\nonumber\\ &&+g_s\sum_{\lambda_l,\lambda_s}{\langle J,\lambda_j|1,\lambda_l;1,\lambda_s\rangle \left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]}\nonumber\\ &&\times\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[18]}){\cal A}^{(n,0)}_{\left\{[8],1,0,1\right\}}(r^{i\backslash})\,. \end{align}\] The amplitude square in the soft limit is \[\begin{align} &&\lim_{k_i\rightarrow 0}{{\mathbb{M}}^{(n,0)}(\dot{r})}=\lim_{k_i\rightarrow 0}{{\cal M}^{(n+1,0)}_{\left\{[1],1,1,J\right\}}(r)}\nonumber\\ &=&g_s^2\mathop{\sum_{k,l=n_{I}}}_{k,l\neq i, k\leq l}^{n_{L}^{(R)}+n_{H}}{\frac{k_k\cdot k_l}{k_k\cdot k_i k_l\cdot k_i}{\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash})}\nonumber\\ &&+g_s^2\mathop{\sum_{k=n_{I}}}_{k\neq i}^{n_{L}^{(R)}+n_{H}}{\left[\frac{k_{k,\mu}}{k_k\cdot k_i K\cdot k_i}-\frac{K\cdot k_k k_{i,\mu}}{k_k\cdot k_i \left(K\cdot k_i\right)^2}\right]{\mathbb{M}}^{(n-1,0)}_{J,k[18]}(\dot{r},\dot{r}_1)^{\mu}}\nonumber\\ &&-g_s^2\left[\frac{g_{\mu\nu}}{\left(K\cdot k_i\right)^2}+\frac{K^2 k_{i,\mu}k_{i,\nu}}{\left(K\cdot k_i\right)^4}\right]C_{{\rm eff}}(Q\bar{Q}^\prime_{[18]}){\mathbb{M}}^{(n-1,0)}_J(\dot{r}_1^{i\backslash})^{\mu\nu}\,, \end{align}\] where \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{3}P^{[1]}_J],\bar{Q}^\prime\backslash\,, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_1^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{3}S^{[8]}_1],\bar{Q}^\prime\backslash\,, \ldots,{\cal I}_{n+3}\right), \end{align}\] and we have used the relations \[\begin{align} K\cdot \varepsilon^{(*)}_{\lambda_l}(K)&=&K\cdot \varepsilon^{(*)}_{\lambda_s}(K)=0\,,\nonumber\\ \sum_{\lambda_l}{\varepsilon^{\mu,*}_{\lambda_l}(K)\varepsilon^{\nu}_{\lambda_l}(K)}&=&-g^{\mu\nu}+\frac{K^\mu K^\nu}{K^2}=\Pi^{\mu\nu}(K)\,,\nonumber\\ \varepsilon_{0,0}^{\mu\nu,*}(K)\varepsilon_{0,0}^{\alpha\beta}(K)&=&\frac{1}{3}\Pi^{\mu\nu}(K)\Pi^{\alpha\beta}(K)\,,\nonumber\\ \sum_{\lambda_j=-1}^{1}{\varepsilon_{1,\lambda_j}^{\mu\nu,*}(K)\varepsilon_{1,\lambda_j}^{\alpha\beta}(K)}&=&\frac{1}{2}\left[\Pi^{\mu\alpha}(K)\Pi^{\nu\beta}(K)-\Pi^{\mu\beta}(K)\Pi^{\nu\alpha}(K)\right]\,,\nonumber\\ \sum_{\lambda_j=-2}^{2}{\varepsilon_{2,\lambda_j}^{\mu\nu,*}(K)\varepsilon_{2,\lambda_j}^{\alpha\beta}(K)}&=&\frac{1}{2}\left[\Pi^{\mu\alpha}(K)\Pi^{\nu\beta}(K)+\Pi^{\mu\beta}(K)\Pi^{\nu\alpha}(K)\right]-\frac{1}{3}\Pi^{\mu\nu}(K)\Pi^{\alpha\beta}(K)\,.\nonumber\\ \end{align}\]

Then, we can easily derive, \[\begin{align} {\mathsf M}^{(n-1,0)}_{{\rm soft}, kj_{Q}}(\dot{r}^{i\backslash})&=&\left[\frac{k_{k,\mu}}{k_k\cdot k_i K\cdot k_i}-\frac{K\cdot k_k k_{i,\mu}}{k_k\cdot k_i \left(K\cdot k_i\right)^2}\right]{\mathbb{M}}^{(n-1,0)}_{J,k[18]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu\,,\nonumber\\ {\mathsf M}^{(n-1,0)}_{{\rm soft}, j_{Q}j_{Q}}(\dot{r}^{i\backslash})&=&-\left[\frac{g_{\mu\nu}}{\left(K\cdot k_i\right)^2}+\frac{K^2 k_{i,\mu}k_{i,\nu}}{\left(K\cdot k_i\right)^4}\right]C_{{\rm eff}}(Q\bar{Q}^\prime_{[18]}){\mathbb{M}}^{(n-1,0)}_J(\dot{r}_1^{i\backslash})^{\mu\nu}\,, \end{align}\] with \(n_{I}\leq k \leq n_{L}^{(R)}+n_{H}\).

3.2.6 Colour-octet spin-triplet P-wave state↩︎

Finally, for the colour-octet spin-triplet P-wave state production, the soft limit of the amplitude from sect. 3.1 (cf. eq.@eq:eq:softOctetP3 ) is \[\begin{align} \lim_{k_i\rightarrow 0}{{\cal A}^{(n+1,0)}_{\left\{[8],1,1,J\right\}}(r)}&=&\mathop{\sum_{j=n_{I}}}_{j\neq i}^{n_{L}^{(R)}+n_{H}}{g_s\frac{k_j\cdot \varepsilon_{\lambda_i}^*(k_i)}{k_j\cdot k_i}\vec{Q}({\cal I}_j){\cal A}^{(n,0)}_{\left\{[8],1,1,J\right\}}(r^{i\backslash})}\nonumber\\ &&+g_s\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}\vec{Q}(Q\bar{Q}^\prime[\bigl.^{3}P^{[8]}_J]){\cal A}^{(n,0)}_{\left\{[8],1,1,J\right\}}(r^{i\backslash})\nonumber\\ &&+g_s\sum_{\lambda_l,\lambda_s}{\langle J,\lambda_j|1,\lambda_l; 1, \lambda_s\rangle\left[\frac{\varepsilon_{\lambda_l}^*(K)\cdot \varepsilon_{\lambda_i}^*(k_i)}{K\cdot k_i}-\frac{K \cdot \varepsilon_{\lambda_i}^*(k_i) k_i\cdot \varepsilon_{\lambda_l}^*(K)}{\left(K\cdot k_i\right)^2}\right]}\nonumber\\ &&\times\left[\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[81]}){\cal A}^{(n,0)}_{\left\{[1],1,0,1\right\}}(r^{i\backslash})+\vec{Q}_{{\rm eff}}(Q\bar{Q}^\prime_{[88]}){\cal A}^{(n,0)}_{\left\{[8],1,0,1\right\}}(r^{i\backslash})\right]. \end{align}\] The amplitude square in the soft limit is \[\begin{align} &&\lim_{k_i\rightarrow 0}{{\mathbb{M}}^{(n,0)}(\dot{r})}=\lim_{k_i\rightarrow 0}{{\cal M}^{(n+1,0)}_{\left\{[8],1,1,J\right\}}(r)}\nonumber\\ &=&g_s^2\mathop{\sum_{k,l=n_{I}}}_{k,l\neq i, k\leq l}^{n_{L}^{(R)}+n_{H}+1}{\frac{\tilde{k}_k\cdot \tilde{k}_l}{\tilde{k}_k\cdot k_i \tilde{k}_l\cdot k_i}{\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash})}\nonumber\\ &&+g_s^2\mathop{\sum_{k=n_{I}}}_{k\neq i}^{n_{L}^{(R)}+n_{H}+1}{\left[\frac{\tilde{k}_{k,\mu}}{\tilde{k}_k\cdot k_i K\cdot k_i}-\frac{K\cdot \tilde{k}_k k_{i,\mu}}{\tilde{k}_k\cdot k_i \left(K\cdot k_i\right)^2}\right]\left[{\mathbb{M}}^{(n-1,0)}_{J,k[88]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu+ {\mathbb{M}}^{(n-1,0)}_{J,k[81]}(\dot{r}^{i\backslash},\dot{r}_2^{i\backslash})^\mu\right]}\nonumber\\ &&-g_s^2\left[\frac{g_{\mu\nu}}{\left(K\cdot k_i\right)^2}+\frac{K^2 k_{i,\mu}k_{i,\nu}}{\left(K\cdot k_i\right)^4}\right]\left[C_{{\rm eff}} (Q\bar{Q}^\prime_{[88]}){\mathbb{M}}^{(n-1,0)}_{J}(\dot{r}_1^{i\backslash})^{\mu\nu}+C_{{\rm eff}} (Q\bar{Q}^\prime_{[81]}){\mathbb{M}}^{(n-1,0)}_{J}(\dot{r}_2^{i\backslash})^{\mu\nu}\right],\nonumber\\ \end{align}\] where \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{3}P^{[8]}_J],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_1^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{3}S^{[8]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_2^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{3}S^{[1]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right), \end{align}\] and \(\tilde{k}_{n_{L}^{(R)}+n_{H}+1}=K\) and \(\tilde{k}_j=k_j\) when \(j\neq j_{Q}, j_{\bar{Q}^\prime}\). Once again, we have used the relations eq.@eq:eq:polrelation1 .

Thus, we can easily get, for \(n_{I}\leq k \leq n_{L}^{(R)}+n_{H}\), \[\begin{align} {\mathsf M}^{(n-1,0)}_{{\rm soft}, kj_{Q}}(\dot{r}^{i\backslash})&=&\frac{k_k\cdot K}{k_k\cdot k_i K\cdot k_i}{\mathbb{M}}^{(n-1,0)}_{kj_{Q}}(\dot{r}^{i\backslash})\nonumber\\ &&+\left[\frac{k_{k,\mu}}{k_k\cdot k_i K\cdot k_i}-\frac{K\cdot k_k k_{i,\mu}}{k_k\cdot k_i \left(K\cdot k_i\right)^2}\right]\left[{\mathbb{M}}^{(n-1,0)}_{J,k[88]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu+ {\mathbb{M}}^{(n-1,0)}_{J,k[81]}(\dot{r}^{i\backslash},\dot{r}_2^{i\backslash})^\mu\right],\nonumber\\ {\mathsf M}^{(n-1,0)}_{{\rm soft}, j_{Q}j_{Q}}(\dot{r}^{i\backslash})&=&-\frac{K^2}{\left(K\cdot k_i\right)^2}C_{A}{\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash})\nonumber\\ &&+\left[\frac{K_{\mu}}{\left(K\cdot k_i\right)^2}-\frac{K^2 k_{i,\mu}}{\left(K\cdot k_i\right)^3}\right]\left[{\mathbb{M}}^{(n-1,0)}_{J,j_{Q}[88]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu+{\mathbb{M}}^{(n-1,0)}_{J,j_{Q}[81]}(\dot{r}^{i\backslash},\dot{r}_2^{i\backslash})^\mu\right]\nonumber\\ &&-\left[\frac{g_{\mu\nu}}{\left(K\cdot k_i\right)^2}+\frac{K^2 k_{i,\mu}k_{i,\nu}}{\left(K\cdot k_i\right)^4}\right]\left[C_{{\rm eff}} (Q\bar{Q}^\prime_{[88]}){\mathbb{M}}^{(n-1,0)}_{J}(\dot{r}_1^{i\backslash})^{\mu\nu}\right.\nonumber\\ &&\left.+C_{{\rm eff}} (Q\bar{Q}^\prime_{[81]}){\mathbb{M}}^{(n-1,0)}_{J}(\dot{r}_2^{i\backslash})^{\mu\nu}\right]\,. \end{align}\]

4 FKS subtraction for single quarkonium production↩︎

4.1 FKS pairs and partition functions↩︎

With the squared real amplitudes in the soft limit at hand, we will now structure our approach using the FKS formalism. To categorise the IR divergences and facilitate their subtraction, we introduce a set of ordered pairs for any given process \(\dot{r}\in \dot{{\cal R}}_{n}\). This set is denoted as the set of FKS pairs \[\begin{align} {\cal P}_{\rm FKS}(\dot{r})&=&\Big\{(i,j)\;\Big|\;3\le i\le n_{L}^{(R)}+2\,, n_{I}\le j\le n_{L}^{(R)}+n_{H}+1\,, i\ne j\,, \nonumber\\*&&\phantom{aaa} {\mathbb{M}}^{(n,0)}(\dot{r})J^{n_{L}^{(B)}}\rightarrow\infty~~{\rm if}~~\tilde{k}_i^0\rightarrow 0~~ {\rm or}~~\tilde{k}_j^0\rightarrow 0~~{\rm or}~~\vec{\tilde{k}}_i\parallel \vec{\tilde{k}}_j\Big\}. \phantom{aaaaa} \label{PFKSdef} \end{align}\tag{31}\] This implies that a pair of particles belongs to the set of FKS pairs if they induce soft or collinear singularities (or both) in the \(n\)-body real matrix elements. It is important to note that \(\tilde{k}_j^0\rightarrow 0\) is irrelevant when \(j=1,2\). In the calculation of an NLO cross section within the FKS formalism, each pair belonging to \({\cal P}_{\rm FKS}\) corresponds to a set of subtractions of soft and collinear singularities.

In the FKS formalism, we multiply the real emission matrix elements with the so-called measurement/partition function \({\cal S}_{ij}\), so that the phase space is partitioned into different kinematic regions where each region contains at most one soft and one collinear singularity. The partitioning is accomplished through the introduction of a set of positive-definite functions \[\begin{align} {\cal S}_{ij}(\dot{r})\,,\;\;\;\;\;\;(i,j)\in{\cal P}_{\rm FKS}(\dot{r})\,, \end{align}\] where the argument \(\dot{r}\in \dot{{\cal R}}_{n}\) means that we can choose different \({\cal S}\) for different processes. As described in ref. [56], \({\cal S}_{ij}\) is defined in different regions as following: \[\begin{align} \!\!\!\sum_{(i,j)\in {\cal P}_{\rm FKS}(\dot{r})}{{\cal S}_{ij}(\dot{r})} &=&1, \nonumber \\ \lim_{\vec{\tilde{k}}_i\parallel\vec{\tilde{k}}_j}{{\cal S}_{ij}(\dot{r})} &=& h_{ij}\left(\frac{\tilde{k}_i^0}{\tilde{k}_i^0+\tilde{k}_j^0}\right),\,~~~~{\rm if}~~\tilde{m}_i=\sqrt{\tilde{k}_i^2}=\tilde{m}_j=\sqrt{\tilde{k}_j^2}=0\,,\nonumber \\ \lim_{\tilde{k}_i^0 \rightarrow 0}{{\cal S}_{ij}(\dot{r})} &=& c_{ij}\,,\phantom{}\; {\rm if}~~~{\cal I}_i=g\,,\,~~~~{\rm with}~~0<c_{ij}\le 1\qquad{\rm and}\,\,\,\,\, \mathop{\sum_{j}}_{(i,j)\in{\cal P}_{\rm FKS}(\dot{r})}\!\!\!\!\!c_{ij}=1\,,\nonumber\\ \lim_{\vec{\tilde{k}}_k\parallel\vec{\tilde{k}}_l}{{\cal S}_{ij}(\dot{r})}&=&0~~\forall\, \{k,l\}\ne\{i,j\}~~{\rm with}~~(k,l)\in{\cal P}_{\rm FKS}(\dot{r})~~{\land}~~ \tilde{m}_k=\tilde{m}_l=0\,,\nonumber\\ \lim_{\tilde{k}_k^0\rightarrow 0}{{\cal S}_{ij}(\dot{r})}&=&0~~\forall k~~{\rm with}~~{\cal I}_k=g ~~{\rm and}~~\exists l~~{\rm with}~~(k,l)\in{\cal P}_{\rm FKS}(\dot{r})~{\lor}~ (l,k)\in{\cal P}_{\rm FKS}(\dot{r}). \nonumber\\&& \label{eq:SfunC} \end{align}\tag{32}\] In other words, \({\cal S}_{ij}\) goes to zero in all regions of the phase space where the real emission matrix elements diverge, except if this involves particle \(i\) being soft, or particles \(i\) and \(j\) being collinear. The functions \(h_{ij}(z)\) introduced in eq.@eq:eq:SfunC are defined in \(0\le z\le 1\), and have the following properties: \[\begin{align} h_{ij}(z)&=&1\,,~~~~~~~~~{\rm if}~~n_{I}\le j\le 2\,, \tag{33} \\ h_{ij}(z)&=&h(z)\,,~~~~~{\rm if}~~3\le j\le n_{L}^{(R)}+n_{H}+1\,, \tag{34} \end{align}\] with \(h(z)\) a positive-definite function such that \[\lim_{z\rightarrow 0}{h(z)}=1\,,\;\;\;\;\;\; \lim_{z\rightarrow 1}{h(z)}=0\,,\;\;\;\;\;\; h(z)+h(1-z)=1\,. \label{eq:hdef3}\tag{35}\] Note that in \(\dot{r}\), we have glued \(Q\) and \(\bar{Q}^\prime\) to a single particle \(Q\bar{Q}^\prime[n]\). This means that the number of strongly interacting massive particles (except \(n\) is a colour-singlet S-wave state) is \(n_{H}-1\) instead of \(n_{H}\).

The real matrix elements can be rewritten \[\begin{align} {\mathbb{M}}^{(n,0)}(\dot{r})&=&\sum_{(i,j)\in {\cal P}_{\rm FKS}(\dot{r})}{{\cal S}_{ij}(\dot{r}){\mathbb{M}}^{(n,0)}(\dot{r})}, \end{align}\] which will only be singular for a given term on the r.h.s if particle \(i\) is soft and/or particles \(i\) and \(j\) are collinear after applying \(J^{n_{L}^{(B)}}\).

4.2 Local subtraction counterterms↩︎

The IR divergence subtracted real cross sections can be formulated in the center of mass frame of the incoming partons: \[\begin{align} k_1=\tilde{k}_1=\frac{\sqrt{s}}{2}(1,0,0,1)\,,\;\;\;\;\; k_2=\tilde{k}_2=\frac{\sqrt{s}}{2}(1,0,0,-1)\,. \end{align}\] In this frame, for each pair \((i,j)\in {\cal P}_{\rm FKS}(\dot{r})\), we can introduce the variables \(\xi_i\) and \(y_{ij}\), where \[\begin{align} \tilde{k}_i^0&=&\frac{\sqrt{s}}{2}\xi_i\,, \tag{36} \\ \vec{\tilde{k}}_i\!\cdot\!\vec{\tilde{k}}_j&=&\left|\vec{\tilde{k}}_i\right|\left|\vec{\tilde{k}}_j\right|y_{ij}\,. \tag{37} \end{align}\] Thus, \(\xi_i\) is the rescaled energy of the FKS parton \(i\), and \(y_{ij}\) is the cosine of the angle between the FKS parton \(i\) and its sister \(j\). The soft and collinear singularities of \({\cal S}_{ij}(\dot{r}){\mathbb{M}}^{(n,0)}(\dot{r})J^{n_{L}^{(B)}}\) correspond to \(\xi_i=0\) and to \(y_{ij}=1\), respectively. The IR-divergence locally subtracted partonic cross section is \[\begin{align} d\hat{\sigma}_{{\rm FKS}}(\dot{r})&=&\sum_{(i,j)\in{\cal P}_{\rm FKS}(\dot{r})}{d\hat{\sigma}_{ij}(\dot{r})}, \end{align}\] where the IR-divergence locally subtracted real partonic cross section is \[d\hat{\sigma}_{ij}(\dot{r})=\left(\frac{1}{\xi_i}\right)_c\left(\frac{1}{1-y_{ij}}\right)_\delta \Big((1-y_{ij})\xi_i^2{\mathbb{M}}^{(n,0)}(\dot{r})\Big) {\cal S}_{ij}(\dot{r})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r})}{\cal G}(\dot{r})\, d\xi_idy_{ij}d\varphi_id\widetilde{\phi}_{n-1}^{ij}(\dot{r})\,. \label{eq:dsigijnpo}\tag{38}\] The variable \(\varphi_i\) is the azimuthal direction of the FKS parton. The quantity \(d\widetilde{\phi}_{n-1}^{ij}\) is the reduced \((n-1)\)-body phase space via the following relation: \[\begin{align} d\phi_{n}(\dot{r})&=&\xi_i^{1-2\epsilon} d\xi_i(1-y_{ij}^2)^{-\epsilon}dy_{ij}d\Omega_{i}^{(2-2\epsilon)} d\widetilde{\phi}_{n-1}^{ij}(\dot{r})\,. \label{eq:tphspdef} \end{align}\tag{39}\] The reduced phase space measure has the following limits: \[\begin{align} \lim_{\xi_i\rightarrow 0}{d\widetilde{\phi}_{n-1}^{ij}(\dot{r})}&=&\frac{s^{1-\epsilon}}{(4\pi)^{3-2\epsilon}}d\phi_{n-1}(\dot{r}^{i\backslash})\,,\,~~~~~~~~~{\rm if}~~\tilde{m}_i=0\,, \tag{40} \\ \lim_{y_{ij}\rightarrow 1}{d\widetilde{\phi}_{n-1}^{ij}(\dot{r})}&=&\frac{s^{1-\epsilon}}{(4\pi)^{3-2\epsilon}}d\phi_{n-1}(\dot{r}^{j\oplus i,i\backslash})\,,~~~~~{\rm if}~~\tilde{m}_i=\tilde{m}_j=0\,. \tag{41} \end{align}\] In \(d=4\) dimensions, we can simply set \(\epsilon=0\) and \(d\Omega_i^{(2-2\epsilon)}=d\varphi_i\). The distributions entering eq.@eq:eq:dsigijnpo are defined as follows, for any test functions \(f()\) and \(g()\): \[\begin{align} \int_0^{\xi_{\rm max}}{d\xi_if(\xi_i)\left(\frac{1}{\xi_i}\right)_c}&=& \int_0^{\xi_{\rm max}}{d\xi_i\frac{f(\xi_i)-f(0)\Theta(\xi_{cut}-\xi_i)}{\xi_i}}\,, \tag{42} \\ \int_{-1}^1{dy_{ij}g(y_{ij})\left(\frac{1}{1-y_{ij}}\right)_\delta}&=& \int_{-1}^1{dy_{ij}\frac{g(y_{ij})-g(1)\Theta(y_{ij}-1+\delta)}{1-y_{ij}}}\,, \tag{43} \end{align}\] where \[\xi_{\rm max}=1-\frac{1}{s}\left(\sum_{k=3}^{n+2}\tilde{m}_k\right)^2\,,\] and \(\Theta()\) is the Heaviside theta function. Note that \(\tilde{k}_j=k_j, \tilde{m}_{j}=\sqrt{\tilde{k}_j^2}=\sqrt{k_j^2}=m_j\) when \(j\leq n_{L}^{(R)}+n_{H}\), \(\tilde{k}_{n_{L}^{(R)}+n_{H}+1}=k_{n_{L}^{(R)}+n_{H}+1}+k_{n_{L}^{(R)}+n_{H}+2}=K,\tilde{m}_{n_{L}^{(R)}+n_{H}+1}=m_{Q}+m_{\bar{Q}^\prime}\), and \(\tilde{k}_j=k_{j+1},\tilde{m}_j=m_{j+1}\) when \(j\geq n_{L}^{(R)}+n_{H}+3\). In eqs.@eq:eq:distrxii and 43 , \(\xi_{cut}\) and \(\delta\) are free parameters, that can be chosen in the ranges \[0<\xi_{cut}\le\xi_{\rm max}\,,\;\;\;\;\; 0<\delta\le 2\,.\] In MadFKS [56], \(\delta=\delta_{I}\) for the initial state collinear singularities \((i,j)\in {\cal P}_{\rm FKS}(\dot{r}), j\leq 2\), and \(\delta=\delta_{O}\) for the final state collinear singularities \((i,j)\in {\cal P}_{\rm FKS}(\dot{r}) , j\geq 3\). Now, we can introduce the quantity in \(d=4\) dimensions \[\begin{align} \Sigma_{ij}(\dot{r}; \xi_i, y_{ij})&=&\Big((1-y_{ij})\xi_i^2{\mathbb{M}}^{(n,0)}(\dot{r})\Big) {\cal S}_{ij}(\dot{r})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r})}{\cal G}(\dot{r})\, d\varphi_id\widetilde{\phi}_{n-1}^{ij}(\dot{r}), \end{align}\] so that \[\begin{align} d\hat{\sigma}_{ij}(\dot{r})=\left(\frac{1}{\xi_i}\right)_c\left(\frac{1}{1-y_{ij}}\right)_\delta\Sigma_{ij}(\dot{r}; \xi_i, y_{ij}) d\xi_idy_{ij}. \end{align}\] If we expand the plus distributions, we have \[\begin{align} d\hat{\sigma}_{ij}(\dot{r})&=& \frac{1}{\xi_i(1-y_{ij})} \Big[\Sigma_{ij}(\dot{r}; \xi_i,y_{ij})-\Sigma_{ij}(\dot{r};\xi_i,1)\Theta(y_{ij}-1+\delta) \nonumber\\*&&-\Sigma_{ij}(\dot{r}; 0,y_{ij})\Theta(\xi_{cut}-\xi_i) +\Sigma_{ij}(\dot{r}; 0,1)\Theta(\xi_{cut}-\xi_i)\Theta(y_{ij}-1+\delta) \Big]d\xi_idy_{ij}. \nonumber\\*&& \label{eq:dsigijnpoE} \end{align}\tag{44}\] The first term in the integrand, called “events", contributes to the initial partonic real cross section, while the remaining three terms are the local subtraction counterterms, which are called”collinear counterevent", “soft counterevent", and”soft-collinear counterevent", respectively. The collinear and soft-collinear local counterterms can be inferred from the collinear limit of the real-emission matrix elements \[\begin{align} \lim_{\vec{\tilde{k}}_i\parallel \vec{\tilde{k}}_j}{(1-y_{ij})\xi_i^2{\mathbb{M}}^{(n,0)}(\dot{r})}&=&\frac{4}{s}g_s^2\mu^{2\epsilon}\xi_iP_{{\cal I}_{j\oplus \bar{i}}{\cal I}_j}^{<}(1-\xi_i,\epsilon){\mathbb{M}}^{(n-1,0)}(\dot{r}^{j\oplus \bar{i},i\backslash})\nonumber\\ &&+\frac{4}{s}g_s^2\mu^{2\epsilon}\xi_i\underbrace{Q_{{\cal I}_{j\oplus \bar{i}}^\star{\cal I}_j}(1-\xi_i)\tilde{{\mathbb{M}}}^{(n-1,0)}_{ij}(\dot{r}^{j\oplus \bar{i},i\backslash})}_{\equiv \Delta_{ij}},\tag{45}\\ &&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~j=n_{I},\ldots,2,\nonumber\\ \lim_{\vec{\tilde{k}}_i\parallel \vec{\tilde{k}}_j}{(1-y_{ij})\xi_i^2{\mathbb{M}}^{(n,0)}(\dot{r}){\cal S}_{ij}(\dot{r})}&=&\frac{4}{s}g_s^2\mu^{2\epsilon}\frac{1-z}{z}h(z)P_{{\cal I}_j{\cal I}_{j\oplus i}}^{<}(z,\epsilon){\mathbb{M}}^{(n-1,0)}(\dot{r}^{j\oplus i,i\backslash})\nonumber\\ &&+\frac{4}{s}g_s^2\mu^{2\epsilon}\frac{1-z}{z}h(z)\underbrace{Q_{{\cal I}_j{\cal I}_{j\oplus i}^\star}(z)\tilde{{\mathbb{M}}}^{(n-1,0)}_{ij}(\dot{r}^{j\oplus i,i\backslash})}_{\equiv \Delta_{ij}},\tag{46}\\ &&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~j=3,\ldots, n_{L}^{(R)}+2, j\neq i,\nonumber \end{align}\] where \(\bar{i}\) denotes the antiparticle of the particle \(i\), \(P_{{\cal I}_k{\cal I}_l}^{<}(z,\epsilon)\) is the unregularised Altarelli-Parisi kernel for \(z<1\) in \(d=4-2\epsilon\) dimensions that can be found in the literature (see, e.g., eqs.(D.15-D.18) in ref. [56]), and \(Q_{{\cal I}_{j\oplus \bar{i}}^\star{\cal I}_j}\) and \(Q_{{\cal I}_j{\cal I}_{j\oplus i}^\star}\) are given in eqs.(D.3-D.10) of ref. [56]. The reduced matrix element is defined as [50], [56] \[\begin{align} \tilde{{\mathbb{M}}}^{(n-1,0)}_{ij}(\dot{r}^{j\oplus i,i\backslash})&=&\frac{1}{2s} \frac{1}{\omega({\cal I}_1)\omega({\cal I}_2)} \Re{\left\{\frac{\langle ij\rangle}{[ij]}\tilde{\mathop{\sum_{\rm colour}}_{\rm spin}}\! {\mathbb{A}}^{(n-1,0)}_+(\dot{r}^{j\oplus i,i\backslash}) {{\mathbb{A}}^{(n-1,0)}_-(\dot{r}^{j\oplus i,i\backslash})}^{\star}\right\}}\,,~~~~~ \end{align}\] where the spinor-helicity formalism takes the conventions of ref. [72], \({\mathbb{A}}^{(n-1,0)}_{\pm}\) represents the helicity amplitude with the helicity of the parton \(j\oplus i\) being \(\pm\), and the sum with a tilde has summed over the colour and spin of the external states except the spin of the particle \(j\oplus i\). The \(\Delta_{ij}\) term vanishes upon the integration of the azimuthal variable \(d\varphi_i\) of the parton \(i\). In eq.@eq:eq:localfinalcoll0 , we have assumed \(\tilde{k}_i=(1-z)\left(\tilde{k}_i+\tilde{k}_j\right)\) and \(\tilde{k}_j=z\left(\tilde{k}_i+\tilde{k}_j\right)\). This essentially gives us for the initial-state collinear singularities (\(n_{I}\leq j\leq 2\)) \[\begin{align} \Sigma_{ij}(\dot{r}; \xi_i, 1)&=&\frac{4}{\left(4\pi\right)^3}g_s^2\xi_i\left[P_{{\cal I}_{j\oplus \bar{i}}{\cal I}_j}^{<}(1-\xi_i,0){\mathbb{M}}^{(n-1,0)}(\dot{r}^{j\oplus \bar{i},i\backslash})+\Delta_{ij}\right]\nonumber\\ &&\times\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r})}{\cal G}(\dot{r})\, d\varphi_id\phi_{n-1}(\dot{r}^{j\oplus \bar{i},i\backslash}), \end{align}\] and, for the final state collinear singularities (\(3\leq j \leq n_{L}^{(R)}+2, j\neq i\)), we have \[\begin{align} \Sigma_{ij}(\dot{r}; \xi_i, 1)&=&\frac{4}{\left(4\pi\right)^3}g_s^2\frac{1-z}{z}h(z)\left[P_{{\cal I}_j{\cal I}_{j\oplus i}}^{<}(z,0){\mathbb{M}}^{(n-1,0)}(\dot{r}^{j\oplus i,i\backslash})+\Delta_{ij}\right]\nonumber\\ &&\times\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r})}{\cal G}(\dot{r})\, d\varphi_id\phi_{n-1}(\dot{r}^{j\oplus i,i\backslash})\,. \end{align}\] For the soft-collinear counterparts, we need to take \(\xi_i\rightarrow 0\) and \(z\rightarrow 1\) on the r.h.s of the above two equations. The soft local counterterm is \[\begin{align} \Sigma_{ij}(\dot{r}; 0, y_{ij})&=&\frac{s}{\left(4\pi\right)^3}(1-y_{ij})\Big(\lim_{\xi_i\rightarrow 0}{\xi_i^2{\mathbb{M}}^{(n,0)}(\dot{r})}\Big) c_{ij}\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r})}{\cal G}(\dot{r})\, d\varphi_id\phi_{n-1}(\dot{r}^{i\backslash}), \end{align}\] where the soft limit of the real matrix element has been derived in sect. 3.2.

4.3 Integrated soft counterterm↩︎

We can split the integrated soft counterterms into two parts: \[\begin{align} d\hat{\sigma}^{(S)}(\dot{r})&=&d\hat{\sigma}^{(S,1)}(\dot{r})+d\hat{\sigma}^{(S,2)}(\dot{r}). \end{align}\]

The first one is simply stemming from the \((4-2\epsilon)\)-dimensional counterpart of the function \(\Sigma_{ij}(\dot{r}; 0, y_{ij})\) in the local subtraction counterterms appearing in the former section. It is \[\begin{align} &&d\hat{\sigma}^{(S,1)}(\dot{r})\nonumber\\ &=&d\phi_{n-1}(\dot{r}^{i\backslash})\left(\frac{\mu^{2}}{s}\right)^\epsilon\frac{s}{(4\pi)^{3-2\epsilon}}\sum_{(i,j)\in{\cal P}_{\rm FKS}(\dot{r})}{c_{ij}\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r})}{\cal G}(\dot{r})}\nonumber\\ &&\times\int{\xi_i^{-1-2\epsilon}(1-y_{ij}^2)^{-\epsilon}\left(\lim_{\xi_i\rightarrow 0}{\xi_i^2{\mathbb{M}}^{(n,0)}(\dot{r})}\right)\Theta(\xi_{cut}-\xi_i)d\xi_idy_{ij}d\Omega_i^{(2-2\epsilon)}}\nonumber\\ &=&\frac{d\phi_{n-1}(\dot{r}^{i\backslash})}{8\pi^2}\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r})}{\cal G}(\dot{r})\mathop{\sum_{i=3}}_{{\cal I}_i=g}^{n_{L}^{(R)}+2}{\left[-\frac{\xi_{cut}^{-2\epsilon}}{2\epsilon}\frac{2^{2\epsilon}}{(2\pi)^{1-2\epsilon}}\left(\frac{s}{\mu^2}\right)^{-\epsilon}\int{\left(\lim_{\xi_i\rightarrow 0}{\left(\tilde{k}_i^0\right)^2{\mathbb{M}}^{(n,0)}(\dot{r})}\right)d\Omega_i}\right]}\nonumber\\ &=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r})}{\cal G}(\dot{r})\nonumber\\ &&\times\mathop{\sum_{i=3}}_{{\cal I}_i=g}^{n_{L}^{(R)}+2}{\mathop{\sum_{k,l=n_{I}}}_{k,l\neq i, k\leq l}^{n_{L}^{(R)}+n_{H}+1}{\left[-\frac{\xi_{cut}^{-2\epsilon}}{2\epsilon}\frac{2^{2\epsilon}}{(2\pi)^{1-2\epsilon}}\left(\frac{s}{\mu^2}\right)^{-\epsilon}\int{\left(\lim_{\xi_i\rightarrow 0}{\left(\tilde{k}_i^0\right)^2{\mathsf M}^{(n-1,0)}_{{\rm soft},kl}(\dot{r}^{i\backslash})}\right)d\Omega_i}\right]}}\nonumber\\ &=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\nonumber\\ &&\times\sum_{k=n_{I}}^{n_{L}^{(B)}+n_{H}+1}{\sum_{l=k}^{n_{L}^{(B)}+n_{H}+1}{\left[-\frac{\xi_{cut}^{-2\epsilon}}{2\epsilon}\frac{2^{2\epsilon}}{(2\pi)^{1-2\epsilon}}\left(\frac{s}{\mu^2}\right)^{-\epsilon}\int{\left(\lim_{\xi_i\rightarrow 0}{\left(\tilde{k}_i^0\right)^2{\mathsf M}^{(n-1,0)}_{{\rm soft},kl}(\dot{r}^{i\backslash})}\right)d\Omega_i}\right]}}\,,\nonumber\\ \end{align}\] where the \(d-1=3-2\epsilon\) dimensional solid angle measure is \[\begin{align} d\Omega_i&=&(1-y_{ij}^2)^{-\epsilon}dy_{ij}d\Omega_i^{(2-2\epsilon)},\label{eq:solidanglemeasure} \end{align}\tag{47}\] and we have used eq.@eq:eq:softampsq0 as well as the fact that \(\lim_{\xi_i\rightarrow 0}{\xi_i^2{\mathbb{M}}^{(n,0)}(\dot{r})}\) is independent of \(\xi_i\). In the last equation, we have removed any final state gluon \({\cal I}_i=g\) and have changed the final state symmetry factor accordingly.

The second one is from the soft singularities of the LDMEs in NRQCD, which we will give explicitly in the following.

4.3.1 Soft counterterms of LDMEs↩︎

The second term originates from the renormalisation of LDMEs and is analogous to the initial state collinear counterterm at hadron colliders or the final state collinear counterterm due to the presence of fragmentation functions [2]. The soft counterterms of the LDMEs can be incorporated into either the real or the virtual matrix elements, depending on the preference. In this paper, we include them as part of the integrated soft subtraction terms for the real matrix element. We have the following perturbative corrections for the LDMEs [63] in the \(\overline{{\rm MS}}\) scheme, by defining \(1/\bar{\epsilon}=1/\epsilon+\log{\left(4\pi\right)}-\gamma_E\), \[\begin{align} \braket{{\mathcal{O}}^H_{\bigl.^{3}S^{[8]}_1}}(\mu)&=&\braket{{\mathcal{O}}^H_{\bigl.^{3}S^{[8]}_1}}+\frac{4\alpha_s}{3\pi m_{Q}m_{\bar{Q}^\prime}}\left(\frac{1}{\bar{\epsilon}}+\log{\frac{\mu^2}{\mu_{{\rm NRQCD}}^2}}\right)\!\left[\frac{C_{F}}{2N_{c}}\!\sum_{J=0}^{2}{\braket{{\mathcal{O}}^H_{\bigl.^{3}P^{[1]}_J}}}\!+\!B_{F}\!\sum_{J=0}^{2}{\braket{{\mathcal{O}}^H_{\bigl.^{3}P^{[8]}_J}}}\right],\nonumber\\ \braket{{\mathcal{O}}^H_{\bigl.^{3}S^{[1]}_1}}(\mu)&=&\braket{{\mathcal{O}}^H_{\bigl.^{3}S^{[1]}_1}}+\frac{4\alpha_s}{3\pi m_{Q}m_{\bar{Q}^\prime}}\left(\frac{1}{\bar{\epsilon}}+\log{\frac{\mu^2}{\mu_{{\rm NRQCD}}^2}}\right)\!\left[\sum_{J=0}^{2}{\braket{{\mathcal{O}}^H_{\bigl.^{3}P^{[8]}_J}}}\right],\nonumber\\ \braket{{\mathcal{O}}^H_{\bigl.^{1}S^{[8]}_0}}(\mu)&=&\braket{{\mathcal{O}}^H_{\bigl.^{1}S^{[8]}_0}}+\frac{4\alpha_s}{3\pi m_{Q}m_{\bar{Q}^\prime}}\left(\frac{1}{\bar{\epsilon}}+\log{\frac{\mu^2}{\mu_{{\rm NRQCD}}^2}}\right)\!\left[\frac{C_{F}}{2N_{c}}\braket{{\mathcal{O}}^H_{\bigl.^{1}P^{[1]}_1}}+B_{F}\braket{{\mathcal{O}}^H_{\bigl.^{1}P^{[8]}_1}}\right],\nonumber\\ \braket{{\mathcal{O}}^H_{\bigl.^{1}S^{[1]}_0}}(\mu)&=&\braket{{\mathcal{O}}^H_{\bigl.^{1}S^{[1]}_0}}+\frac{4\alpha_s}{3\pi m_{Q}m_{\bar{Q}^\prime}}\left(\frac{1}{\bar{\epsilon}}+\log{\frac{\mu^2}{\mu_{{\rm NRQCD}}^2}}\right)\!\left[\braket{{\mathcal{O}}^H_{\bigl.^{1}P^{[8]}_1}}\right],\label{eq:NRQCDCTs} \end{align}\tag{48}\] where \(\gamma_E\) is the Euler-Mascheroni constant, \(\mu_{{\rm NRQCD}}\) is the NRQCD cutoff scale and \(B_{F}=(N_{c}^2-4)/(4N_{c})=5/12\). If we use the Binoth Les Houches Accord convention by taking out the global prefactor \((4\pi)^\epsilon/\Gamma(1-\epsilon)\), \(1/\bar{\epsilon}\) simply becomes \(1/\epsilon\) in eq.@eq:eq:NRQCDCTs . Therefore, we have the additional integrated soft counterterms for P-waves:

• \(Q\bar{Q}^\prime[\bigl.^{1}P^{[1]}_1]\): \[\begin{align} d\hat{\sigma}^{(S,2)}(\dot{r})&=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\left[\frac{1}{2N_{c}}{\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash}_1)\right]\nonumber\\ &&\times\frac{8}{m_{Q}m_{\bar{Q}^\prime}}\left(\frac{1}{\bar{\epsilon}}+\log{\frac{\mu^2}{\mu_{{\rm NRQCD}}^2}}\right)\,, \end{align}\] where \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{1}P^{[1]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_1^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{1}S^{[8]}_0],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right), \end{align}\] and \({\cal I}_i=g\).

• \(Q\bar{Q}^\prime[\bigl.^{1}P^{[8]}_1]\): \[\begin{align} d\hat{\sigma}^{(S,2)}(\dot{r})&=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\left[B_{F}{\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash}_1)+C_{F}{\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash}_2)\right]\nonumber\\ &&\times\frac{8}{m_{Q}m_{\bar{Q}^\prime}}\left(\frac{1}{\bar{\epsilon}}+\log{\frac{\mu^2}{\mu_{{\rm NRQCD}}^2}}\right), \end{align}\] where \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{1}P^{[8]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_1^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{1}S^{[8]}_0],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_2^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{1}S^{[1]}_0],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right). \end{align}\]

• \(Q\bar{Q}^\prime[\bigl.^{3}P^{[1]}_J]\): \[\begin{align} d\hat{\sigma}^{(S,2)}(\dot{r})&=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\left[\frac{2J+1}{2N_{c}}{\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash}_1)\right]\nonumber\\ &&\times\frac{8}{9 m_{Q}m_{\bar{Q}^\prime}}\left(\frac{1}{\bar{\epsilon}}+\log{\frac{\mu^2}{\mu_{{\rm NRQCD}}^2}}\right)\,, \end{align}\] where \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{3}P^{[1]}_J],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_1^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{3}S^{[8]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right). \end{align}\]

• \(Q\bar{Q}^\prime[\bigl.^{3}P^{[8]}_J]\): \[\begin{align} d\hat{\sigma}^{(S,2)}(\dot{r})&=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\left[B_{F}{\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash}_1)+C_{F}{\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash}_2)\right]\nonumber\\ &&\times\frac{8(2J+1)}{9 m_{Q}m_{\bar{Q}^\prime}}\left(\frac{1}{\bar{\epsilon}}+\log{\frac{\mu^2}{\mu_{{\rm NRQCD}}^2}}\right), \end{align}\] where \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{3}P^{[8]}_J],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_1^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{3}S^{[8]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_2^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,{\cal I}_{n_{L}^{(R)}+n_{H}},Q\bar{Q}^\prime[\bigl.^{3}S^{[1]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right). \end{align}\]

It is clear that for S-waves \(d\hat{\sigma}^{(S,2)}(\dot{r})=0\).

4.3.2 Colour-singlet S-wave state↩︎

The soft integrated partonic cross section for a colour-singlet S-wave quarkonium production is \[\begin{align} d\hat{\sigma}^{(S)}(\dot{r})&=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\sum_{k=n_{I}}^{n_{L}^{(B)}+n_{H}}{\sum_{l= k}^{n_{L}^{(B)}+n_{H}}{\bar{\mathcal{E}}(\{1,1\},\{k_k,k_l\}){\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash})}}\,,\nonumber\\ \end{align}\] where the eikonal integrals \(\bar{\mathcal{E}}()\) can be found in app. 7, and \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[1]}_J],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right)\label{eq:BornCSSwave} \end{align}\tag{49}\] with \({\cal I}_i\) being a final state gluon. Now, the momenta \(k_k, k_l\) are defined in the reduced Born process \(\dot{r}^{i\backslash}\in \dot{{\cal R}}_{n-1}\) with \(k_{n_{L}^{(B)}+n_{H}+1}=K\) (the quarkonium momentum).

4.3.3 Colour-octet S-wave state↩︎

The integrated soft counterterm for a colour-octet S-wave quarkonium is \[\begin{align} d\hat{\sigma}^{(S)}(\dot{r})&=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\sum_{k=n_{I}}^{n_{L}^{(B)}+n_{H}+1}{\sum_{l= k}^{n_{L}^{(B)}+n_{H}+1}{\bar{\mathcal{E}}(\{1,1\},\{k_k, k_l\}){\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash})}}\,,\nonumber\\ \end{align}\] where the eikonal integrals can be found in app. 7, and \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,Q\bar{Q}^\prime[\bigl.^{2S+1}S^{[8]}_J],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right)\label{eq:BornCOSwave} \end{align}\tag{50}\] with \({\cal I}_i\) being a gluon. The momenta \(k_k, k_l\) are defined in the Born process \(\dot{r}^{i\backslash}\in \dot{{\cal R}}_{n-1}\) with \(k_{n_{L}^{(B)}+n_{H}+1}=K\) (the quarkonium momentum).

4.3.4 Colour-singlet spin-singlet P-wave state↩︎

For a colour-singlet spin-singlet P-wave state, the integrated soft counterterm is \[\begin{align} d\hat{\sigma}^{(S)}(\dot{r})&=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\Bigg[\sum_{k=n_{I}}^{n_{L}^{(B)}+n_{H}}{\sum_{l= k}^{n_{L}^{(B)}+n_{H}}{\bar{\mathcal{E}}(\{1,1\},\{k_k, k_l\}){\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash})}}\nonumber\\ &&+\sum_{k=n_{I}}^{n_{L}^{(B)}+n_{H}}{\left(\bar{\mathcal{E}}(\{1,1\},\{k_k,K\})\frac{k_{k,\mu}}{K\cdot k_k}-\bar{\mathcal{E}}_{\mu}(\{1,2\},\{k_k,K\})\right){\mathbb{M}}^{(n-1,0)}_{k[18]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu}\nonumber\\ &&+\left(\frac{1}{2N_{c}}\frac{8}{m_{Q}m_{\bar{Q}^\prime}}\left(\frac{1}{\bar{\epsilon}}+\log{\frac{\mu^2}{\mu^2_{\rm NRQCD}}}\right)-\frac{2\epsilon-2}{K^2}\bar{\mathcal{E}}(\{1,1\},\{K,K\})C_{{\rm eff}}(Q\bar{Q}^\prime_{[18]})\right)\nonumber\\ &&\times{\mathbb{M}}^{(n-1,0)}(\dot{r}_1^{i\backslash})\Bigg]\,, \end{align}\] where the eikonal integrals can be found in app. 7, and \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,Q\bar{Q}^\prime[\bigl.^{1}P^{[1]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_1^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,Q\bar{Q}^\prime[\bigl.^{1}S^{[8]}_0],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\label{eq:BornCS1Pwave} \end{align}\tag{51}\] with \({\cal I}_i=g\). Once again, we remind the readers that the momenta \(k_k, k_l\) are defined in the reduced Born process \(\dot{r}^{i\backslash}\in \dot{{\cal R}}_{n-1}\) with the quarkonium momentum \(k_{n_{L}^{(B)}+n_{H}+1}=K\). Now, besides the usual scalar eikonal integrals, we also need to introduce the rank-\(1\) tensor integrals \(\bar{\mathcal{E}}_{\mu}(\{1,2\},\{k_k,K\})\).

4.3.5 Colour-octet spin-singlet P-wave state↩︎

The soft integrated counterterm for a colour-octet spin-singlet P-wave quarkonium is \[\begin{align} d\hat{\sigma}^{(S)}(\dot{r})&=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\left[\sum_{k=n_{I}}^{n_{L}^{(B)}+n_{H}+1}{\sum_{l= k}^{n_{L}^{(B)}+n_{H}+1}{\bar{\mathcal{E}}(\{1,1\},\{k_k, k_l\}){\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash})}}\right.\nonumber\\ &&+\sum_{k=n_{I}}^{n_{L}^{(B)}+n_{H}+1}{\left(\bar{\mathcal{E}}(\{1,1\},\{k_k,K\})\frac{k_{k,\mu}}{K\cdot k_k}-\bar{\mathcal{E}}_{\mu}(\{1,2\},\{k_k,K\})\right)}\nonumber\\ &&\times\left({\mathbb{M}}^{(n-1,0)}_{k[88]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu+{\mathbb{M}}^{(n-1,0)}_{k[81]}(\dot{r}^{i\backslash},\dot{r}_2^{i\backslash})^\mu\right)\nonumber\\ &&+\left(B_{F}\frac{8}{m_{Q}m_{\bar{Q}^\prime}}\left(\frac{1}{\bar{\epsilon}}+\log{\frac{\mu^2}{\mu^2_{\rm NRQCD}}}\right)-\frac{2\epsilon-2}{K^2}\bar{\mathcal{E}}(\{1,1\},\{K,K\})C_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})\right)\nonumber\\ &&\times{\mathbb{M}}^{(n-1,0)}(\dot{r}_1^{i\backslash})\nonumber\\ &&+\left(C_{F}\frac{8}{m_{Q}m_{\bar{Q}^\prime}}\left(\frac{1}{\bar{\epsilon}}+\log{\frac{\mu^2}{\mu^2_{\rm NRQCD}}}\right)-\frac{2\epsilon-2}{K^2}\bar{\mathcal{E}}(\{1,1\},\{K,K\})C_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})\right)\nonumber\\ &&\times{\mathbb{M}}^{(n-1,0)}(\dot{r}_2^{i\backslash})\Bigg]\,, \end{align}\] where the expressions of the eikonal integrals can be found in app. 7, and \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,Q\bar{Q}^\prime[\bigl.^{1}P^{[8]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_1^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,Q\bar{Q}^\prime[\bigl.^{1}S^{[8]}_0],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_2^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,Q\bar{Q}^\prime[\bigl.^{1}S^{[1]}_0],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\label{eq:BornCO1Pwave} \end{align}\tag{52}\] and \({\cal I}_i\) is a final state gluon.

4.3.6 Colour-singlet spin-triplet P-wave state↩︎

The integrated soft partonic cross section of a process with a colour-singlet spin-triplet P-wave quarkonium is \[\begin{align} d\hat{\sigma}^{(S)}(\dot{r})&=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\left[\sum_{k=n_{I}}^{n_{L}^{(B)}+n_{H}}{\sum_{l= k}^{n_{L}^{(B)}+n_{H}}{\bar{\mathcal{E}}(\{1,1\},\{k_k, k_l\}){\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash})}}\right.\nonumber\\ &&+\sum_{k=n_{I}}^{n_{L}^{(B)}+n_{H}}{\left(\bar{\mathcal{E}}(\{1,1\},\{k_k,K\})\frac{k_{k,\mu}}{K\cdot k_k}-\bar{\mathcal{E}}_{\mu}(\{1,2\},\{k_k,K\})\right){\mathbb{M}}^{(n-1,0)}_{J,k[18]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu}\nonumber\\ &&-\left(\frac{g_{\mu\nu}}{K^2}\bar{\mathcal{E}}(\{1,1\},\{K,K\})+\bar{\mathcal{E}}_{\mu\nu}(\{2,2\},\{K,K\})\right)C_{{\rm eff}}(Q\bar{Q}^\prime_{[18]}){\mathbb{M}}^{(n-1,0)}_J(\dot{r}_1^{i\backslash})^{\mu\nu}\nonumber\\ &&\left.+\frac{2J+1}{2N_{c}}\frac{8}{9m_{Q}m_{\bar{Q}^\prime}}\left(\frac{1}{\bar{\epsilon}}+\log{\frac{\mu^2}{\mu^2_{\rm NRQCD}}}\right){\mathbb{M}}^{(n-1,0)}(\dot{r}_1^{i\backslash})\right], \end{align}\] where the eikonal integrals can be found in app. 7, and \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,Q\bar{Q}^\prime[\bigl.^{3}P^{[1]}_J],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_1^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,Q\bar{Q}^\prime[\bigl.^{3}S^{[8]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\label{eq:BornCS3Pwave} \end{align}\tag{53}\] and \({\cal I}_i\) is a final state gluon. In this case, we also need the rank-\(2\) eikonal tensor integrals \(\bar{\mathcal{E}}_{\mu\nu}(\{2,2\},\{K,K\})\).

4.3.7 Colour-octet spin-triplet P-wave state↩︎

Finally, the colour-octet spin-triplet P-wave quarkonium production has the following soft integrated counterterm \[\begin{align} d\hat{\sigma}^{(S)}(\dot{r})&=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\left[\sum_{k=n_{I}}^{n_{L}^{(B)}+n_{H}+1}{\sum_{l= k}^{n_{L}^{(B)}+n_{H}+1}{\bar{\mathcal{E}}(\{1,1\},\{k_k, k_l\}){\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash})}}\right.\nonumber\\ &&+\sum_{k=n_{I}}^{n_{L}^{(B)}+n_{H}+1}{\left(\bar{\mathcal{E}}(\{1,1\},\{k_k,K\})\frac{k_{k,\mu}}{K\cdot k_k}-\bar{\mathcal{E}}_{\mu}(\{1,2\},\{k_k,K\})\right)}\nonumber\\ &&\times\left({\mathbb{M}}^{(n-1,0)}_{J,k[88]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu+{\mathbb{M}}^{(n-1,0)}_{J,k[81]}(\dot{r}^{i\backslash},\dot{r}_2^{i\backslash})^\mu\right)\nonumber\\ &&-\left(\frac{g_{\mu\nu}}{K^2}\bar{\mathcal{E}}(\{1,1\},\{K,K\})+\bar{\mathcal{E}}_{\mu\nu}(\{2,2\},\{K,K\})\right)\nonumber\\ &&\times\left(C_{{\rm eff}}(Q\bar{Q}^\prime_{[88]}){\mathbb{M}}^{(n-1,0)}_J(\dot{r}_1^{i\backslash})^{\mu\nu}+C_{{\rm eff}}(Q\bar{Q}^\prime_{[81]}){\mathbb{M}}^{(n-1,0)}_J(\dot{r}_2^{i\backslash})^{\mu\nu}\right)\nonumber\\ &&\left.+\frac{8(2J+1)}{9m_{Q}m_{\bar{Q}^\prime}}\!\left(\frac{1}{\bar{\epsilon}}+\log{\frac{\mu^2}{\mu^2_{\rm NRQCD}}}\right)\!\!\left(B_{F}{\mathbb{M}}^{(n-1,0)}(\dot{r}_1^{i\backslash})+C_{F}{\mathbb{M}}^{(n-1,0)}(\dot{r}_2^{i\backslash})\right)\right], \end{align}\] where the eikonal integrals can be found in app. 7, and \[\begin{align} \dot{r}^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,Q\bar{Q}^\prime[\bigl.^{3}P^{[8]}_J],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_1^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,Q\bar{Q}^\prime[\bigl.^{3}S^{[8]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\nonumber\\ \dot{r}_2^{i\backslash}&=&\left({\cal I}_1,\ldots,{\cal I}\backslash_i,\ldots,{\cal I}_j,\ldots,Q\bar{Q}^\prime[\bigl.^{3}S^{[1]}_1],\bar{Q}^\prime\backslash, \ldots,{\cal I}_{n+3}\right),\label{eq:BornCO3Pwave} \end{align}\tag{54}\] with \({\cal I}_i\) being a final state gluon.

4.4 Integrated collinear and soft-collinear counterterms↩︎

As explained earlier, the integrated collinear and soft-collinear counterterms should be identical to the original case since quarkonia are massive and do not exhibit any collinear divergences. However, for the completeness, we still discuss them here and adjust to our notations with the risk that they are actually well understood.

Following the explanation of ref. [50], it is straightforward to show that the sum of the integrated counterterms for collinear and soft-collinear emissions from an incoming initial state parton takes the form \[\begin{align} d\hat{\sigma}^{(in)}(\dot{r})&=&-\frac{\alpha_s}{2\pi}\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\sum_{k=n_{I}}^{2}\int_0^{\xi_{\rm max}} d\xi_i\left(\frac{1}{\bar{\epsilon}}-\log\!\left(\frac{s \delta_{I}}{2\mu^2}\right)\right)\left[\left(\frac{1}{\xi_i}\right)_c-2\epsilon\left(\frac{\log(\xi_i)}{\xi_i}\right)_c\right]\nonumber\\ &&\times\xi_iP_{{\cal I}_{k\oplus \bar{i}}{\cal I}_k}^{<}(1-\xi_i,\epsilon){\mathbb{M}}^{(n-1,0)}(\dot{r}^{k\oplus \bar{i},i\backslash})d\phi_{n-1}(\dot{r}^{k\oplus \bar{i},i\backslash})\,, \label{eq:dsiginitial} \end{align}\tag{55}\] and that for emissions from an outgoing final state parton for a given underlying Born process \(\dot{r}^{i\backslash}\) the counterterms can be written as \[\begin{align} d\hat{\sigma}^{(out)}(\dot{r})&=&\frac{\alpha_s}{2\pi}\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\sum_{k=3}^{n_{L}^{(B)}+2}\bigg\lbrace\frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)} \left(\frac{\mu^2}{Q_{\rm ES}^2}\right)^\epsilon\frac{1}{\epsilon}\left[\gamma({\cal I}_k)+C({\cal I}_k)\log\!\left(\frac{\xi_{cut}^2s}{4E_k^2}\right)\right]\nonumber\\ &&+\left[\gamma^\prime({\cal I}_k)-\log\!\left(\displaystyle\frac{s\delta_{O}}{2Q_{\rm ES}^2}\right)\left(\gamma({\cal I}_k)-2C({\cal I}_k)\log\!\left(\displaystyle\frac{2E_k}{\xi_{cut}\sqrt{s}}\right)\right)\right.\nonumber\\ &&\left.+2C({\cal I}_k)\left(\log^2\!\left(\displaystyle\frac{2E_k}{\sqrt{s}}\right)-\log^2(\xi_{cut})\right)-2\gamma({\cal I}_k)\log\!\left(\displaystyle\frac{2E_k}{\sqrt{s}}\right)\right]\bigg\rbrace\nonumber\\ &&\times{\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash})d\phi_{n-1}(\dot{r}^{i\backslash})\,, \label{eq:dsigfinal} \end{align}\tag{56}\] where \(Q_{\rm ES}\) is the Ellis-Sexton scale [73], the Casimir factors are \[\begin{align} C({\cal I})&=\left\{\begin{array}{ll} C_F, &~~{\rm if}~~{\cal I}\in \mathbf{3}^{\count 0=0 \loop \ifnum\count 0>0 \advance\count 0 by -1 \prime\repeat}, \settoheight{\irrepbarheight}{\mathbf{3}} \settowidth{\irrepwidth}{\mathbf{3}} \makebox[0pt][l]{\mathbf{3}} \rule[1.2\irrepbarheight]{\irrepwidth}{\irrepbarthickness}^{\count 0=0 \loop \ifnum\count 0>0 \advance\count 0 by -1 \prime\repeat}\\ C_A, &~~{\rm if}~~{\cal I}\in \mathbf{8}^{\count 0=0 \loop \ifnum\count 0>0 \advance\count 0 by -1 \prime\repeat}\\ \end{array}\right.. \end{align}\] The collinear anomalous dimensions are \[\begin{align} \gamma({\cal I})&=\left\{\begin{array}{ll} \frac{3}{2}C_F\,,\hfill&~~{\rm if}~~{\cal I}=q,\bar{q}\\ \frac{11}{6}C_A-\frac{2}{3}T_Fn_f\,,&~~{\rm if}~~{\cal I}=g\\ \end{array}\right., \end{align}\] with \(T_F=\frac{1}{2}\) and \(n_f\) being the number of massless quark flavours, and \[\begin{align} \gamma^\prime({\cal I})&=\left\{\begin{array}{ll} \left(\frac{13}{2}-\frac{2\pi^2}{3}\right)C_F\,,\hfill&~~{\rm if}~~{\cal I}=q,\bar{q}\\ \left(\frac{67}{9}-\frac{2\pi^2}{3}\right)C_A-\frac{23}{9}T_Fn_f\,,&~~{\rm if}~~{\cal I}=g\\ \end{array}\right.. \end{align}\] These integrated counterterms have to be combined with the counterterms stemming from the collinear renormalisation of PDFs which read \[\begin{align} d\hat{\sigma}^{(cnt)}(\dot{r})&=&\frac{\alpha_s}{2\pi}\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\sum_{k=n_{I}}^{2}\int_{1-\xi_{\rm max}}^1 dz \left(\frac{1}{\bar{\epsilon}}P_{{\cal I}_{k\oplus \bar{i}}{\cal I}_k}(z,0)-K_{{\cal I}_{k\oplus \bar{i}}{\cal I}_k}(z)\right)\nonumber\\ &&\times{\mathbb{M}}^{(n-1,0)}(\dot{r}^{k\oplus \bar{i},i\backslash})d\phi_{n-1}(\dot{r}^{k\oplus \bar{i},i\backslash}) \end{align}\] with the four dimensional regularised Altarelli-Parisi splitting kernels \[P_{ab}(z,0) = \frac{(1-z)P_{ab}^{<}(z,0)}{(1-z)_{+}} + \gamma(a)\delta_{ab}\,\delta(1-z)\] and the subtraction scheme dependent function \(K_{ab}(z)\) ⁸. A change of the integration variable, \(z\rightarrow 1-\xi_i\), results in ⁹ \[\begin{align} d\hat{\sigma}^{(cnt)}(\dot{r})&=&\frac{\alpha_s}{2\pi}\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\sum_{k=n_{I}}^{2}\left\lbrace\frac{1}{\bar{\epsilon}}\left[\gamma({\cal I}_k)+C({\cal I}_k)\log\!\left(\frac{\xi_{cut}^2s}{4 E_k^2}\right)\right]{\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash})d\phi_{n-1}(\dot{r}^{i\backslash})\right.\nonumber\\ &&+\int_0^{\xi_{\rm max}} d\xi_i\left[\frac{1}{\bar{\epsilon}}\left(\frac{1}{\xi_i}\right)_c\xi_iP_{{\cal I}_{k\oplus \bar{i}}{\cal I}_k}^{<}(1-\xi_i,0)-K_{{\cal I}_{k\oplus \bar{i}}{\cal I}_k}(1-\xi_i)\right]{\mathbb{M}}^{(n-1,0)}(\dot{r}^{k\oplus \bar{i},i\backslash})\nonumber\\ &&\left.\times d\phi_{n-1}(\dot{r}^{k\oplus \bar{i},i\backslash})\right\rbrace\,, \label{eq:dsigpdf} \end{align}\tag{57}\] where we made use of the identity \[\delta(\xi_i)\xi_iP_{ab}^{<}(1-\xi_i,0) = 2C(a)\delta_{ab}\,\delta(\xi_i)\,.\] If we sum up eqs.(55 ), (56 ), and (57 ), we obtain the sum of the integrated collinear and soft-collinear counterterms \[\begin{align} d\hat{\sigma}^{(C)}(\dot{r})&=&d\hat{\sigma}^{(in)}(\dot{r})+d\hat{\sigma}^{(out)}(\dot{r})+d\hat{\sigma}^{(cnt)}(\dot{r})\nonumber\\ &=&d\hat{\sigma}^{(C)}_{\rm sing.}(\dot{r})+d\hat{\sigma}^{(C,n-1)}_{\rm FIN}(\dot{r})+d\hat{\sigma}^{(C,n)}_{\rm FIN}(\dot{r})\,, \end{align}\] which we have decomposed into a singular term \[\begin{align} d\hat{\sigma}^{(C)}_{\rm sing.}(\dot{r})&=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\sum_{k=n_{I}}^{n_{L}^{(B)}+2}\frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)} \left(\frac{\mu^2}{Q_{\rm ES}^2}\right)^\epsilon\frac{1}{\epsilon}\nonumber\\ &&\times\left[\gamma({\cal I}_k)+C({\cal I}_k)\log\!\left(\frac{\xi_{cut}^2s}{4E_k^2}\right)\right]{\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash})\,, \end{align}\] a finite \((n-1)\)-body contribution (equivalent to eqs.(4.5-4.6) in ref. [56])

\[\begin{align} d\hat{\sigma}^{(C,n-1)}_{\rm FIN}(\dot{r})&=&\frac{\alpha_s}{2\pi}d\phi_{n-1}(\dot{r}^{i\backslash})\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\Bigg\lbrace-\log\!\left(\displaystyle\frac{\mu^2}{Q_{\rm ES}^2}\right)\sum_{k=n_{I}}^{2}\Big(\gamma({\cal I}_k)+2C({\cal I}_k)\log(\xi_{cut})\Big)\nonumber\\ &&+\sum_{k= 3}^{n_{L}^{(B)}+2}\left[\gamma^\prime({\cal I}_k)-\log\!\left(\displaystyle\frac{s\delta_{O}}{2Q_{\rm ES}^2}\right)\left(\gamma({\cal I}_k)-2C({\cal I}_k)\log\!\left(\displaystyle\frac{2E_k}{\xi_{cut}\sqrt{s}}\right)\right)\right.\nonumber\\ &&\left.+2C({\cal I}_k)\left(\log^2\!\left(\displaystyle\frac{2E_k}{\sqrt{s}}\right)-\log^2(\xi_{cut})\right)-2\gamma({\cal I}_k)\log\!\left(\displaystyle\frac{2E_k}{\sqrt{s}}\right)\right]\Bigg\rbrace{\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash})\,,\nonumber\\ \end{align}\] and a finite contribution with a degenerated \(n\)-body phase space (equivalent to eqs.(4.40-4.42) in ref. [56]), where the angular dependence of the collinear emitted parton is integrated out, but its energy integration remains left, \[\begin{align} d\hat{\sigma}^{(C,n)}_{\rm FIN}(\dot{r})&=&\frac{\alpha_s}{2\pi}\frac{J^{n_{L}^{(B)}}}{{\cal N}(\dot{r}^{i\backslash})}{\cal G}(\dot{r}^{i\backslash})\sum_{k=n_{I}}^{2}\int_0^{\xi_{\rm max}} d\xi_i\left\lbrace\left[\left(\frac{1}{\xi_i}\right)_c\log\!\left(\frac{s\delta_{I}}{2\mu^2}\right)+2\left(\frac{\log(\xi_i)}{\xi_i}\right)_c\right]\right.\nonumber\\ &&\left.\times\xi_iP_{{\cal I}_{k\oplus \bar{i}}{\cal I}_k}^{<}(1-\xi_i,0)-\left(\frac{1}{\xi_i}\right)_c\xi_iP_{{\cal I}_{k\oplus \bar{i}}{\cal I}_k}^{\,\prime<}(1-\xi_i,0)-K_{{\cal I}_{k\oplus \bar{i}}{\cal I}_k}(1-\xi_i)\right\rbrace\nonumber\\ &&\times{\mathbb{M}}^{(n-1,0)}(\dot{r}^{k\oplus \bar{i},i\backslash})d\phi_{n-1}(\dot{r}^{k\oplus \bar{i},i\backslash})\,. \end{align}\] \(P_{ab}^{\,\prime<}(z,0)\) is the \(\epsilon\) part of the \(d=4-2\epsilon\) dimensional (unregularised) Altarelli-Parisi splitting kernels \[\begin{align} P_{ab}^{<}(z,\epsilon)&=&P_{ab}^{<}(z,0)+\epsilon P_{ab}^{\,\prime<}(z,0). \end{align}\]

5 Validations↩︎

We present a few cross-checks of our new formalism derived in the last sections. Concerning the local FKS counterterms, we have verified the initial collinear and soft limits of all \(2\rightarrow 1\) partonic Born processes \(gg\rightarrow c\bar{c}[n]\) and \(q\bar{q}\rightarrow c\bar{c}[n]\) with \(n\) covering all S- and P-wave charmonium states, including both colour-singlet and colour-octet. One of the most non-trivial checks we performed involves the \(2\rightarrow 2\) Born processes illustrated in fig. 1. We consider the real processes \(gg\rightarrow c\bar{c}[n]gg\), \(q\bar{q}\rightarrow c\bar{c}[n]gg\) and \(bg\rightarrow c\bar{c}[n]bg\), where both the charm quark \(c\) and the bottom quark \(b\) are massive, and the light quark \(q\) is massless. The real matrix elements are numerically evaluated with HELAC-Onia [9], [10], while the (reduced) Born and colour-linked Born matrix elements are computed analytically. The relative difference of the matrix elements is defined as \[\begin{align} \Delta_{\rm soft}&=&\left|\displaystyle\frac{{\mathbb{M}}^{(n,0)}(\dot{r})-\lim_{k_i\rightarrow 0}{{\mathbb{M}}^{(n,0)}(\dot{r})}}{{\mathbb{M}}^{(n,0)}(\dot{r})}\right|\propto\xi_i+\mathcal{O}(\xi_i^2), \end{align}\] where the expression of \(\lim_{k_i\rightarrow 0}{{\mathbb{M}}^{(n,0)}(\dot{r})}\) corresponds to the local soft counterterm of the real-emission counterpart \({\mathbb{M}}^{(n,0)}(\dot{r})\) defined on the r.h.s of eq.@eq:eq:softampsq0 . In fig. 1, we demonstrate the soft gluon limit as \(\xi_5\rightarrow 0\) for all \(8\) P-wave Fock states. The relative differences \(\Delta_{\rm soft}\) linearly vanish with \(\xi_5\) (the rescaled final gluon energy defined in eq.@eq:eq:xiidef ) asymptotically approaching zero. For S-wave states, similar soft limit tests were conducted, and the same behaviour was observed. These non-trivial checks confirm the correctness of our local soft counterterms defined in sect. 3.1. Additionally, we have performed initial and final collinear tests for some single quarkonium production processes.

In order to test the integrated counterterms, we have verified the universal IR poles of the one-loop virtual matrix elements outlined in app. 8, stemming from the derived integrated counterterms due to the KLN theorem. These checks are carried out explicitly by computing the analytic UV-renormalised one-loop matrix elements for the processes \(q\bar{q}\rightarrow c\bar{c}[n]g\) with all eight P-wave and four S-wave Fock states \(c\bar{c}[n]\). The IR poles of the virtual matrix elements for these \(12\) processes cancel perfectly with the formulas given in app. 8 at the analytic level.

Finally, we consider the NLO QCD corrections to the \(8\) physical hadronic cross sections of \(2\rightarrow 1\) processes, as shown in fig. 2. The NLO partonic cross sections ¹⁰ implemented in HELAC-Oniaare derived analytically, and additionally, we have computed the NLO cross sections within the FKS subtraction approach numerically. The relative differences between them are below \(10^{-5}\), demonstrating perfect agreement within numerical errors.

6 Summary↩︎

In this paper, we have generalised the FKS subtraction formalism to processes involving S- or P-wave quarkonium and elementary particles at NLO QCD in NRQCD factorisation. Our main new results are the local soft counterterms in sect. 3.2 and the integrated soft counterterms in sect. 4.3. As a byproduct, we have also derived the universal IR poles of one-loop matrix elements in app. 8. This serves as a cornerstone to achieve the NLO automation of cross section computations of quarkonium production processes ¹¹. Our next steps would be to implement the new formulas in the MadGraph5_aMC@NLOframework, and to generalise the formulation presented here to processes involving more-than-one quarkonia. The latter case, however, involves additional subtleties, such as the breakdown of NRQCD factorisation for processes with two P-wave bound states elucidated in ref. [75].

7 Eikonal tensor integrals↩︎

From the integrated soft counterterms, we need to solve the general eikonal tensor integral \[\begin{align} \bar{\mathcal{E}}^{\alpha_1\ldots \alpha_{n_1-1} \beta_1\ldots \beta_{n_2-1}}(\{n_1,n_2\},\{k_k,k_l\})&=&-\frac{\xi_{{\rm cut}}^{-2\epsilon}}{2\epsilon}\frac{2^{2\epsilon}}{(2\pi)^{1-2\epsilon}}\left(\frac{s}{\mu^2}\right)^{-\epsilon}\left(k_i^0\right)^2k_k\cdot k_l\nonumber\\ &&\times\underbrace{\int{d\Omega_i \frac{k_i^{\alpha_1}\ldots k_i^{\alpha_{n_1-1}}k_i^{\beta_1}\ldots k_i^{\beta_{n_2-1}}}{(k_k\cdot k_i)^{n_1} (k_l\cdot k_i)^{n_2}}}}_{\equiv I^{\alpha_1\ldots \alpha_{n_1-1} \beta_1\ldots \beta_{n_2-1}}(\{n_1,n_2\},\{k_k,k_l\})},\quad n_1,n_2 \geq 1,\nonumber\\ \end{align}\] where we have \(m_k^2=k_k^2\) and \(m_l^2=k_l^2\). We have followed the conventions in ref. [56]. The energy of the parton \(i\) and the measure (in \(d-1=3-2\epsilon\) dimensions) over its angular variables \(d\Omega_i\) (cf. eq.@eq:eq:solidanglemeasure ) are defined in the center of mass frame of the colliding partons. We can rewrite the eikonal integral into the phase-space measure as \[\begin{align} \bar{\mathcal{E}}^{\alpha_1\ldots \alpha_{n_1-1} \beta_1\ldots \beta_{n_2-1}}(\{n_1,n_2\},\{k_k,k_l\})&=&8\pi^2\mu^{2\epsilon}k_k\cdot k_l\nonumber\\ &&\times\int{\frac{d^{3-2\epsilon}{\boldsymbol{k}}_i}{(2\pi)^{3-2\epsilon}2k_i^0} \frac{k_i^{\alpha_1}\ldots k_i^{\alpha_{n_1-1}}k_i^{\beta_1}\ldots k_i^{\beta_{n_2-1}}}{(k_k\cdot k_i)^{n_1} (k_l\cdot k_i)^{n_2}}}\nonumber\\ &&{\times\Theta(\xi_{{\rm cut}}-\xi_i)}\,. \end{align}\] The introduction of \(\Theta(\xi_{{\rm cut}}-\xi_i)\) spoils the Lorentz covariance.

In order to solve these tensor integrals, we are based on the observation that ¹² \[\begin{align} &&\frac{\partial }{\partial k_{k, \alpha_{n_1}}}I^{\alpha_1\ldots \alpha_{n_1-1} \beta_1\ldots \beta_{n_2-1}}(\{n_1,n_2\},\{k_k,k_l\})\notag\\ &=&-n_1I^{\alpha_1\ldots \alpha_{n_1-1} \alpha_{n_1} \beta_1\ldots \beta_{n_2-1}}(\{n_1+1,n_2\},\{k_k,k_l\})\,,\nonumber\\ &&\frac{\partial }{\partial k_{l, \beta_{n_2}}}I^{\alpha_1\ldots \alpha_{n_1-1} \beta_1\ldots \beta_{n_2-1}}(\{n_1,n_2\},\{k_k,k_l\})\notag\\ &=&-n_2 I^{\alpha_1\ldots \alpha_{n_1-1} \beta_1\ldots \beta_{n_2-1}\beta_{n_2}}(\{n_1,n_2+1\},\{k_k,k_l\})\,. \end{align}\] Then, we can derive the tensor integrals from the scalar integral via \[\begin{align} I^{\alpha_1}(\{2,1\},\{k_k,k_l\})&=&-\frac{\partial}{\partial k_{k, \alpha_1}}I(\{1,1\},\{k_k,k_l\}),\nonumber\\ I^{\beta_1}(\{1,2\},\{k_k,k_l\})&=&-\frac{\partial}{\partial k_{l, \beta_1}}I(\{1,1\},\{k_k,k_l\}),\nonumber\\ I^{\alpha_1\beta_1}(\{2,2\},\{k_k,k_l\})&=&\frac{\partial^2}{\partial k_{k,\alpha_1}\partial k_{l, \beta_1}}I(\{1,1\},\{k_k,k_l\}). \end{align}\] The analytic expressions for \(I(\{1,1\},\{k_k,k_l\})\) are known in the literature and are given in appendix A of ref. [56].

In the following, we will present the concrete analytic expressions of the eikonal tensor integrals that we need. In general, it would be convenient to further split the tensor integrals into the pole and finite parts: \[\begin{align} \bar{\mathcal{E}}^{\alpha_1\ldots \alpha_{n_1-1} \beta_1\ldots \beta_{n_2-1}}(\{n_1,n_2\},\{k_k,k_l\})&=&\hat{\mathcal{E}}^{\alpha_1\ldots \alpha_{n_1-1} \beta_1\ldots \beta_{n_2-1}}(\{n_1,n_2\},\{k_k,k_l\}) \nonumber\\ &&+\mathcal{E}^{\alpha_1\ldots \alpha_{n_1-1} \beta_1\ldots \beta_{n_2-1}}(\{n_1,n_2\},\{k_k,k_l\}), \end{align}\] where \(\hat{\mathcal{E}}\) and \(\mathcal{E}\) represent the IR poles and the finite component of \(\bar{\mathcal{E}}\), respectively. We discuss their expressions in four categories.

7.1 Two-massless case↩︎

We first consider the case of two massless external legs (\(m_k=m_l=0\)). Since neither \({\cal I}_k\) nor \({\cal I}_l\) can be a quarkonium, we only need to consider the scalar integral (\(n_1=n_2=1\)), which has been known (e.g.eqs.(A.5-A.6) in ref. [56]). In this case, if \(l=k\), the tensor integrals are zero. If \(l\neq k\), for the completeness, the expressions are \[\begin{align} \hat{\mathcal{E}}(\{1,1\},\{k_k,k_l\})&=&\frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)} \left(\frac{\mu^2}{Q_{\rm ES}^2}\right)^\epsilon\left[ \frac{1}{\epsilon^2}-\frac{1}{\epsilon} \left(\log\!\left(\frac{2 k_k \!\cdot\!k_l}{Q_{\rm ES}^2}\right)- \log\!\left(\frac{4E_kE_l}{\xi_{cut}^2 s}\right)\right)\right]\,,~~~~~ \\ \mathcal{E}(\{1,1\},\{k_k,k_l\})&=& \displaystyle\frac{1}{2}\log^2\!\left(\frac{\xi_{cut}^2 s}{Q_{\rm ES}^2}\right)+ \log\!\left(\frac{\xi_{cut}^2 s}{Q_{\rm ES}^2}\right)\log\!\left(\frac{k_k \!\cdot\!k_l}{2E_kE_l}\right) -{\rm Li}_2\!\left(\frac{k_k \!\cdot\!k_l}{2E_kE_l}\right) \nonumber\\*&&+ \displaystyle\frac{1}{2}\log^2\!\left(\frac{k_k \!\cdot\!k_l}{2E_kE_l}\right) -\log\!\left(1-\frac{k_k \!\cdot\!k_l}{2E_kE_l}\right)\log\!\left(\frac{k_k \!\cdot\!k_l}{2E_kE_l}\right)\,, \end{align}\] where \({\rm Li}_2()\) is the dilogarithm.

7.2 One-massive-one-massless case↩︎

The second case we consider involves one-massive and one-massless external leg (\(m_k=0\) and \(m_l\neq 0\)). The scalar case has been given in eqs.(A.7-A.8) in ref. [56]: \[\begin{align} \hat{\mathcal{E}}(\{1,1\},\{k_k,k_l\})&=&\frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)} \left(\frac{\mu^2}{Q_{\rm ES}^2}\right)^\epsilon\left[ \frac{1}{2\epsilon^2}-\frac{1}{\epsilon} \left(\log\!\left(\frac{2k_k \!\cdot\!k_l}{Q_{\rm ES}^2}\right)- \displaystyle\frac{1}{2}\log\!\left(\frac{4m_l^2E_k^2}{\xi_{cut}^2 s Q_{\rm ES}^2}\right)\right)\right]\,,\nonumber\\ \\ \mathcal{E}(\{1,1\},\{k_k,k_l\})&=& \log\!\left(\xi_{cut}\right)\left(\log\!\left(\frac{\xi_{cut}s}{Q_{\rm ES}^2}\right)+2\log\!\left(\frac{k_k \!\cdot\!k_l}{m_lE_k}\right)\right) -\frac{\pi^2}{12}+\displaystyle\frac{1}{4}\log^2\!\left(\frac{s}{Q_{\rm ES}^2}\right) \nonumber\\*&&- \displaystyle\frac{1}{4}\log^2\!\left(\frac{1+\beta_l}{1-\beta_l}\right) +\displaystyle\frac{1}{2}\log^2\!\left(\frac{k_k \!\cdot\!k_l}{(1-\beta_l)E_k E_l}\right) +\log\!\left(\frac{s}{Q_{\rm ES}^2}\right)\log\!\left(\frac{k_k \!\cdot\!k_l}{m_lE_k}\right) \nonumber\\*&&- {\rm Li}_2\!\left(1-\frac{(1+\beta_l)E_kE_l}{k_k \!\cdot\!k_l}\right) +{\rm Li}_2\!\left(1-\frac{k_k \!\cdot\!k_l}{(1-\beta_l)E_kE_l}\right)\,, \end{align}\] where \[\beta_l=\sqrt{1-\frac{m_l^2}{E_l^2}}.\] The IR pole part of the rank-1 tensor integral needed in the P-wave soft integrated counterterms is \[\begin{align} &&\hat{\mathcal{E}}^{\mu}(\{1,2\},\{k_k,k_l\}) \nonumber\\ &=&\displaystyle\frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)}\left(\displaystyle\frac{\mu^2}{Q_{\rm ES}^2}\right)^\epsilon\!\left\lbrace \displaystyle\frac{1}{2\epsilon^2} \displaystyle\frac{k_k^{\mu }}{k_k \!\cdot\!k_l} + \displaystyle\frac{1}{\epsilon} \left[\displaystyle\frac{k_k^{\mu }}{k_k \!\cdot\!k_l}\! \left(1-\log\!\left(\displaystyle\frac{k_k \!\cdot\!k_l}{E_k m_l}\right)-\displaystyle\frac{1}{2}\log\!\left(\displaystyle\frac{\xi_{cut}^2 s}{Q_{\rm ES}^2}\right)\right)-\displaystyle\frac{k_l^{\mu}}{m_l^2} \right]\right\rbrace\,,\nonumber\\ \end{align}\] while its finite part can be decomposed into 3 tensor structures \[\mathcal{E}^{\mu}(\{1,2\},\{k_k,k_l\})= \sum_{i=1}^3{T_{i,12}^{(0,m_l),\mu}},\] where \[\begin{align} T_{1,12}^{(0,m_l),\mu}&=& \displaystyle\frac{k_k^{\mu}}{k_k \!\cdot\!k_l}\left\lbrace-\displaystyle\frac{\pi ^2}{12}-\log\! \left(\displaystyle\frac{\xi_{cut}^2 s}{Q_{\rm ES}^2}\right)+\displaystyle\frac{1}{4} \log^2\!\left(\displaystyle\frac{\xi_{cut}^2 s}{Q_{\rm ES}^2}\right)+\log\! \left(\displaystyle\frac{k_k \!\cdot\!k_l}{E_k m_l}\right)\log\! \left(\displaystyle\frac{\xi_{cut}^2 s}{Q_{\rm ES}^2}\right)\right.\nonumber\\ &&-\log\! \left(\displaystyle\frac{k_k \!\cdot\!k_l}{E_k E_l \left(1-\beta_l\right)}\right)+\displaystyle\frac{1}{2} \log\! ^2\left(\displaystyle\frac{k_k \!\cdot\!k_l}{E_k E_l \left(1-\beta_l\right)}\right)-\displaystyle\frac{1}{4} \log^2\!\left(\displaystyle\frac{1+\beta_l}{1-\beta_l}\right)\nonumber\\ &&+\mathrm{Li}_2\!\left(1-\displaystyle\frac{k_k \!\cdot\!k_l}{E_k E_l \left(1-\beta _l\right)}\right)-\mathrm{Li}_2\!\left(1-\displaystyle\frac{E_k E_l \left(1+\beta_l\right)}{k_k \!\cdot\!k_l}\right)\nonumber\\ &&+\displaystyle\frac{k_k \!\cdot\!k_l}{k_k \!\cdot\!k_l-E_k E_l(1-\beta_l)}\log\! \left(\displaystyle\frac{k_k \!\cdot\!k_l}{E_k E_l \left(1-\beta_l\right)}\right)+ \displaystyle\frac{E_k E_l(1+\beta_l)}{k_k \!\cdot\!k_l-E_k E_l(1+\beta_l)}\nonumber\\ &&\left.\times\log\! \left(\displaystyle\frac{k_k \!\cdot\!k_l}{E_k E_l \left(1+\beta_l\right)}\right)\right\rbrace,\nonumber \end{align}\] \[\begin{align} T_{2,12}^{(0,m_l),\mu} &=& \displaystyle\frac{k_l^{\mu}}{m_l^2}\left\lbrace \log\!\left(\displaystyle\frac{\xi_{cut}^2 s}{Q_{\rm ES}^2}\right) + \displaystyle\frac{m_l^2}{\beta_l E_l^2}\left[-\displaystyle\frac{1}{1-\beta_l^2}\log\!\left(\displaystyle\frac{1+\beta_l}{1-\beta_l}\right) -\displaystyle\frac{E_k E_l}{k_k \!\cdot\!k_l - E_k E_l (1-\beta_l)} \right.\right.\nonumber\\ &&\left.\left.\times\log\!\left(\displaystyle\frac{k_k \!\cdot\!k_l}{E_k E_l(1-\beta_l)}\right)+\displaystyle\frac{E_kE_l}{k_k \!\cdot\!k_l - E_k E_l (1+\beta_l)}\log\!\left(\displaystyle\frac{k_k \!\cdot\!k_l}{E_k E_l(1+\beta_l)}\right)\right] \right\rbrace,\nonumber\\ T_{3,12}^{(0,m_l),\mu}&=& \displaystyle\frac{\delta^{\mu0}}{k_k \!\cdot\!k_l-E_k E_l \left(1-\beta_l\right)}\left\lbrace E_k \left\lbrack\log\! \left(\displaystyle\frac{E_k E_l \left(1+\beta _l\right)}{k_k \!\cdot\!k_l}\right)-\log\! \left(\displaystyle\frac{k_k \!\cdot\!k_l}{E_k E_l \left(1-\beta _l\right)}\right)\right\rbrack\right.\nonumber\\ &&+\displaystyle\frac{k_k \!\cdot\!k_l \, m_l^2}{E_l^3 \beta _l \left(1-\beta _l^2\right)}\log\! \left(\displaystyle\frac{1+\beta _l}{1-\beta _l}\right) +\displaystyle\frac{E_k \left(E_l^2 \beta _l \left(1+\beta _l\right)+m_l^2\right)}{E_l^2 \left(k_k \!\cdot\!k_l-E_k E_l \left(1+\beta _l\right)\right)}\nonumber\\ &&\left.\times\left\lbrack E_k E_l\! \left(\log\! \left(\displaystyle\frac{E_k E_l \left(1+\beta _l\right)}{k_k \!\cdot\!k_l}\right)\!-\!\log\! \left(\displaystyle\frac{k_k \!\cdot\!k_l}{E_k E_l \left(1-\beta _l\right)}\right)\right)+\displaystyle\frac{k_k \!\cdot\!k_l}{1+\beta _l}\log\! \left(\displaystyle\frac{1+\beta _l}{1-\beta _l}\right)\right\rbrack\right\rbrace.\nonumber\\ \end{align}\]

7.3 Massive self-eikonal case↩︎

In the case of the massive self-eikonal integrals, we have \(l=k\) and \(m_k=m_l\neq 0\). The scalar integral corresponds to eqs.(A.9-A.10) in ref. [56]: \[\begin{align} \hat{\mathcal{E}}(\{1,1\},\{k_k,k_k\})&=&\frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)} \left(\frac{\mu^2}{Q_{\rm ES}^2}\right)^\epsilon\left(-\frac{1}{\epsilon}\right)\,,\\ \mathcal{E}(\{1,1\},\{k_k,k_k\})&=& \log\!\left(\frac{\xi_{cut}^2 s}{Q_{\rm ES}^2}\right) -\frac{1}{\beta_k}\log\!\left(\frac{1+\beta_k}{1-\beta_k}\right)\,. \end{align}\] For the rank-1 tensor integral, the poles are \[\hat{\mathcal{E}}^{\mu}(\{1,2\},\{k_k,k_k\})= \displaystyle\frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)}\left(\displaystyle\frac{\mu^2}{Q_{\rm ES}^2}\right)^\epsilon\displaystyle\frac{k_k^{\mu }}{m_k^2} \left(-\displaystyle\frac{1}{\epsilon} \right),\] and the finite term is \[\begin{align} \mathcal{E}^{\mu}(\{1,2\},\{k_k,k_k\})&=&\displaystyle\frac{k_k^{\mu }}{2\beta_k^3 m_k^2} \left\lbrack \log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right) -2 \beta_k -3\log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right)\beta_k^2 + 2 \log\!\left(\displaystyle\frac{\xi_{cut}^2 s}{Q_{\rm ES}^2}\right)\beta_k^3 \right\rbrack\nonumber\\ && + \displaystyle\frac{\delta^{\mu0}}{2\beta_k^3 E_k}\left\lbrack-\log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right)+2 \beta_k+\log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right)\beta_k^2\right\rbrack. \end{align}\] The pole of the rank-2 integral is \[\hat{\mathcal{E}}^{\mu\nu}(\{2,2\},\{k_k,k_k\})= \displaystyle\frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)}\left(\displaystyle\frac{\mu^2}{Q_{\rm ES}^2}\right)^\epsilon\displaystyle\frac{1}{3 m_k^2} \left\lbrack -\displaystyle\frac{1}{\epsilon} \left(-g^{\mu\nu}+\displaystyle\frac{4 k_k^{\mu}k_k^{\nu}}{m_k^2}\right)\right\rbrack.\] The finite piece is \[\mathcal{E}^{\mu\nu}(\{2,2\},\{k_k,k_k\})= \sum_{i=1}^4{T^{(m_k,m_k),\mu\nu}_{i,22}},\] where \[\begin{align} T_{1,22}^{(m_k,m_k),\mu\nu}&=&\displaystyle\frac{\delta^{\mu0}\delta^{\nu0}}{6\beta_k^5 m_k^2}(1-\beta_k^2)\left\lbrack - 3\log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right) + 6\beta_k + 3 \log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right)\beta_k^2 - 4 \beta_k^3 \right\rbrack,\nonumber\\ T_{2,22}^{(m_k,m_k),\mu\nu} &=& \displaystyle\frac{g^{\mu\nu}}{6\beta_k^3 m_k^2}\left\lbrack -\log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right) + 2\beta_k + 3\log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right)\beta_k^2 - 2\log\!\left(\displaystyle\frac{\xi_{cut}^2 s}{Q^2_{\rm ES}}\right)\beta_k^3 \right\rbrack,\nonumber\\ T_{3,22}^{(m_k,m_k),\mu\nu}&=&\displaystyle\frac{\delta^{\mu0}k_k^{\nu}+k_k^{\mu} \delta^{\nu0}}{6\beta_k^5 E_k m_k^2}\left\lbrack 3\log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right) - 6\beta_k - 6\log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right)\beta_k^2 \right.\nonumber\\ &&\left.+ 10\beta_k^3 + 3\log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right)\beta_k^4 \right\rbrack,\nonumber\\ T_{4,22}^{(m_k,m_k),\mu\nu}&=&\displaystyle\frac{k_k^{\mu}k_k^{\nu}}{6\beta_k^5 m_k^4}\left\lbrack -3\log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right) + 6\beta_k + 10\log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right)\beta_k^2 \right.\nonumber\\ &&\left.- 18\beta_k^3 - 15\log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right)\beta_k^4 + 8\log\!\left(\displaystyle\frac{\xi_{cut}^2 s}{Q_{\rm ES}^2}\right)\beta_k^5 \right\rbrack. \end{align}\]

7.4 Two-massive case↩︎

Finally, let us consider the most complicated case with \(l\neq k\), \(m_k\neq0\) and \(m_l\neq 0\). The scalar integral case can be referred to eqs.(A.11-A.12) in ref. [56]. Its expression is \[\begin{align} \hat{\mathcal{E}}(\{1,1\},\{k_k,k_l\})&=&\frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)} \left(\frac{\mu^2}{Q_{\rm ES}^2}\right)^\epsilon\left( -\frac{1}{2\epsilon}\frac{1}{v_{kl}}\log\!\left(\frac{1+v_{kl}}{1-v_{kl}}\right)\right)\,, \\ \mathcal{E}(\{1,1\},\{k_k,k_l\})&=& \frac{1}{2v_{kl}}\log\!\left(\frac{1+v_{kl}}{1-v_{kl}}\right) \log\!\left(\frac{\xi_{cut}^2 s}{Q_{\rm ES}^2}\right) \nonumber\\*&&+ \frac{(1+v_{kl})(k_k \!\cdot\!k_l)^2}{2m_k^2} \left({\rm J}^{(A)}\left(\alpha_{kl} E_k,\alpha_{kl} E_k\beta_k\right) -{\rm J}^{(A)}\left(E_l,E_l\beta_l\right)\right)\,,\nonumber\\ \end{align}\] where the introduced function is \[\begin{align} {\rm J}^{(A)}\left(x,y\right)&=&\frac{1}{2\lambda\nu}\left\lbrack\log^2\!\left(\frac{x-y}{x+y}\right)+4\mathrm{Li}_2\left(1-\frac{x+y}{\nu}\right)+4\mathrm{Li}_2\left(1-\frac{x-y}{\nu}\right)\right\rbrack \end{align}\] and \[\begin{align} v_{kl}&=&\sqrt{1-\left(\displaystyle\frac{m_km_l}{k_k\!\cdot\!k_l}\right)^2},\quad \alpha_{kl}=\displaystyle\frac{1+v_{kl}}{m_k^2}k_k\!\cdot\!k_l,\nonumber\\ \lambda&=&\alpha_{kl}E_k-E_l,\quad \nu=\displaystyle\frac{\alpha_{kl}^2m_k^2-m_l^2}{2\lambda}. \end{align}\]

For the two-massive rank-1 eikonal tensor integral, its pole is \[\begin{align} \hat{\mathcal{E}}^{\mu}(\{1,2\},\{k_k,k_l\}) &=& \displaystyle\frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)}\left(\displaystyle\frac{\mu^2}{Q_{\rm ES}^2}\right)^\epsilon\displaystyle\frac{1}{2 \epsilon}\displaystyle\frac{m_k^2}{(k_k \!\cdot\!k_l)^2 v_{kl}^3} \left\lbrace -k_k^{\mu}\left\lbrack \displaystyle\frac{k_k \!\cdot\!k_l}{m_k^2} \log\!\left(\displaystyle\frac{1+v_{kl}}{1-v_{kl}}\right)v_{kl}^2 \right.\right.\nonumber\\ &&\left.\left. + \displaystyle\frac{m_l^2}{k_k \!\cdot\!k_l}\left(\log\!\left(\displaystyle\frac{1+v_{kl}}{1-v_{kl}}\right) - \displaystyle\frac{2v_{kl}}{1-v_{kl}^2} \right)\right\rbrack + k_l^{\mu} \left(\log\!\left(\displaystyle\frac{1+v_{kl}}{1-v_{kl}}\right) - \displaystyle\frac{2v_{kl}}{1-v_{kl}^2} \right) \right\rbrace.\nonumber\\ \end{align}\] The finite term can be decomposed into \(8\) tensorial components \[\mathcal{E}^{\mu}(\{1,2\},\{k_k,k_l\})= \sum_{i=1}^8{T_{i,12}^{(m_k,m_l),\mu}},\] where \[\begin{align} T_{1,12}^{(m_k,m_l),\mu} &=& k_k^{\mu} \left\lbrack \displaystyle\frac{1}{2(k_k \!\cdot\!k_l)v_{kl}}\log\!\left(\displaystyle\frac{\xi_{cut}^2 s}{Q_{\rm ES}^2}\right) \log\!\left(\displaystyle\frac{1+v_{kl}}{1-v_{kl}}\right) \right.\nonumber\\ &&\left. - \displaystyle\frac{(1+v_{kl})(k_k \!\cdot\!k_l)}{2 m_k^2} \left( \mathrm{J}^{(A)}(\alpha_{kl}E_k,\alpha_{kl} \beta_k E_k) - \mathrm{J}^{(A)}(E_l,\beta_l E_l) \right) \right\rbrack,\nonumber\\ T_{2,12}^{(m_k,m_l),\mu} &=& \displaystyle\frac{1}{2v_{kl}}\left(k_k^{\mu}\displaystyle\frac{m^2_l}{(k_k \!\cdot\!k_l)} - k_l^{\mu}\right) \left\lbrace \displaystyle\frac{m_k^2}{(k_k \!\cdot\!k_l)^2 v_{kl}^2} \log\!\left(\displaystyle\frac{\xi_{cut}^2 s}{Q_{\rm ES}^2}\right) \left\lbrack \log\!\left(\displaystyle\frac{1+v_{kl}}{1-v_{kl}}\right) - \displaystyle\frac{2v_{kl}}{1-v_{kl}^2} \right\rbrack \right.\nonumber\\ &&\left. - \left( \mathrm{J}^{(A)}(\alpha_{kl}E_k,\alpha_{kl} \beta_k E_k) - \mathrm{J}^{(A)}(E_l,\beta_l E_l) \right) \right\rbrace,\nonumber\\ T_{3,12}^{(m_k,m_l),\mu} &=& \displaystyle\frac{(1+v_{kl})E_k (k_k \!\cdot\!k_l)}{\lambda\nu \left\lbrack(\nu-\alpha_{kl}E_k)^2-(\alpha_{kl}\beta_k E_k)^2\right\rbrack v_{kl} m_k^2}\left(-k_k^{\mu}\left(\alpha_{kl}v_{kl}+\displaystyle\frac{m_l^2}{k_k \!\cdot\!k_l}\right) + k_l^{\mu}\right)\nonumber\\ &&\left\lbrace \left(\nu-\alpha_{kl}E_k(1-\beta_k^2)\right)\left\lbrack \log\!\left(\displaystyle\frac{\alpha_{kl}E_k(1-\beta_{k})}{\nu}\right) + \log\!\left(\displaystyle\frac{\alpha_{kl}E_k(1+\beta_{k})}{\nu}\right) \right\rbrack \right.\nonumber\\ &&\left.+ \nu \beta_k \log\!\left(\displaystyle\frac{1+\beta_k}{1-\beta_k}\right) \right\rbrace,\nonumber\\ T_{4,12}^{(m_k,m_l),\mu}&=&\displaystyle\frac{1+v_{kl}}{2 \lambda^2 \nu^2 v_{kl}m_k^2}\left\lbrack-k_k^{\mu}\left(\nu E_k - \alpha_{kl}m_k^2\right)\left(\alpha_{kl}v_{kl}(k_k \!\cdot\!k_l)+m_l^2\right)\right.\nonumber\\ && \left.+k_l^{\mu}(k_k \!\cdot\!k_l)(\nu E_k-v_{kl}(k_k \!\cdot\!k_l)-\alpha_{kl} m_k^2) + \delta^{\mu0} \nu v_{kl} (k_k \!\cdot\!k_l)^2 \right\rbrack \nonumber\\ && \times \left\lbrack \displaystyle\frac{2\alpha_{kl}E_k(1-\beta_k)}{\nu-\alpha_{kl}E_k(1-\beta_k)}\log\!\left(\displaystyle\frac{\alpha_{kl}E_k(1-\beta_k)}{\nu}\right) + \displaystyle\frac{2\alpha_{kl}E_k(1+\beta_k)}{\nu-\alpha_{kl}E_k(1+\beta_k)} \right.\nonumber\\ &&\left. \times \log\!\left(\displaystyle\frac{\alpha_{kl}E_k(1+\beta_k)}{\nu}\right)+ \lambda \nu \mathrm{J}^{(A)}(\alpha_{kl}E_k,\alpha_{kl}\beta_k E_k) \right\rbrack,\nonumber\\ T_{5,12}^{(m_k,m_l),\mu} &=&\displaystyle\frac{1+v_{kl}}{2\lambda v_{kl} m_k^2} \left\lbrack k_k^{\mu}E_k(\alpha_{kl}v_{kl}(k_k \!\cdot\!k_l)+m_l^2) - k_l^{\mu}E_k(k_k \!\cdot\!k_l) - \delta^{\mu0} v_{kl} (k_k \!\cdot\!k_l)^2 \right\rbrack \nonumber\\ &&\times\left\lbrack \mathrm{J}^{(A)}(\alpha_{kl}E_k,\alpha_{kl} \beta_k E_k) - \mathrm{J}^{(A)}(E_l,\beta_l E_l) \right\rbrack,\nonumber\\ T_{6,12}^{(m_k,m_l),\mu} &=& \displaystyle\frac{(1+v_{kl})(k_k \!\cdot\!k_l)^2}{\lambda \nu m_k^2\left\lbrack (\nu-E_l)^2-(\beta_l E_l)^2\right\rbrack} \delta^{\mu0} \left\lbrace \left( \nu - E_l(1-\beta_l^2) \right) \left\lbrack \log\!\left(\displaystyle\frac{E_l(1-\beta_l)}{\nu}\right)\right. \right.\nonumber\\ &&\left.\left.+ \log\!\left(\displaystyle\frac{E_l(1+\beta_l)}{\nu}\right) \right\rbrack + \nu \beta_l \log\!\left(\displaystyle\frac{1+\beta_l}{1-\beta_l}\right) \right\rbrace,\nonumber \end{align}\] \[\begin{align} T_{7,12}^{(m_k,m_l),\mu}&=&\displaystyle\frac{(1+v_{kl})(k_k \!\cdot\!k_l)^2}{\lambda \nu \beta_l E_l^2 m_k^2} \bigg(k_l^{\mu} E_l - \delta^{\mu0} m_l^2\bigg) \left\lbrace \displaystyle\frac{1}{\nu - E_l(1-\beta_l)}\log\!\left(\displaystyle\frac{E_l(1-\beta_l)}{\nu}\right) \right.\nonumber\\ &&\left. - \displaystyle\frac{1}{\nu - E_l(1+\beta_l)}\log\!\left(\displaystyle\frac{E_l(1+\beta_l)}{\nu}\right) - \displaystyle\frac{1}{E_l(1-\beta_l^2)}\log\!\left(\displaystyle\frac{1+\beta_l}{1-\beta_l}\right) \right\rbrace,\nonumber\\ T_{8,12}^{(m_k,m_l),\mu}&=& \displaystyle\frac{1+v_{kl}}{2 \lambda^2 \nu^2 v_{kl}m_k^2}\left\lbrack k_k^{\mu}\left(\nu E_k - \alpha_{kl}m_k^2\right)\left(\alpha_{kl}v_{kl}(k_k \!\cdot\!k_l)+m_l^2\right) \right.\nonumber\\ && \left.- k_l^{\mu}(k_k \!\cdot\!k_l)(\nu E_k-v_{kl}(k_k \!\cdot\!k_l)-\alpha_{kl} m_k^2) - \delta^{\mu0} \nu v_{kl} (k_k \!\cdot\!k_l)^2 \right\rbrack \nonumber\\ &&\times \left\lbrack \displaystyle\frac{2 E_l(1-\beta_l)}{\nu-E_l(1-\beta_l)}\log\!\left(\displaystyle\frac{E_l(1-\beta_l)}{\nu}\right) + \displaystyle\frac{2 E_l(1+\beta_l)}{\nu-E_l(1+\beta_l)}\log\!\left(\displaystyle\frac{E_l(1+\beta_l)}{\nu}\right) \right.\nonumber\\ &&\left. + \lambda \nu \mathrm{J}^{(A)}(E_l,\beta_l E_l) \right\rbrack. \end{align}\] We have verified that if we take \(k=l\), the two-massive eikonal integrals are reduced to the massive self-eikonal case.

8 Infrared poles of one-loop matrix elements↩︎

We present here the general IR poles of the UV renormalised one-loop matrix element at NLO QCD for arbitrary process that involves a quarkonium and elementary particles by assuming the validity of the KLN theorem. The one-loop virtual matrix element can be written as \[\begin{align} {\mathbb{M}}^{(n-1,1)}(\dot{r}^{i\backslash})=\frac{\alpha_s}{2\pi}\frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)} \left(\frac{\mu^2}{Q_{\rm ES}^2}\right)^\epsilon{\mathbb{V}}(\dot{r}^{i\backslash})\,. \end{align}\] We can derive the IR poles of \({\mathbb{V}}(\dot{r}^{i\backslash})\) generally from our integrated counterterms with its finite term denoted as \({\mathbb{V}}^{(n-1,1)}_{\rm FIN}(\dot{r}^{i\backslash})\). We remind readers that in the process \(\dot{r}^{i\backslash}\), the particles with their indices from \(n_{I}\) to \(n_{L}^{(B)}+2\) are coloured and massless, while those from \(n_{L}^{(B)}+3\) to \(n_{L}^{(B)}+n_{H}\) (\(n_{H}\geq 2\)) are massive coloured elementary particles. The index of the quarkonium is \(n_{L}^{(B)}+n_{H}+1\).

8.1 Colour-singlet S-wave state↩︎

For a colour-singlet S-wave state with \(\dot{r}^{i\backslash}\) being eq.@eq:eq:BornCSSwave , the IR poles can be obtained by taking the bound state as a colour-singlet elementary particle, which amount to \[\begin{align} {\mathbb{V}}(\dot{r}^{i\backslash})&=&-\Bigg( \frac{1}{\epsilon^2}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}C({\cal I}_k) +\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}\gamma({\cal I}_k) +\frac{1}{\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}}C({\cal I}_k) \Bigg){\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash}) \nonumber\\ &&+\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2} \sum_{l=k+1}^{n_{L}^{(B)}+n_{H}}\log\!\left(\frac{2k_k\!\cdot\!k_l}{Q_{\rm ES}^2}\right) {\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\ &&+\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}-1} \sum_{l=k+1}^{n_{L}^{(B)}+n_{H}} \frac{1}{v_{kl}}\log\!\left(\frac{1+v_{kl}}{1-v_{kl}}\right) {\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\ &&-\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}} \log\!\left(\frac{m_k^2}{Q_{\rm ES}^2}\right) \sum_{l=n_{I}}^{n_{L}^{(B)}+2}{\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) +{\mathbb{V}}^{(n-1,1)}_{\rm FIN}(\dot{r}^{i\backslash})\,. \end{align}\] This equation is equivalent to eq.(B.2) in ref. [56].

8.2 Colour-octet S-wave state↩︎

The colour-octet S-wave case with \(\dot{r}^{i\backslash}\) being eq.@eq:eq:BornCOSwave has the IR poles as follows: \[\begin{align} {\mathbb{V}}(\dot{r}^{i\backslash})&=&-\Bigg( \frac{1}{\epsilon^2}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}C({\cal I}_k) +\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}\gamma({\cal I}_k) +\frac{1}{\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}+1}C({\cal I}_k) \Bigg){\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash}) \nonumber\\ &&+\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2} \sum_{l=k+1}^{n_{L}^{(B)}+n_{H}+1}\log\!\left(\frac{2k_k\!\cdot\!k_l}{Q_{\rm ES}^2}\right) {\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\ &&+\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}} \sum_{l=k+1}^{n_{L}^{(B)}+n_{H}+1} \frac{1}{v_{kl}}\log\!\left(\frac{1+v_{kl}}{1-v_{kl}}\right) {\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\ &&-\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}+1} \log\!\left(\frac{m_k^2}{Q_{\rm ES}^2}\right) \sum_{l=n_{I}}^{n_{L}^{(B)}+2}{\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) +{\mathbb{V}}^{(n-1,1)}_{\rm FIN}(\dot{r}^{i\backslash})\,. \end{align}\]

8.3 Colour-singlet spin-singlet P-wave state↩︎

The IR poles of the virtual matrix elements for a colour-singlet spin-singlet P-wave quarkonium production are \[\begin{align} {\mathbb{V}}(\dot{r}^{i\backslash})&=&-\Bigg( \frac{1}{\epsilon^2}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}C({\cal I}_k) +\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}\gamma({\cal I}_k) +\frac{1}{\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}}C({\cal I}_k) \Bigg){\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash}) \nonumber\\ &&+\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2} \sum_{l=k+1}^{n_{L}^{(B)}+n_{H}}\log\!\left(\frac{2k_k\!\cdot\!k_l}{Q_{\rm ES}^2}\right) {\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\ &&+\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}-1} \sum_{l=k+1}^{n_{L}^{(B)}+n_{H}} \frac{1}{v_{kl}}\log\!\left(\frac{1+v_{kl}}{1-v_{kl}}\right) {\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\ &&-\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}} \log\!\left(\frac{m_k^2}{Q_{\rm ES}^2}\right) \sum_{l=n_{I}}^{n_{L}^{(B)}+2}{\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\ &&+\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}\displaystyle\frac{k_{k,\mu}}{K\!\cdot\!k_k}{\mathbb{M}}^{(n-1,0)}_{k[18]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu\nonumber\\ &&+\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}}\displaystyle\frac{k_{k,\mu}}{K\!\cdot\!k_k} \frac{m_k^2K^2}{ v_{k}^3(K\!\cdot\!k_k)^2}\left[\frac{2v_k}{1-v_k^2}-\log\!\left(\frac{1+v_k}{1-v_k}\right)\right]{\mathbb{M}}^{(n-1,0)}_{k[18]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu\nonumber\\ &&-\frac{1}{\epsilon}\left(\frac{1}{2N_{c}}\frac{8}{m_{Q}m_{\bar{Q}^\prime}}-\frac{2}{K^2}C_{{\rm eff}}(Q\bar{Q}^\prime_{[18]})\right){\mathbb{M}}^{(n-1,0)}(\dot{r}_1^{i\backslash})^{\mu\nu}+{\mathbb{V}}^{(n-1,1)}_{\rm FIN}(\dot{r}^{i\backslash})\,, \end{align}\] where we have defined the abbreviation for \(v_{kl}\) with \(l=n_{L}^{(B)}+n_{H}+1\) (quarkonium) \[\begin{align} v_k&=&\sqrt{1-\frac{k_k^2K^2}{(k_k\cdot K)^2}}. \end{align}\] Note that, the two reduced Born processes \(\dot{r}^{i\backslash}\) and \(\dot{r}_1^{i\backslash}\) have been given in eq.@eq:eq:BornCS1Pwave .

8.4 Colour-octet spin-singlet P-wave state↩︎

For the colour-octet spin-singlet P-wave state, the IR poles of the virtual corrections are \[\begin{align} {\mathbb{V}}(\dot{r}^{i\backslash})&=&-\Bigg( \frac{1}{\epsilon^2}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}C({\cal I}_k) +\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}\gamma({\cal I}_k) +\frac{1}{\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}+1}C({\cal I}_k) \Bigg){\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash}) \nonumber\\*&& +\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2} \sum_{l=k+1}^{n_{L}^{(B)}+n_{H}+1}\log\!\left(\frac{2k_k\!\cdot\!k_l}{Q_{\rm ES}^2}\right) {\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\*&& +\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}} \sum_{l=k+1}^{n_{L}^{(B)}+n_{H}+1} \frac{1}{v_{kl}}\log\!\left(\frac{1+v_{kl}}{1-v_{kl}}\right) {\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\*&& -\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}+1} \log\!\left(\frac{m_k^2}{Q_{\rm ES}^2}\right) \sum_{l=n_{I}}^{n_{L}^{(B)}+2}{\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\*&& +\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}\displaystyle\frac{k_{k,\mu}}{K\!\cdot\!k_k}\left({\mathbb{M}}^{(n-1,0)}_{k[88]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu+{\mathbb{M}}^{(n-1,0)}_{k[81]}(\dot{r}^{i\backslash},\dot{r}_2^{i\backslash})^\mu\right)\nonumber\\ &&+\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}}\displaystyle\frac{k_{k,\mu}}{K\!\cdot\!k_k} \frac{m_k^2K^2}{ v_{k}^3(K\!\cdot\!k_k)^2}\left[\frac{2v_k}{1-v_k^2}-\log\!\left(\frac{1+v_k}{1-v_k}\right)\right]\nonumber\\ &&\times\left({\mathbb{M}}^{(n-1,0)}_{k[88]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu+{\mathbb{M}}^{(n-1,0)}_{k[81]}(\dot{r}^{i\backslash},\dot{r}_2^{i\backslash})^\mu\right)\nonumber\\ &&-\frac{1}{\epsilon}\left(B_{F}\frac{8}{m_{Q}m_{\bar{Q}^\prime}}-\frac{2}{K^2}C_{{\rm eff}}(Q\bar{Q}^\prime_{[88]})\right){\mathbb{M}}^{(n-1,0)}(\dot{r}_1^{i\backslash})\nonumber\\ &&-\frac{1}{\epsilon}\left(C_{F}\frac{8}{m_{Q}m_{\bar{Q}^\prime}}-\frac{2}{K^2}C_{{\rm eff}}(Q\bar{Q}^\prime_{[81]})\right){\mathbb{M}}^{(n-1,0)}(\dot{r}_2^{i\backslash})+{\mathbb{V}}^{(n-1,1)}_{\rm FIN}(\dot{r}^{i\backslash})\,. \end{align}\] The three reduced Born processes \(\dot{r}^{i\backslash}\), \(\dot{r}_1^{i\backslash}\) and \(\dot{r}_2^{i\backslash}\) have been given in eq.@eq:eq:BornCO1Pwave .

8.5 Colour-singlet spin-triplet P-wave state↩︎

The IR poles of the one-loop virtual corrections for the processes involving a colour-singlet spin-triple P-wave quarkonium are \[\begin{align} {\mathbb{V}}(\dot{r}^{i\backslash})&=&-\Bigg( \frac{1}{\epsilon^2}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}C({\cal I}_k) +\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}\gamma({\cal I}_k) +\frac{1}{\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}}C({\cal I}_k) \Bigg){\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash}) \nonumber\\ &&+\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2} \sum_{l=k+1}^{n_{L}^{(B)}+n_{H}}\log\!\left(\frac{2k_k\!\cdot\!k_l}{Q_{\rm ES}^2}\right) {\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\ &&+\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}-1} \sum_{l=k+1}^{n_{L}^{(B)}+n_{H}} \frac{1}{v_{kl}}\log\!\left(\frac{1+v_{kl}}{1-v_{kl}}\right) {\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\ &&-\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}} \log\!\left(\frac{m_k^2}{Q_{\rm ES}^2}\right) \sum_{l=n_{I}}^{n_{L}^{(B)}+2}{\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\ &&+\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}\displaystyle\frac{k_{k,\mu}}{K\!\cdot\!k_k}{\mathbb{M}}^{(n-1,0)}_{J,k[18]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu\nonumber\\ &&+\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}}\displaystyle\frac{k_{k,\mu}}{K\!\cdot\!k_k} \frac{m_k^2K^2}{ v_{k}^3(K\!\cdot\!k_k)^2}\left[\frac{2v_k}{1-v_k^2}-\log\!\left(\frac{1+v_k}{1-v_k}\right)\right]{\mathbb{M}}^{(n-1,0)}_{J,k[18]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu\nonumber\\ &&-\frac{1}{\epsilon}\frac{2}{3K^2}g_{\mu\nu}C_{{\rm eff}}(Q\bar{Q}^\prime_{[18]}){\mathbb{M}}^{(n-1,0)}_J(\dot{r}_1^{i\backslash})^{\mu\nu}-\frac{1}{\epsilon}\frac{2J+1}{2N_{c}}\frac{8}{9m_{Q}m_{\bar{Q}^\prime}}{\mathbb{M}}^{(n-1,0)}(\dot{r}_1^{i\backslash})\nonumber\\ &&+{\mathbb{V}}^{(n-1,1)}_{\rm FIN}(\dot{r}^{i\backslash})\,, \end{align}\] where the processes \(\dot{r}^{i\backslash}\) and \(\dot{r}_1^{i\backslash}\) are eq.@eq:eq:BornCS3Pwave .

8.6 Colour-octet spin-triplet P-wave state↩︎

Similarly, the IR poles of the virtual corrections for the colour-octet spin-triplet P-wave quarkonium are \[\begin{align} {\mathbb{V}}(\dot{r}^{i\backslash})&=&-\Bigg( \frac{1}{\epsilon^2}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}C({\cal I}_k) +\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}\gamma({\cal I}_k) +\frac{1}{\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}+1}C({\cal I}_k) \Bigg){\mathbb{M}}^{(n-1,0)}(\dot{r}^{i\backslash}) \nonumber\\*&& +\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2} \sum_{l=k+1}^{n_{L}^{(B)}+n_{H}+1}\log\!\left(\frac{2k_k\!\cdot\!k_l}{Q_{\rm ES}^2}\right) {\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\*&& +\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}} \sum_{l=k+1}^{n_{L}^{(B)}+n_{H}+1} \frac{1}{v_{kl}}\log\!\left(\frac{1+v_{kl}}{1-v_{kl}}\right) {\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\*&& -\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}+1} \log\!\left(\frac{m_k^2}{Q_{\rm ES}^2}\right) \sum_{l=n_{I}}^{n_{L}^{(B)}+2}{\mathbb{M}}^{(n-1,0)}_{kl}(\dot{r}^{i\backslash}) \nonumber\\*&& +\frac{1}{\epsilon}\sum_{k=n_{I}}^{n_{L}^{(B)}+2}\displaystyle\frac{k_{k,\mu}}{K\!\cdot\!k_k}\left({\mathbb{M}}^{(n-1,0)}_{J,k[88]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu+{\mathbb{M}}^{(n-1,0)}_{J,k[81]}(\dot{r}^{i\backslash},\dot{r}_2^{i\backslash})^\mu\right)\nonumber\\ &&+\frac{1}{2\epsilon}\sum_{k=n_{L}^{(B)}+3}^{n_{L}^{(B)}+n_{H}}\displaystyle\frac{k_{k,\mu}}{K\!\cdot\!k_k} \frac{m_k^2K^2}{ v_{k}^3(K\!\cdot\!k_k)^2}\left[\frac{2v_k}{1-v_k^2}-\log\!\left(\frac{1+v_k}{1-v_k}\right)\right]\nonumber\\ &&\times\left({\mathbb{M}}^{(n-1,0)}_{J,k[88]}(\dot{r}^{i\backslash},\dot{r}_1^{i\backslash})^\mu+{\mathbb{M}}^{(n-1,0)}_{J,k[81]}(\dot{r}^{i\backslash},\dot{r}_2^{i\backslash})^\mu\right)\nonumber\\ &&-\frac{1}{\epsilon}\frac{2}{3K^2}g_{\mu\nu}\left(C_{{\rm eff}}(Q\bar{Q}^\prime_{[88]}){\mathbb{M}}^{(n-1,0)}_J(\dot{r}_1^{i\backslash})^{\mu\nu}+C_{{\rm eff}}(Q\bar{Q}^\prime_{[81]}){\mathbb{M}}^{(n-1,0)}_J(\dot{r}_2^{i\backslash})^{\mu\nu}\right)\nonumber\\ &&-\frac{1}{\epsilon}\frac{8(2J+1)}{9m_{Q}m_{\bar{Q}^\prime}}\left(B_{F}{\mathbb{M}}^{(n-1,0)}(\dot{r}_1^{i\backslash})+C_{F}{\mathbb{M}}^{(n-1,0)}(\dot{r}_2^{i\backslash})\right)+{\mathbb{V}}^{(n-1,1)}_{\rm FIN}(\dot{r}^{i\backslash})\,. \end{align}\] The processes \(\dot{r}^{i\backslash}\), \(\dot{r}_1^{i\backslash}\) and \(\dot{r}_2^{i\backslash}\) are eq.@eq:eq:BornCO3Pwave .

References↩︎

---

1. Strictly speaking, this statement only applies to processes with not-too-high particle multiplicity, given the limitations of computing resources. It is also possible that special issues may arise for particular problems that have not yet been addressed in the automated codes. Additionally, this consideration does not cover processes involving loops at their lowest order. In the context of this paper, when we refer to a similar statement, we mean it in a loose sense.↩︎

2. Thanks to recent advancements in parton showers [11], there is a chance that the problem of giant \(K\) factors can be addressed through the matching and merging of matrix elements and parton showers.↩︎

3. In simple cases, such as \(2\rightarrow 1\) or \(1\rightarrow 2\) underlying Born processes, fully analytical \(d\)-dimensional phase space integrations of real emissions are possible, similar to the approach taken in ref. [63].↩︎

4. The formalism was derived under the circumstance of a \(2\rightarrow 2\) underlying Born process.↩︎

5. Note that the flavours of the two constituent heavy quarks do not have to be identical. For instance, the constituent quarks of \(B_c^+\) are a charm quark and a bottom antiquark.↩︎

6. We do not restrict the numbers of \(Q\) and \(\bar{Q}^\prime\) appearing in the process.↩︎

7. The earlier analysis of the soft limit tailored for specific processes can be found in the literature, such as in sect. 4 of ref. [63] for \(2\rightarrow 1\) processes.↩︎

8. \(K_{ab}(z)\) is trivially zero for \(\overline{\rm MS}\) PDFs.↩︎

9. Note that \(E_1=E_2=\sqrt{s}/2\) here.↩︎

10. The analytic results for majority processes, excluding \(gg\rightarrow c\bar{c}[{\bigl.^1P^{[8]}_1}]+X\), can also be found in ref. [63].↩︎

11. Note that the factorisation-breaking effect due to colour transfer between a quarkonium and a heavy quark [74] only occurs in a particular phase-space corner, i.e.  the mass threshold region, where the relative velocity of the heavy quark and the quarkonium tends to zero. Therefore, it would be harmless for our general purpose unless we are probing that particular phase-space region, where, in any case, we need additional theoretical care.↩︎

12. Note that in the case \(l=k\) the derivative acts on both momenta.↩︎