1 Introduction↩︎

Actual and planned CERN LHC operation opens the possibility to perform a large number of accurate measurements in high energy physics. It is clear that for many of them the overall experimental uncertainty will be much smaller than the theory uncertainty estimated based on next-to-leading order QCD corrections. Then, including higher order corrections turns out to be mandatory for a meaningful comparison among theory predictions and experimental data. Thus, next-to-next-to-leading order (NNLO) computation has received considerable attention and several approaches for performing these calculations have been proposed. These include the \(q_T\) [1] and \(N\)-jettiness [2], [3] slicing methods, as well as subtraction schemes, like antenna [4]–[12], CoLoRFulNNLO [13]–[22], residue-improved [23]–[26], nested soft-collinear [27]–[31] and projection-to-Born [32] with yet other approaches under development [33]–[35].

From the mathematical point of view, computations at NNLO are more elaborate than ones at NLO and for this reason the level of automation is still much less advanced. On the one hand, difficulties lie in the computation of the double virtual amplitudes for processes with many particles in the final state and with masses. On the other hand, there are still no fully satisfactory, complete and general algorithms for the regularization of the infrared divergences for fully differential NNLO computation as there are for the NLO case. Here by general, we refer to a scheme that applies to any kind of singularity coming from both initial and massive or massless final state particles. By complete, we mean that the set of subtractions defining a scheme is given in full detail together with the complete integrated versions. This last point is of course mandatory to allow for independent applications, validations or just implementations for analysis purposes.

Here we consider the CoLoRFulNNLO method that has been formulated in ref. [13]–[22]. The basic ideas of the method are quite general and apply to final state as well as initial state singularities. We focus on the final state infrared singularities and in particular, on the production of a massive quark-antiquark pair from a colourless initial state. NNLO corrections to such processes have been computed in the literature previously [36]–[40] but still, a complete procedure for the local subtraction of all infrared singularities, supplemented by a compact analytic expression for the sum of the integrated counterterms has been missing. In the present paper we fill this gap by providing a complete scheme for the pair creation of a heavy quark-antiquark pair out of a colourless initial state.

The paper is organized as follows. In section 2 we define the problem we want to address and set up our notation. In section 3, we compute the NLO correction to heavy quark production from a general colourless initial state. We introduce the set of NNLO counterterms in section 4, where we also present a complete and compact analytic expression for the integrated form of the sum of all counterterms. In section 5, we present the application of our subtraction scheme to the production of heavy quarks in the decay of the Standard Model Higgs boson. Finally, we draw our conclusions in section 6.

2 Heavy quark-antiquark pair production from colourless initial states↩︎

Let us consider the production of a heavy quark-antiquark pair from a generic colourless initial state \(X\). The list of relevant subprocesses up to NNLO accuracy in QCD is given by

In the above list the heavy quarks are denoted by \(Q\) and \(\bar{Q}\), while radiated gluons and light quarks are denoted by \(g\) and \(q{\bar q}\) respectively. The matrix elements for all of these partonic processes are well-known for scalar, pseudo-scalar, as well as vector and axial currents. In particular the two-loop form factors were first computed in [41]–[44] up to finite terms in the parameter of dimensional regularization \(\epsilon\). Recently these results have been extended up to \(\mathrm O(\epsilon^{1})\) terms in ref. [45].

For amplitudes, we use the colour-space notation of [46] where \(|{\cal M}_{}^{}\rangle\) is an abstract vector in colour- and spin-space, such that the matrix element summed over colours and spins can be written as \[|{\cal M}_{}^{}{}|^2 = \langle{\cal M}_{}^{}|{\cal M}_{}^{}\rangle\,.\] Then, the insertion of colour charge operators, such as \({\boldsymbol{T}}_i {\boldsymbol{T}}_k\) or the anti-commutator \(\{{\boldsymbol{T}}_i {\boldsymbol{T}}_k, {\boldsymbol{T}}_j {\boldsymbol{T}}_l\}\), as well as polarization-dependent tensors such as the Altarelli–Parisi splitting kernels \(\hat{P}_{f_i f_r}\) (to be specified below), into a given amplitude interference will be denoted as _i \(T\)_k |M_^|^2 & _^| \(T\)_i \(T\)_k |M_^ ,
{\(T\)_i \(T\)_k, \(T\)_j \(T\)_l} |M_^|^2 & _^| {\(T\)_i \(T\)_k, \(T\)_j \(T\)_l} |M_^ ,
_f_i f_r |M_^|^2 & _^| _f_i f_r |M_^ .

The amplitude \(|{\cal M}_{}^{}\rangle\) has the formal loop expansion \[|{\cal M}_{}^{}\rangle = |{\cal M}_{}^{(0)}\rangle + |{\cal M}_{}^{(1)}\rangle + |{\cal M}_{}^{(2)}\rangle + \ldots\,.\] We consider amplitudes computed in conventional dimensional regualrization with the strong coupling \(\alpha_{\mathrm{s}}\) and gluon wave function renormalized in the \(\overline{\mathrm{MS}}\) scheme, while the heavy quark mass and wave function are renormalized on-shell. Furthermore, assuming \(n_l\) light flavours, plus a single heavy flavour \(Q\), we implement the transition that allows us to use \(n_l+1\) active flavours in the running strong coupling [47]. In particular, \[\alpha_{\mathrm{s}}^B \mu_0^{2\epsilon} = \frac{\alpha_{\mathrm{s}}\mu_R^{2\epsilon}}{C(\epsilon)} \left[1 - \frac{\alpha_{\mathrm{s}}}{4\pi} \frac{\beta_0}{\epsilon} + \mathrm O(\alpha_{\mathrm{s}}^{2})\right]\,, \label{eq:aS-ren}\tag{1}\] where \(\alpha_{\mathrm{s}}^B\) is the bare coupling, while \(\alpha_{\mathrm{s}}\) denotes the renormalized coupling which will appear in all subsequent equations. Furthermore, the beta-function coefficient \(\beta_0\) reads \[\beta_0 = \frac{11}{3}C_{\mathrm{A}}- \frac{4}{3}T_{\mathrm{R}}(n_l+1)\] and \[C(\epsilon) = (4\pi)^{\epsilon}\Gamma(1+\epsilon)\,. \label{eq:Cep-def}\tag{2}\] Note that \(C(\epsilon)\) as defined above is different form the usual expression in the \(\overline{\mathrm{MS}}\) scheme of \(S_\epsilon^{\overline{\mathrm{MS}}} = (4\pi)^\epsilon \exp(-\epsilon\gamma_E)\) and agrees with the convention of refs. [41]–[44]. However, for the processes considered here the perturbative expansion starts at \(\alpha_{\mathrm{s}}^0\), therefore the inclusion of the NNLO corrections implies just one-loop renormalization for the coupling constant. For this reason, the \(\mathrm O(\epsilon^{2})\) difference between these two conventions turns out to have no impact on the physical result.

Throughout we denote by \(P\) the total incoming momentum of the process and make use of the following definitions \[s_{ij} \equiv 2p_i \cdot p_j\,, \qquad y_{ij} \equiv \frac{s_{ij}}{P^2} \qquadand\qquad y_{(ij)k} \equiv y_{ik} + y_{jk}\,.\] We also employ this notation for mapped momenta (see eqs. (7 ) and (11 ) below), so that e.g., \[s_{\widehat{\,i\,}\widehat{\,k\,}} \equiv 2\widehat{\,p\,}_i\cdot \widehat{\,p\,}_k\,, \qquad s_{\widehat{\,i\,}r} \equiv 2\widehat{\,p\,}_i\cdot p_r \qquadand\qquad s_{\widetilde{\,i\,}\widetilde{\,k\,}} \equiv 2\widetilde{\,p\,}_i\cdot \widetilde{\,p\,}_k\,.\] The phase space of \(n\) outgoing particles of total momentum \(P\) is defined as \[{\mathrm{d}}\phi_{n}(p_1,\ldots,p_n;P) = \prod_{i=1}^{n} \frac{{\mathrm{d}}^d p_i}{(2\pi)^{d-1}} \delta_{+}(p_i^2 - m_i^2) (2\pi)^d \delta^{(4)}(p_1+\ldots+p_n-P)\,.\]

Finally, integrated subtractions are given in terms of multiple polylogarithms (\(G\)), which can be defined recursively by the iterated integral [48], [49] \[G_{a_1,\ldots,a_n}(y) = \int_0^y {\mathrm{d}}t\, \frac{1}{t - a_1}G_{a_2,\ldots,a_n}(t) \label{eq:G-def}\tag{3}\] with \(G(y)=1\). In the special case where all the \(a_i\)’s are zero, the multiple polylogarithm is defined as \[G_{0_1,\ldots,0_n}(y) = \frac{1}{n!}\ln^n y\,,\] which is consistent with \(G(y)=1\) for \(n=0\). For completeness we note that for \(a_i\in\{0,\pm 1\}\), the \(G\)’s are related to the harmonic polylogarithms (\(H\)) of [50] by the relation \[G_{a_1,\ldots,a_n}(y) = (-1)^p H_{a_1,\ldots,a_n}(y)\,, \qquad a_i\in\{0,\pm 1\}\,\] where \(p\) denotes the number of \(a_i\)’s that are equal to \(+1\).

3 Subtractions at NLO↩︎

Consider the NLO correction to the production of the heavy quark-antiquark pair, which is the sum of the real contribution involving the emission of an extra gluon and the virtual contribution containing the one-loop correction, \[\Gamma^{\mathrm{NLO}}[J] = \int_3 {\mathrm{d}}\Gamma^{{\rm R}} J_3 + \int_2 {\mathrm{d}}\Gamma^{{\rm V}} J_2\,. \label{eq:Gamma-NLO}\tag{4}\] Here \(\int_3\) and \(\int_2\) denote the integration over the \(Q\bar{Q}g\) and \(Q\bar{Q}\) phase space, while \(J_3\) and \(J_2\) are the values of some infrared-safe observable \(J\) computed with the corresponding 3- and 2-parton kinematics. By introducing an appropriate local subtraction term to regulate infrared divergences, we rewrite eq. (4 ) as a sum of two finite terms \[\Gamma^{\mathrm{NLO}}[J] = \int_3 {\mathrm{d}}\Gamma^{{\rm NLO}}_3 + \int_2 {\mathrm{d}}\Gamma^{{\rm NLO}}_2\,,\] with the regularized real and regularized virtual contributions¹ given by _3 &= ^R J_3 - ^R,A_ J_2 ,
^NLO_2 &= J_2 . Above \(\int_{[1]}\) denotes the integration of the subtraction terms over the radiation variables of the extra gluon.

3.1 Regularized real contribution↩︎

Let us consider first the real emission process, \(X(P) \to Q(p_1) + \bar{Q}(p_2) + g(p_3)\). Denoting by \(F\) the flux factor², we have \[{\mathrm{d}}\Gamma^{{\rm R}} = \frac{1}{F} {\mathrm{d}}\phi_{3}(p_1,p_2,p_3;P) |{\cal M}_{Q\bar{Q} g}^{(0)}{}|^2\] and \[{\mathrm{d}}\Gamma^{{\rm R,A}_{\scriptscriptstyle 1}} = \frac{1}{F} {\mathrm{d}}\phi_{3}(p_1,p_2,p_3;P) {\cal A}_{1} |{\cal M}_{Q\bar{Q} g}^{(0)}{}|^2\,.\] We note that throughout this subsection, \(p_1\), \(p_2\) and \(p_3\) refer to the momenta of the heavy quark \(Q\), the heavy antiquark \(\bar{Q}\) and the gluon \(g\) in the three-particle real emission phase space. The matrix element is singular only in the \(p_3^\mu \to 0\) soft gluon limit, thus the structure of the approximate matrix element is very simple and we have just one subtraction term, \[{\cal A}_1 |{\cal M}_{Q\bar{Q} g}^{(0)}{}|^2 = {\cal S}_{g_3}^{(0)}\,.\]

3.1.0.1 Single soft subtraction.

The subtraction term follows the structure of the general formula for the approximation of a tree-level \(n+1\)-parton matrix element in the \(p_r^\mu \to 0\) soft limit [51] and is given by \[{\cal S}_{g_r}^{(0)}(p_1,p_2,p_3) \equiv -8\pi\frac{\alpha_{\mathrm{s}}\mu^{2\epsilon}}{C(\epsilon)} \sum_{\widehat{\,i\,},\widehat{\,k\,}} \frac{1}{2}S_{\widehat{\,i\,}\widehat{\,k\,}}(r)\, {\boldsymbol{T}}_{\widehat{\,i\,}} {\boldsymbol{T}}_{\widehat{\,k\,}} \otimes |{\cal M}_{n}^{(0)}{(\widehat{\,p\,}_1,\widehat{\,p\,}_2)}|^2 \,, \label{eq:S3g}\tag{5}\] where the summation indices \(\widehat{\,i\,}\) and \(\widehat{\,k\,}\) run over the labels of the hard momenta that appear in the factorized matrix element (i.e., \(\widehat{\,1\,}\) and \(\widehat{\,2\,}\), see below), and the eikonal factor is also computed using these momenta, \[S_{\widehat{\,i\,}\widehat{\,k\,}}(r) = \frac{2 s_{\widehat{\,i\,}\widehat{\,k\,}}}{s_{\widehat{\,i\,}r}s_{\widehat{\,k\,}r}} = \frac{(\widehat{\,p\,}_i\cdot \widehat{\,p\,}_k)}{(\widehat{\,p\,}_i\cdot p_r)(\widehat{\,p\,}_k\cdot p_r)} \quadand\quad S_{\widehat{\,i\,}\widehat{\,i\,}}(r) = \frac{2 s_{\widehat{\,i\,}\widehat{\,i\,}}}{s_{\widehat{\,i\,}r}^2} = \frac{\widehat{\,p\,}_i^2}{(\widehat{\,p\,}_i\cdot p_r)^2} \,. \label{eq:SiHrkHr}\tag{6}\]

The \(3\to 2\) mapping \(\{p_1,p_2,p_3\} \to \{\widehat{\,p\,}_1,\widehat{\,p\,}_2\}\) that specifies the hatted momenta which enter the factorized matrix element in eq. (5 ) is defined as follows (to be clear, the heavy quark \(Q\) and anitquark \(\bar{Q}\) carry momenta \(\widehat{\,p\,}_1\) and \(\widehat{\,p\,}_2\) in the two-parton matrix element on the right hand side of eq. (5 )), \[\esp{p}_{1}^\mu &= \Lambda^{\mu}_{\nu}(\widehat{\,K\,},K) \frac{1}{\Delta} (p_1^\mu - \gamma_1 p_{3}^\mu)\,, \\ \widehat{\,p\,}_{2}^\mu &= \Lambda^{\mu}_{\nu}(\widehat{\,K\,},K) \frac{1}{\Delta} (p_2^\mu - \gamma_2 p_{3}^\mu)\,, \esp \label{eq:p1p2-Cirmap}\tag{7}\] where \(\Delta\) and \(\gamma_{1,2}\) are \[\Delta = \sqrt{1 - \frac{s_{13} s_{23} s_{3P}}{P^2 s_{13} s_{23} - m_Q^2 s_{3P}^2}}\,,\] and \[\gamma_1 = \frac{m_Q^2 s_{23} s_{3P}}{P^2 s_{13} s_{23} - m_Q^2 s_{3P}^2}\,, \qquad \gamma_2 = \frac{m_Q^2 s_{13} s_{3P}}{P^2 s_{13} s_{23} - m_Q^2 s_{3P}^2}\,.\] With these definitions we have \(\widehat{\,p\,}_1^2=\widehat{\,p\,}_2^2=m_Q^2\) and \(\widehat{\,K\,}^2 = K^2\), where \[K^\mu = P^\mu \qquadand\qquad \widehat{\,K\,}^\mu = \frac{1}{\Delta}[p_1^\mu + p_2^\mu - (\gamma_1+\gamma_2)p_3^\mu]\,.\] Finally, \(\Lambda^{\mu}_{\nu}(\widehat{\,K\,},K)\) is a (proper) Lorentz transformation that takes \(K^\mu\) into \(\widehat{\,K\,}^\mu\), its explicit form can be chosen e.g., as \[\Lambda^{\mu}_{\nu}(K,\widehat{\,K\,}) = g^{\mu}_{\nu} - \frac{2(K+\widehat{\,K\,})^\mu (K+\widehat{\,K\,})_\nu}{(K+\widehat{\,K\,})^2} + \frac{2 \widehat{\,K\,}^\mu K_\nu}{K^2}\,. \label{eq:Lambda}\tag{8}\]

3.2 Regularized virtual contribution↩︎

Turning to the two terms in eq. ([eq:dGammaNLO2]), the virtual contribution \({\mathrm{d}}\Gamma^{{\rm V}}\) involves the one-loop correction to the process \(X(P) \to Q(p_1) + \bar{Q}(p_2)\) and we have \[{\mathrm{d}}\Gamma^{{\rm V}} = \frac{1}{F} {\mathrm{d}}\phi_{2}(p_1,p_2;P) 2\Re\langle{{\cal M}}_{Q\bar{Q}}^{(0)}|{\cal M}_{Q\bar{Q}}^{(1)}\rangle\] with \[\int_{[1]} {\mathrm{d}}\Gamma^{{\rm R,A}_{\scriptscriptstyle 1}} = \frac{1}{F} {\mathrm{d}}\phi_{2}(p_1,p_2;P)\, {\boldsymbol{I}}_{1}(p_1,p_2;\epsilon) \otimes |{\cal M}_{Q\bar{Q}}^{(0)}{}|^2\] where the \({\boldsymbol{I}}_{1}(p_1,p_2;\epsilon)\) operator corresponds to the integral of the only subtraction term, \({\cal S}_{g_3}^{(0)}\). We remark that throughout this subsection, \(p_1\) and \(p_2\) denote the momenta of the heavy quark \(Q\) and antiquark \(\bar{Q}\) in the two-body phase space.

After resolving the summation in eq. (5 ) using the colour algebra relations \({\boldsymbol{T}}_1^2 = {\boldsymbol{T}}_2^2 = C_{\mathrm{F}}\) and \({\boldsymbol{T}}_1{\boldsymbol{T}}_2 = -C_{\mathrm{F}}\), the insertion operator can be computed in a straightforward way by integrating \({\cal S}_{g_r}^{(0)}\) over the full three body phase space and simply divide by the volume of the (massive) two particle phase space, which is a constant. The result of the integration in can be cast in the form \[{\boldsymbol{I}}_{1}(p_1,p_2;\epsilon) = \frac{\alpha_{\mathrm{s}}}{2\pi} \left(\frac{\mu^2}{P^2}\right)^\epsilon C_{\mathrm{F}}\left(\frac{1}{\epsilon} a_{-1} + a_0 + \epsilon\, a_1 + \mathrm O(\epsilon^{2})\right)\,, \label{eq:I10-a-coeffs}\tag{9}\] where the coefficients of the Laurent expansion are functions of the customary variable \[y \equiv \frac{\sqrt{P^2} - \sqrt{P^2 - 4m_Q^2}}{\sqrt{P^2} + \sqrt{P^2 - 4m_Q^2}}\,, \label{eq:y-def}\tag{10}\] which is real for the physical decay process. The coefficients appearing in eq. (9 ) above can be expressed in terms of multiple polylogarithms defined in eq. (3 ), always with argument \(y\). Setting \(G_{a_1,\ldots,a_n}(y) = G_{a_1,\ldots,a_n}\) for ease of notation, we have _-1 &= 2 + G_0 ,
a_0 &= 4 -8 G_1 -G_0 +(4 G_-1,0-G_0,0-4 G_0,1+2 G_1,0-4 _2) ,
a_1 &= 8 -16 G_1 +32 G_1,1 -  . For the NLO cross section only the first two terms in the expansion of \({\boldsymbol{I}}_{1}(p_1,p_2;\epsilon)\) are needed. Nevertheless, the order \(\epsilon\) term will enter the integrated subtraction for the single unresolved limit of real-virtual emission at NNLO and so we present it here.

Note that in eq. (9 ) we have expanded the factor of \(C(\epsilon)\) that appears in the denominator of eq. (5 ) which cancels, as usual, terms of \(\gamma_E\) and \(\ln(4\pi)\) coming from phase space integration. If the strong coupling is defined with a different \(\epsilon\)-dependent prefactor in eq. (1 ), the explicit forms of the expansion coefficients change accordingly. In particular, adopting the standard \(\overline{\mathrm{MS}}\) factor of \(S_\epsilon^{\overline{\mathrm{MS}}} = (4\pi)^\epsilon\exp(-\epsilon\gamma_E)\), we would have \[a_1^{\overline{\mathrm{MS}}} = a_1 + \frac{\zeta_{2}}{2} a_{-1}\,,\] while the lower order expansion coefficients remain unchanged.

4 Subtractions at NNLO↩︎

The NNLO correction is the sum the double real, real-virtual and double virtual parts, \[\Gamma^{\mathrm{NNLO}}[J] = \int_4 {\mathrm{d}}\Gamma^{{\rm RR}} J_4 + \int_3 {\mathrm{d}}\Gamma^{{\rm RV}} J_3 + \int_2 {\mathrm{d}}\Gamma^{{\rm VV}} J_2\,,\] which we rearrange into three finite contributions by including appropriate subtraction terms, \[\Gamma^{\mathrm{NNLO}}[J] = \int_4 {\mathrm{d}}\Gamma^{{\rm NNLO}}_4 + \int_3 {\mathrm{d}}\Gamma^{{\rm NNLO}}_3 + \int_2 {\mathrm{d}}\Gamma^{{\rm NNLO}}_2\,.\] The regularized double real, regularized real-virtual and regularized double virtual contributions are _4 &= ^RR J_4 - ^RR,A_ J_3 - ^RR,A_ J_2 + ^RR,A_ J_2
^NNLO_3 &= J_3 - J_2
^NNLO_2 &= J_2 . Above \(\int_{[1]}\) and \(\int_{[2]}\) denote the integration of subtraction terms over the radiation variables of one and two extra partons.

4.1 Regularized double real contribution↩︎

Considering all possible subprocesses with one heavy flavour \(Q\) and \(n_l\) massless flavours \(q\), the sum of all such contributions reads \[{\mathrm{d}}\Gamma^{{\rm RR}} = \frac{1}{F} {\mathrm{d}}\phi_{4}(p_1,p_2,p_3,p_4;P) \left[\frac{1}{2} |{\cal M}_{Q\bar{Q} gg}^{(0)}{}|^2 +n_l |{\cal M}_{Q\bar{Q} q{\bar q}}^{(0)}{}|^2 +\frac{1}{4} |{\cal M}_{Q\bar{Q} Q\bar{Q}}^{(0)}{}|^2\right]\,,\] where we have explicitly reported the statistical factors in front of the matrix elements for the production of two gluons and the one for the production of two heavy quark-antiquark pairs. We label the particles such that the heavy quark \(Q\) and heavy antiquark \(\bar{Q}\) always carry momenta \(p_1\) and \(p_2\), while \(p_3\) and \(p_4\) are the momenta associated with the extra emissions (either two gluons \(gg\), a light quark-antiquark pair \(q{\bar q}\), or one more heavy quark and antiquark, \(Q\bar{Q}\)). We emphasize that \(p_1\), \(p_2\), \(p_3\) and \(p_4\) denote momenta in the four-particle double real emission phase space throughout this subsection. The subtraction terms introduced in eq. ([eq:dGNNLO4]) are³ &= _4(p_1,p_2,p_3,p_4;P)  ,
^RR,A_ &= _4(p_1,p_2,p_3,p_4;P)  ,
^RR,A_ &= _4(p_1,p_2,p_3,p_4;P)  .

Starting with the \(X(P) \to Q(p_1) + Q(p_2) + g(p_3) + g(p_4)\) subprocess, the \(|{\cal M}_{Q\bar{Q} gg}^{(0)}{}|^2\) matrix element requires regularization by subtraction only in the following infrared limits:

1. the \(p_3^\mu || p_4^\mu\) single unresolved collinear limit

2. the \(p_3^\mu \to 0\) soft limit

3. the \(p_4^\mu \to 0\) soft limit

4. the \(p_3^\mu \to 0\), \(p_4^\mu \to 0\) double soft gluon limit

Hence, the single, double and iterated single unresolved approximate matrix elements for this subprocess have the following structure, _1 |M_Q|Q gg^(0)|^2 &= C_g_3 g_4^(0) F^C_34 + (S_g_3^(0) - C_g_3 g_4^S_g_3^(0)) F^S_3 + (S_g_4^(0) - C_g_3 g_4^S_g_4^(0)) F^S_4 ,
A_2 |M_Q|Q gg^(0)|^2 &= S_g_3 g_4^(0) ,
A_12 |M_Q|Q gg^(0)|^2 &= C_g_3 g_4^S_g_3 g_4^(0) F^C_34 + (S_g_3^S_g_3 g_4^(0) - C_g_3 g_4^S_g_3^S_g_3 g_4^(0)) F^S_3 + (S_g_4^S_g_3 g_4^(0) - C_g_3 g_4^S_g_4^S_g_3 g_4^(0)) F^S_4 . In eq. ([eq:A1bbgg]), the subtraction terms \({\cal C}_{g_3 g_4}^{}{\cal S}_{g_3}^{(0)}\) and \({\cal C}_{g_3 g_4}^{}{\cal S}_{g_4}^{(0)}\) are included to avoid double subtraction over those regions of phase space where the collinear and soft limits overlap. Moreover, note that formally \({\cal A}_{12} = {\cal A}_{1} {\cal A}_{2}\), i.e., the form of the iterated single unresolved approximate matrix element agrees with that of the single unresolved approximate matrix element.

Although the structure of eqs. ([eq:A1bbgg])–([eq:A12bbgg]) is dictated by the types of infrared limits which require regularization, the explicit definition of the subtraction terms is obviously not unique. Different choices can have various advantages and drawbacks (e.g., locality of subtractions versus full analytic control over the integrated subtraction terms). In particular, a general issue for any subtraction scheme at NNLO concerns the integration of counterterms, which can turn out to be a very elaborate task. Thus, on practical grounds, once the general structure of the counterterms is defined and momentum conservation has been implemented, one may seek to exploit the freedom in the definitions of counterterms to simplify the integration. This consideration motivates the inclusion of the collinear factor \({\cal F}^{C}_{34}\) and the soft factors \({\cal F}^{S}_{3}\) and \({\cal F}^{S}_{4}\) in the above formulae⁴. Clearly, the collinear and soft factors must go to the identity in the corresponding limit. Furthermore, to preserve the structure of cancellations among the subtraction terms in all limits, we find that in our construction the soft-collinear overlap must be multipiled with the soft factor, while \({\cal A}_{12}|{\cal M}_{Q\bar{Q} gg}^{(0)}{}|^2\) inherits the pattern of modifications of \({\cal A}_1|{\cal M}_{Q\bar{Q} gg}^{(0)}{}|^2\). In the following, we present a concrete example of a constructive procedure for obtaining factors that lead to a fully analytic result for the sum of all integrated subtraction terms which is very compact, see section 4.3.

Turning to the \(X(P) \to Q(p_1) + \bar{Q}(p_2) + q(p_3) + {\bar q}(p_4)\) subprocess, the only infrared limits of \(|{\cal M}_{Q\bar{Q} q{\bar q}}^{(0)}{}|^2\) that require regularization by subtraction are:

1. the \(p_3^\mu || p_4^\mu\) single collinear limit

2. the \(p_3^\mu \to 0\), \(p_4^\mu \to 0\) double soft quark-antiquark limit

Correspondingly, the structure of the subtractions is very simple and each approximate matrix element is built form a single term, _1 |M_Q|Q q|q^(0)|^2 &= C_q_3 |q_4^(0) F^C_34 ,
A_2 |M_Q|Q q|q^(0)|^2 &= S_q_3 |q_4^(0) ,
A_12 |M_Q|Q q|q^(0)|^2 &= C_q_3 |q_4^S_q_3 |q_4^(0) F^C_34 , As previously, \({\cal A}_{12} = {\cal A}_{1} {\cal A}_{2}\) formally.

Before presenting the explicit expressions of each subtraction term, let us first discuss the kinematics and in particular the momentum mappings used to enforce exact phase space factorization. The definition of subtraction terms involves the specification of functions which map the double real emission phase space into phase spaces of lower multiplicity plus radiation variables. In particular, we find that all single unresolved subtraction terms can be defined using just one \(4\to 3\) momentum mapping. The mapping appropriate to the double unresolved and iterated single unresolved subtractions is then obtained simply by applying the \(3\to 2\) mapping discussed in section 3.1 to the output of the \(4\to 3\) mapping presented below. Given momenta \(\{p_1,p_2,p_3,p_4\}\) where \(p_1^2=p_2^2=m_Q^2\) and \(p_3^2=p_4^2=0\), to be mapped to \(\{\widehat{\,p\,}_1,\widehat{\,p\,}_2,\widehat{\,p\,}_{34}\}\) with \(\widehat{\,p\,}_1^2=\widehat{\,p\,}_2^2=m_Q^2\) and \(\widehat{\,p\,}_{34}^2=0\), we set \[\esp{p}_{34}^\mu &= 1/\beta (p_3^\mu + p_4^\mu - \alpha P^\mu)\,, \\ \widehat{\,p\,}_n^\mu &= \Lambda^{\mu}_{\nu}(K,\widehat{\,K\,}) p_{n}^{\nu}\,, \qquad\qquad\qquad n=1,2\,, \label{eq:4to3map} \esp\tag{11}\] where \(\alpha\) and \(\beta\) are \[\alpha = \frac{1}{2} \left[y_{(34)P} - \sqrt{y_{(34)P}^2 - 4y_{34}}\,\right] \qquadand\qquad \beta = \frac{\sqrt{y_{(34)P}^2 - 4y_{34}}}{y_{(34)P} - y_{34}}\,. \label{eq:al-be-def}\tag{12}\] With these definitions \(\widehat{\,p\,}_{34}\) is massless and the momenta \[K^\mu = P^\mu - p_3^\mu - p_4^\mu \qquadand\qquad \widehat{\,K\,}^\mu = P^\mu - \widehat{\,p\,}_{34}^\mu\] have the same mass, \(K^2 = \widehat{\,K\,}^2\). Hence, \(\Lambda^{\mu}_{\nu}(\widehat{\,K\,},K)\) is a (proper) Lorentz transformation that takes \(K^\mu\) into \(\widehat{\,K\,}^\mu\), whose explicit form can be chosen as in eq. (8 ). We note that this mapping is equivalent to the final state mapping presented in ref. [53].

The momentum mapping introduced above leads to the exact factorization of the four particle phase space in the following form, \[{\mathrm{d}}\phi_{4}(p_1,p_2,p_3,p_4;P) = {\mathrm{d}}\phi_{3}(\widehat{\,p\,}_1,\widehat{\,p\,}_2,\widehat{\,p\,}_{34};P)[{\mathrm{d}}p]\,,\] where the measure for the factorized radiation variables \([{\mathrm{d}}p]\) reads \[[{\mathrm{d}}p] = x^{-1+2\epsilon} \frac{P^2}{2\pi} \int_{\alpha_{\rm min}}^{\alpha_{\rm max}} {\mathrm{d}}\alpha\, (1-\alpha)^{-3+2\epsilon} (x - 2\alpha+ \alpha^2)^{2-2\epsilon} {\mathrm{d}}\phi_{2}(p_3,p_4; \alpha P + \beta\widehat{\,p\,}_{34})\,, \label{eq:dp}\tag{13}\] with \[x \equiv y_{\widehat{34}P} = \frac{2\widehat{\,p\,}_{34} \cdot P}{P^2}\,, \label{eq:x-def}\tag{14}\] and expressing \(\beta\) of eq. (12 ) in terms of \(\alpha\) and \(x\) we find \[\beta= \frac{x - 2\alpha+ \alpha^2}{x(1-\alpha)}\,. \label{eq:beta-alpha}\tag{15}\] Since \(\widehat{\,p\,}_{34}^2=0\) we have \[\alpha_{\rm min} = 0 \qquadand\qquad \alpha_{\rm max} = 1 - \sqrt{1-x}\,. \label{eq:al-min-max}\tag{16}\]

Before going on, let us anticipate some difficulties which appear when integrating the single unresolved subtraction terms over the measure in eq. (13 ). First, the definition of the collinear subtraction term involves the specification of a momentum fraction \(z_{r}\) (\(r=3,4\)) associated with the splitting. However, a natural candidate for this variable, \(z_{r} = \frac{p_{r}\cdot P}{(p_3+p_4)\cdot P}\), turns out to be a somewhat complicated function of the radiation variables. Second, the soft subtraction term involves the eikonal factor with the hard momenta \(\widehat{\,p\,}_1\) and \(\widehat{\,p\,}_2\). In addition, the measure \([dp]\) evidently depends also on \(\widehat{\,p\,}_{34}\). Thus, the result of the integration will depend on all independent dot-products between these three vectors in a very complicated way. Last, the upper limit in eq. (16 ) is a square root function of the invariant \(x\), which implies that the integrated counterterms will also be function of this square root. Regardless of the first two issues, this last point alone leads to difficulties when computing the iterated single unresolved subtraction terms.

However, exploiting the freedom in the definition of the subtraction terms, one can devise a strategy to tackle the above mentioned difficulties. First, a more convenient choice of momentum fractions can be made upon examination of the explicit form of the collinear integral without affecting the structure of singularities. Second, as anticipated in eq. (5 ) we multiply the soft integral by an appropriate function \({\cal F}^S_r\) that reduces to one in the unresolved limit in \(d\) dimensions. This function can be chosen in such a way that it cancels regular factors, effectively reducing the multiple angular dependence of the integrand. Last, one can restrict the phase space of the subtraction⁵. This restriction can be implemented by adopting an appropriate functional form of the upper limit of integration with respect to \(x\) (e.g., linear) which avoids the dependence on square roots of invariants.

4.1.0.1 Single collinear subtraction.

In order to define the collinear subtraction term, we start from the well-known approximation to the matrix element in this limit [57] \[|{\cal M}_{Q\bar{Q} f_3 f_4}^{(0)}{(p_1,p_2,p_3,p_4)}|^2 \simeq 8\pi\frac{\alpha_{\mathrm{s}}\mu^{2\epsilon}}{C(\epsilon)} \frac{1}{s_{34}} \hat{P}_{f_3 f_4}(z_3, k_{\perp};\epsilon) \otimes |{\cal M}_{Q\bar{Q} g}^{(0)}{(p_1, p_2, p_3+p_4)}|^2 \label{eq:C34-limit}\tag{17}\] where \(\hat{P}_{f_3 f_4}(z_3, k_{\perp};\epsilon)\) is the \(d\)-dimensional Altarelli–Parisi splitting kernel for the \(f_{(34)} \to f_3+f_4\) splitting (here \(f\) denotes the parton flavour) that are functions of the momentum fraction (\(z_3\)) and the transverse momentum (\(k_{\perp}\)) of the splitting. For our calculation only gluon splitting is relevant, for which the kernels are given explicitly by ^_gg(z_3, k_) &= 2 C_ ,
P^_q|q(z_3, k_) &= T_ .

To build a proper subtraction counterterm from the above limit formula, as usual we need to evaluate the factorized matrix element on the right hand side with mapped momenta, that respect momentum conservation and the mass shell conditions. Furthermore, the momentum fractions and \(k_{\perp}\) must be properly defined over the full phase space. A straightforward choice for \(z_3\) would read \[z_{3} = \frac{p_3 \cdot P}{(p_3+p_4)\cdot P} \qquadand\qquad 1-z_{3} = \frac{p_4 \cdot P}{(p_3+p_4)\cdot P}\,. \label{eq:Z34Z43}\tag{18}\] Although it is a simple exercise to construct the subtraction term in this way, using the momentum mapping of eq. (11 ), it turns out that the integrated form of this subtraction is rather cumbersome. In order to exhibit the reasons behind this, we recall the measure for the radiation variables, eq. (13 ), and note that the two-particle phase space \({\mathrm{d}}\phi_{2}(p_3,p_4; \alpha P + \beta\widehat{\,p\,}_{34})\) appearing there can be parametrised as follows, \[{\mathrm{d}}\phi_{2}(p_3,p_4;\alpha P + \beta\widehat{\,p\,}_{34}) = \frac{1}{8\pi} \frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)} (P^2)^{-\epsilon} {\mathrm{d}}v\, \alpha^{-\epsilon} (\alpha+ \beta x)^{-\epsilon} v^{-\epsilon} (1-v)^{-\epsilon} \Theta(v)\Theta(1-v)\,, \label{eq:dphi2}\tag{19}\] where \(\beta\) is given in eq. (15 ), while \(v\) is defined implicitly by the following relation, \[z_{3} = \frac{\alpha(1-\alpha) + (x - 2\alpha+ \alpha^2) v}{x - \alpha^2}\,.\] Notice that \(1-z_{3}\) is obtained by \(v \to 1-v\) in the above expression. Furthermore, in this parametrization the two-paricle invariant \(y_{34}\) reads \[y_{34} = \frac{\alpha(x - \alpha)}{1-\alpha}\,,\] so the colliner limit (when \(y_{34}\to 0\)) corresponds to \(\alpha\to 0\) (note that \(\alpha_{\rm max} = 1 - \sqrt{1-x}\) so \(x-\alpha> 0\), since \(0<x<1\)), and in the limit \(v\) is simply the momentum fraction of the splitting.

Examining the explicit forms of the Altarelli–Parisi splitting kernels, it is clear that integrals involving \(1/z_3\) must be evaluated, and the quadratic expression in the numerator of \(z_3\) appears in the denominator of the integrand, causing the presence of square root functions of \(x\) in the integrated expressions. Because of this, integrating these expressions further, as is necessary when computing the integrals of iterated single unresolved subtraction terms, becomes extremely complicated.

In order to avoid such complications, let us drop all terms that are not linear in \(\alpha\) and in \(v\) in the numerator of \(z_{3}\), so that this numerator simply reads \(\alpha+ x v\). Enforcing the correct collinear limit (\(z_3 \to v\) as \(\alpha\to 0\)) as well as the correspondence between \((v,z_3) \leftrightarrow (1-v,1-z_3)\), we define the new variable⁶ \[\hat{z}_{3} = \frac{\alpha+ x v}{2\alpha+ x}\,. \label{eq:z3hat-def}\tag{20}\] The last requirement ensures the preservation of the symmetry between the daughter partons in the splitting at the integrand level. Thus our choice for the collinear subtraction term is \[{\cal C}_{f_3 f_4}^{(0)}(p_1,p_2,p_3,p_4) \equiv 8\pi\frac{\alpha_{\mathrm{s}}\mu^{2\epsilon}}{C(\epsilon)} \frac{1}{s_{34}} \hat{P}_{f_3 f_4}(\hat{z}_{3}, \hat{k}_{\perp};\epsilon) \otimes |{\cal M}_{Q\bar{Q} g}^{(0)}{(\widehat{\,p\,}_1, \widehat{\,p\,}_2, \widehat{\,p\,}_{34})}|^2 \,, \label{eq:C34-no-F34}\tag{21}\] where the hatted momenta appearing in the factorized matrix element are given in eq. (11 ). We note that \(\hat{z}_{3}\) can be expressed in terms of \(z_3\) of eq. (18 ) as follows, \[\esp{3} &= \frac{2 x (x + y_{34}) z_3 - (x - r)(x + y_{34} - r)}{2 r (2x + y_{34} - r)}\,, \\ 1-\hat{z}_{3} &= \frac{2 x (x + y_{34}) (1-z_3) - (x - r)(x + y_{34} - r)}{2 r (2x + y_{34} - r)}\,, \esp \label{eq:hz3}\tag{22}\] where \[r = \sqrt{(x + y_{34})^2 - 4 y_{34}}\,.\] The definition of the transverse momentum \(\hat{k}_{\perp}\) that enters the Altarelli–Parisi splitting kernel reads \[\hat{k}_{\perp}^\mu = \zeta_{3} p_4^\mu - \zeta_{4} p_3^\mu + \zeta_{34} \widehat{\,p\,}_{34}^\mu\,, \label{eq:hkT}\tag{23}\] where \[\zeta_{3} = z_{3} - \frac{y_{34}}{\alpha y_{(34)P}}\,, \qquad \zeta_{4} = (1-z_{3}) - \frac{y_{34}}{\alpha y_{(34)P}} \qquadand\qquad \zeta_{34} = \frac{y_{34}}{\alpha x} [(1 - z_3) - z_{3}]\,.\] Notice that in the above equation, we have made use of \(z_3\) of eq. (18 ). With this definition \(\hat{k}_{\perp}^\mu\) is perpendicular to the parent momentum \(\widehat{\,p\,}_{34}\) and also \(\hat{k}_{\perp}^\mu \to 0\) in the \(p_3^\mu || p_4^\mu\) collinear limit.

However, even after introducing the new variable in eq. (20 ), the issue of square root functions of invariants in the integrated form of \({\cal C}_{f_3 f_4}^{(0)}\) is still present due to the appearance of the factor \((x-2\alpha+\alpha^2)^{2-2\epsilon}\) in eq. (13 ). We deal with this factor by exploiting the freedom to multiply the subtraction term with a suitable regular function \({\cal F}^C_{34}\), see eq. ([eq:A1bbgg]). To make an optimal choice, we take the occasion to collect all factors coming from the factorized measure and the factor of \(1/s_{34}\) which is common to all collinear integrals. Inserting the explicit expression for \({\mathrm{d}}\phi_{2}(p_3,p_4;\alpha P + \beta\widehat{\,p\,}_{34})\) from eq. (19 ) into eq. (13 ), we find \[[{\mathrm{d}}p]\frac{1}{s_{34}} = \frac{x^{-1+2\epsilon}}{(4\pi)^2} \frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)} (P^2)^{-\epsilon} {\mathrm{d}}v\, \alpha^{-1-\epsilon} v^{-\epsilon} (1-v)^{-\epsilon} (1-\alpha)^{-2+3\epsilon}(x-\alpha)^{-1-\epsilon} (x-2\alpha+\alpha^2)^{2-2 \epsilon}\,.\] We note that the product of the last three factors, \[{\mathcal{G}}(\alpha,x;\epsilon) \equiv (1-\alpha)^{-2+3\epsilon}(x-\alpha)^{-1-\epsilon} (x-2\alpha+\alpha^2)^{2-2 \epsilon}\,,\] does not play a role in regularizing any divergent behaviour, hence the integrand may be simplified (without changing the pole structure of the integral) by multiplying with \[\frac{{\displaystyle \lim_{\alpha\to 0}} {\mathcal{G}}(\alpha,x;\epsilon)}{{\mathcal{G}}(\alpha,x;\epsilon)} = x^{1-3\epsilon}(1-\alpha)^{2-3\epsilon}(x-\alpha)^{1+\epsilon} (x-2\alpha+\alpha^2)^{-2+2 \epsilon}\,.\]

The final source of square roots in the integral is the upper limit of integration in eq. (16 ). Since we are free to restrict the action of the counterterm to a region of phase space around the singular limit, we choose an upper limit \(\alpha_0(x) \le \alpha_{\mathrm{max}}\) such as to avoid the presence of square roots. One simple choice is \[\alpha_0(x) = C\cdot\frac{x}{2} < 1 - \sqrt{1-x}\,,\qquad C,x\in (0,1]\,.\] As the final physical results cannot depend on the constant \(C\), varying its value gives a strong check on the correct implementation of the subtraction scheme. Thus the final form of the regular function \({\cal F}^C_{34}\) is given by \[{\cal F}^C_{34} \equiv x^{1-3\epsilon}(1-\alpha)^{2-3\epsilon}(x-\alpha)^{1+\epsilon}(x-2\alpha+\alpha^2)^{-2+2 \epsilon} \Theta[\alpha_0(x)-\alpha]\,. \label{eq:F34-def}\tag{24}\]

4.1.0.2 Single soft subtraction.

The single soft subtraction to the double real contribution is structurally identical to the NLO soft subtraction term given in eq. (5 ) and we have, \[{\cal S}_{g_r}^{(0)}(p_1,p_2,p_3,p_4) \equiv -8\pi\frac{\alpha_{\mathrm{s}}\mu^{2\epsilon}}{C(\epsilon)} \sum_{\widehat{\,i\,},\widehat{\,k\,}} \frac{1}{2}S_{\widehat{\,i\,}\widehat{\,k\,}}(r)\, {\boldsymbol{T}}_{\widehat{\,i\,}} {\boldsymbol{T}}_{\widehat{\,k\,}} \otimes |{\cal M}_{Q\bar{Q} g}^{(0)}{(\widehat{\,p\,}_1, \widehat{\,p\,}_2, \widehat{\,p\,}_{34})}|^2\,. \label{eq:S3r}\tag{25}\] The mapped momenta that appear in the 3-parton factorized matrix element above can be chosen to coincide with those used to define the collinear subtraction and are given in eq. (11 ). We recall that the summation indices \(\widehat{\,i\,}\) and \(\widehat{\,k\,}\) in eq. (25 ) run over the labels of the mapped momenta that enter the factorized matrix element (i.e., \(\widehat{\,i\,},\widehat{\,k\,} = \widehat{\,1\,}, \widehat{\,2\,}, \widehat{\,34\,}\)).

The integration of the soft counterterm is plagued by similar difficulties as the collinear case discussed above. In particular, the \((x-2\alpha+\alpha^2)^{2-2\epsilon}\) factor in eq. (13 ) is present, as well as the square root in the upper limit of integration. As with the collinear subtraction, we can overcome these problems by a suitable choice of the \({\cal F}^S_r\) function that appears in eq. ([eq:A1bbgg]). In order to obtain this factor, consider the most elaborate soft integral, which involves the eikonal factor \(\frac{s_{\widehat{\,1\,}\widehat{\,2\,}}}{s_{\widehat{\,1\,}r}s_{\widehat{\,2\,}r}}\). It is convenient to write this integral in the rest frame of \(P\), oriented such that \(\widehat{\,p\,}_{34}^\mu\) lies along the \(z\)-axis, \[\esp^\mu &= \sqrt{P^2}(1,\ldots)\,, \\ \widehat{\,p\,}_{34}^\mu &= \hat{E}_{34}(1,\ldots,1)\,, \\ \widehat{\,p\,}_1^\mu &= \hat{E}_1(1,\vec{\beta}_1)\,, \\ \widehat{\,p\,}_2^\mu &= \hat{E}_2(1,\vec{\beta}_2)\,, \esp\] where \(\ldots\) denote components that vanish⁷. In this frame \(p_r\) reads \[p_r^\mu = E_r(1,\vec{n}_r) = E_r(1,\ldots\mathrm{angles}\ldots,\sin\varphi\sin\vartheta_r, \sin\varphi\cos\vartheta_r,\cos\vartheta_r)\,,\] where “\(\ldots\mathrm{angles}\ldots\)” are angular components on which the integrand does not depend. In this frame, \({\mathrm{d}}\phi_{2}(p_3,p_4;\alpha P + \beta\widehat{\,p\,}_{34})\) can be written in the following form \[\esp{2}(p_3,p_4;\alpha P + \beta\widehat{\,p\,}_{34}) &= 2^{-3-2\epsilon}\pi^{-2+\epsilon}\frac{\Gamma(1 - \epsilon)}{\Gamma(1 - 2 \epsilon)} (P^2)^{-\epsilon} {\mathrm{d}}\xi\, {\mathrm{d}}\eta\, \alpha^{1-2\epsilon} \xi^{-\epsilon} (1-\xi)^{-\epsilon} \eta^{-\frac{1}{2}-\epsilon}(1-\eta)^{-\frac{1}{2}-\epsilon} \\ &\times (\alpha+ \beta x)^{1-2\epsilon} (\alpha+\beta x \xi)^{-2+2\epsilon} \Theta(\xi)\Theta(1-\xi)\Theta(\eta)\Theta(1-\eta)\,, \label{eq:dphi2-2} \esp\tag{26}\] where \[\cos\vartheta_r=1-2\xi \qquadand\qquad \cos\varphi=1-2\eta\,.\] Then we find \[\frac{s_{\widehat{\,1\,}\widehat{\,2\,}}}{s_{\widehat{\,1\,}3}s_{\widehat{\,2\,}3}} = \frac{2\widehat{\,p\,}_1\cdot \widehat{\,p\,}_2}{(2\widehat{\,p\,}_1\cdot p_r)(2\widehat{\,p\,}_1\cdot p_r)} = \frac{2\widehat{\,p\,}_1\cdot \widehat{\,p\,}_2}{[2\hat{E}_1 E_r (1-\vec{\beta}_1\cdot\vec{n}_r)] [2\hat{E}_2 E_r (1-\vec{\beta}_2\cdot\vec{n}_r)]}\] while the energy \(E_r\) takes the form \[E_r = \frac{\alpha(\alpha+\beta x)}{\alpha+\beta x \xi} \sqrt{P^2}\,.\] Hence (using \(\hat{E}_1 = y_{\widehat{\,1\,}P} P^2/2\) and \(\hat{E}_2 = y_{\widehat{\,2\,}P} P^2/2\)) \[\esp{\mathrm{d}}p] \frac{s_{\widehat{\,1\,}\widehat{\,2\,}}}{s_{\widehat{\,1\,}3}s_{\widehat{\,2\,}3}} &= 2^{-4 - 2 \epsilon} \pi^{-3 + \epsilon} \frac{\Gamma(1 - \epsilon)}{\Gamma(1 - 2 \epsilon)} (P^2)^{-\epsilon} \frac{y_{\widehat{\,1\,}\widehat{\,2\,}}}{y_{\widehat{\,1\,}P} y_{\widehat{\,2\,}P}} {\mathrm{d}}\xi\, {\mathrm{d}}\eta\, \alpha^{-1 - 2 \epsilon} \frac{\xi^{-\epsilon} (1 - \xi)^{-\epsilon} \eta^{-\frac{1}{2} - \epsilon} (1 - \eta)^{-\frac{1}{2} - \epsilon}}{(1-\vec{\beta}_1\cdot\vec{n}_r) (1-\vec{\beta}_2\cdot\vec{n}_r)} \\ &\times x^{-1 + 2 \epsilon} (1 - \alpha)^{-2 + 2 \epsilon} (x - \alpha)^{-1 - 2 \epsilon} (x -2 \alpha+ \alpha^2)^{2 - 2 \epsilon} [\alpha(1-\alpha) + (x - 2 \alpha+ \alpha^2)\xi]^{2 \epsilon}\,. \esp\] We note that the product of factors on the second line, \[{\mathcal{G}}(\alpha,x,\xi;\epsilon) \equiv x^{-1 + 2 \epsilon} (1 - \alpha)^{-2 + 2 \epsilon} (x - \alpha)^{-1 - 2 \epsilon} (x -2 \alpha+ \alpha^2)^{2 - 2 \epsilon} [\alpha(1-\alpha) + (x - 2 \alpha+ \alpha^2)\xi]^{2 \epsilon}\,,\] does not play a role in regularizing any divergent behaviour, hence the integrand may be simplified (without altering its pole structure) if we multiply it with \[\esp{{\displaystyle \lim_{\alpha\to 0}} {\mathcal{G}}(\alpha,x,\xi;\epsilon)} {{\mathcal{G}}(\alpha,x,\xi;\epsilon)} &= \xi^{2\epsilon} x^{1-2\epsilon} (1-\alpha)^{2-2\epsilon} (x - \alpha)^{1+2\epsilon} (x - 2\alpha+ \alpha^2)^{-2+2\epsilon} \\ &\times [\alpha(1-\alpha) + (x - 2\alpha+ \alpha^2) \xi]^{-2\epsilon}\,. \esp\]

As was the case with the collinear subtraction term, the upper limit of integration again leads to the appearance of square roots in the integral. Following the same strategy as in the case of the collinear subtraction, we arrive at the following formula for \({\cal F}^S_r\)  \[\esp^S_{r} &\equiv \xi^{2\epsilon} x^{1-2\epsilon} (1-\alpha)^{2-2\epsilon} (x - \alpha)^{1+2\epsilon} (x - 2\alpha+ \alpha^2)^{-2+2\epsilon} \\ &\times [\alpha(1-\alpha) + (x - 2\alpha+ \alpha^2) \xi]^{-2\epsilon} \Theta[\alpha_0(x)-\alpha]\,. \esp \label{eq:Fr-def}\tag{27}\] With the above choice of \({\cal F}^S_r\), the soft integral can be performed to yield a fully analytic and reasonably compact expression which is suitable for further integration, as is necessary when computing the integrated forms of the iterated single unresolved counterterms.

4.1.0.3 Single soft-collinear overlap.

The only single unresolved subtraction term in eq. ([eq:A1bbgg]) that we have not yet specified is the soft-collinear overlap \({\cal C}_{g_3 g_4}^{}{\cal S}_{g_r}^{(0)}\). Our choice is _g_3^(0)(p_1,p_2,p_3,p_4)   C_  |M_Q|Q g^(0)(_1, _2, _34)|^2 ,
C_g_3 g_4^S_g_4^(0)(p_1,p_2,p_3,p_4)   C_  |M_Q|Q g^(0)(_1, _2, _34)|^2 . Note that in this subtraction term, the momentum fractions must be evaluated with hatted momenta, so that they match the soft subtraction in the collinear limit. Hence the momentum fractions \(z_{3,\widehat{\,34\,}}\) and \(z_{4,\widehat{\,34\,}}\) are defined as \[z_{3,\widehat{\,34\,}} = \frac{p_3 \cdot P}{(\widehat{\,p\,}_{34} + p_3) \cdot P} \qquadand\qquad z_{4,\widehat{\,34\,}} = \frac{p_4 \cdot P}{(\widehat{\,p\,}_{34} + p_4) \cdot P}\,. \label{eq:z334z434}\tag{28}\] The mapped momenta entering the factorized matrix elements in eqs. ([eq:C34S3]) and ([eq:C34S4]) are once again given by eq. (11 ).

We can now clarify the reason that the soft-collinear overlap terms in eq. ([eq:A1bbgg]) have to be multiplied with the same \({\cal F}^{S}_r\) functions as the soft subtractions. In the \(p_r^\mu \to 0\) soft limit both \({\cal F}^{C}_{34} \to 1\) as well as \({\cal F}^{S}_r \to 1\). Thus \({\cal C}_{g_3 g_4}^{}{\cal S}_{g_r}^{(0)}\) properly regularizes \({\cal C}_{g_3 g_4}^{(0)}\) in the soft limit. On the other hand, in the \(p_3^\mu || p_4^\mu\) collinear limit \({\cal F}^{S}_r \not\to 1\) in \(d\) dimensions. So to insure the proper cancellation of \({\cal S}_{g_r}^{(0)}\) with \({\cal C}_{g_3 g_4}^{}{\cal S}_{g_r}^{(0)}\) in the collinear limit, the latter must be multiplied by the same factor of \({\cal F}^{S}_r\) as the former.

4.1.0.4 Double soft subtraction.

Turning to the double unresolved subtraction, we recall that only the double soft limit requires regularization by subtraction. We choose to define the subtraction term for this limit as follows. For double soft gluon emission we define \[\esp{g_3 g_4}{(0)}(p_1,p_2,p_3,p_4) &\equiv \left[8\pi\frac{\alpha_{\mathrm{s}}\mu^{2\epsilon}}{C(\epsilon)}\right]^2 \bigg\{ \sum_{i,j,k,l=1,2} \frac{1}{8}S_{ik}(3) S_{jl}(4)\, \{{\boldsymbol{T}}_{\widetilde{\,i\,}} {\boldsymbol{T}}_{\widetilde{\,k\,}} , {\boldsymbol{T}}_{\widetilde{\,j\,}} {\boldsymbol{T}}_{\widetilde{\,l\,}}\} \\ &- \frac{1}{4} C_{\mathrm{A}} \sum_{i,k=1,2} \Big[S_{ik}(3,4) - S_{ik}^{\rm mass}(3,4)\Big] {\boldsymbol{T}}_{\widetilde{\,i\,}} {\boldsymbol{T}}_{\widetilde{\,k\,}} \bigg\} \otimes |{\cal M}_{Q\bar{Q} }^{(0)}{(\widetilde{\,p\,}_1,\widetilde{\,p\,}_2)}|^2\,, \label{eq:S3g4g} \esp\tag{29}\] while for a soft quark-antiquark pair we set \[\esp{q_3 {\bar q}_4}{(0)}(p_1,p_2,p_3,p_4) &\equiv \left[8\pi\frac{\alpha_{\mathrm{s}}\mu^{2\epsilon}}{C(\epsilon)}\right]^2 \frac{1}{s_{34}^2} T_{\mathrm{R}} \\ &\times \sum_{i,k=1,2} \frac{s_{i3} s_{k4} + s_{i4} s_{k3} - s_{ik} s_{34}}{(s_{i3}+s_{i4}+s_{34})(s_{k3}+s_{k4}+s_{34})} {\boldsymbol{T}}_{\widetilde{\,i\,}} {\boldsymbol{T}}_{\widetilde{\,k\,}} \otimes |{\cal M}_{Q\bar{Q} }^{(0)}{(\widetilde{\,p\,}_1,\widetilde{\,p\,}_2)}|^2\,. \label{eq:S3q4q} \esp\tag{30}\] The eikonal factors \(S_{ik}(r)\) and \(S_{jl}(r)\) read \[S_{ik}(3) = \frac{2s_{ik}}{s_{i3} s_{k3}} = \frac{(p_i\cdot p_k)}{(p_i\cdot p_3)(p_k\cdot p_3)} \qquadand\qquad S_{jl}(4) = \frac{2s_{jl}}{s_{j4} s_{l4}} = \frac{(p_j\cdot p_l)}{(p_j\cdot p_4)(p_l\cdot p_4)}\,.\] Furthermore, for \(S_{ik}(3,4)\) we have \[S_{ik}(3,4) = S_{ik}^{\rm (s.o.)}(3,4) + 4\frac{s_{i3} s_{k4} + s_{i4} s_{k3}}{(s_{i3} + s_{i4} + s_{34})(s_{k3} + s_{k4} + s_{34})} \bigg[\frac{1-\epsilon}{s_{34}^2} - \frac{1}{8} S_{ik}^{\rm (s.o.)}(3,4) \bigg] - \frac{4}{s_{34}} S_{ik}(34)\,, \label{eq:Sikrs}\tag{31}\] where the \(S_{ik}^{\rm (s.o.)}(3,4)\) is the strongly-ordered limit of this expression in either the \(p_3^\mu\to 0\) or \(p_4^\mu\to 0\) limit (it is symmetric in \(3\) and \(4\)), \[S_{ik}^{\rm (s.o.)}(3,4) = S_{ik}(4)\Big(S_{i4}(3) + S_{k4}(3) - S_{ik}(3)\Big) = \frac{4s_{ik}}{s_{i3} s_{k4} s_{34}} + \frac{4s_{ik}}{s_{i4} s_{k3} s_{34}} - \frac{4s_{ik}^2}{s_{i3} s_{i4} s_{k3} s_{k4}}\] and \[S_{ik}(34) = \frac{2s_{ik}}{(s_{i3} + s_{i4} + s_{34})(s_{k3} + s_{k4} + s_{34})}\,.\] Last, \(S_{ik}^{\rm mass}(3,4)\) is directly proportional to the square of the heavy quark mass, \[\esp_{ik}^{\rm mass}(3,4) &= \frac{s_{i3}}{s_{i3} + s_{i4} + s_{34}} S_{ik}^{\rm mass, (s.o.)}(3,4) + \frac{s_{k4}}{s_{k3} + s_{k4} + s_{34}} S_{ki}^{\rm mass, (s.o.)}(4,3) \\ & -\frac{2}{s_{34}} \frac{s_{i3} s_{k4} + s_{i4} s_{k3}}{(s_{i3} + s_{i4} + s_{34})(s_{k3} + s_{k4} + s_{34})} \bigg[\frac{s_{ii}}{s_{i3} s_{i4}} + \frac{s_{kk}}{s_{k3} s_{k4}}\bigg]\,, \esp\] where \(S_{ik}^{\rm mass, (s.o.)}(3,4)\) and \(S_{ki}^{\rm mass, (s.o.)}(4,3)\) are the strongly-ordered limits of \(S_{ik}^{\rm mass}(3,4)\) in the \(p_4^\mu \to 0\) and \(p_3^\mu \to 0\) limit, _ik^mass, (s.o.)(3,4) &= S_ii(3) (S_k3(4) - S_ik(4)) ,
S_ki^mass, (s.o.)(4,3) &= S_kk(4) (S_i4(3) - S_ik(3)) . Summations over all indices in eq. (29 ) run over \(i,j,k,l=1,2\) and the equivalence of any and all indices is allowed.

We remark that contrary to the choice in eqs. (5 ) and (25 ), here the hard momenta \(p_i\), \(p_k\), \(p_j\) and \(p_l\) that appear in the various functions just defined are simply the original momenta of the heavy quarks in the four-particle phase space and not the mapped momenta. This choice is quite convenient for the present calculation, since it allows us to use known results for massive four particle phase space integrals [58], [59] to compute the integrated subtraction term. For similar reasons, we prefer to define the subtractions in eqs. (29 ) and (30 ) by retaining the subleading (in the double soft limit) \(s_{34}\) term in denominators of the form \((s_{i3}+s_{i4}+s_{34})\) and \((s_{k3}+s_{k4}+s_{34})\) throughout. Thus, our subtraction terms differ by these subleading terms from the double soft limit formulae of [25], [60].

To complete the definition of the subtraction term, we must specify the momenta \(\widetilde{\,p\,}_1\) and \(\widetilde{\,p\,}_2\) that enter the factorized matrix element. Starting from the four momenta of the double real emission phase space, we apply the \(4 \to 3\) mapping of eq. (11 ), followed by the \(3 \to 2\) mapping presented in eq. (7 ) in order to obtain \(\widetilde{\,p\,}_1\) and \(\widetilde{\,p\,}_2\).

Finally, as remarked above, all master integrals that are needed to compute the integrated double soft subtraction term are known in the literature [58], [59], and so we find it most convenient to not include any additional factors with the double soft subtraction, see eqs. ([eq:A2bbgg]) and ([eq:A2bbqq]).

4.1.0.5 Single collinear–double soft subtraction.

In order to cancel the singularities of the double soft subtraction term in the single collinear limit, as well as the singularities of the single collienar subtraction term in the double soft limit, we introduce the iterated single unresolved counterterm \[\esp{f_3 f_4}{}{\cal S}_{f_4 f_4}^{(0)}(p_1,p_2,p_3,p_4) &\equiv \left[8\pi\frac{\alpha_{\mathrm{s}}\mu^{2\epsilon}}{C(\epsilon)}\right]^2 \frac{1}{s_{34}} \sum_{\widehat{\,i\,},\widehat{\,k\,} = \widehat{\,1\,},\widehat{\,2\,}} \frac{2 \widehat{\,p\,}_{i,\mu} \widehat{\,p\,}_{k,\nu}}{s_{\widehat{\,i\,}\widehat{\,34\,}} s_{\widehat{\,k\,}\widehat{\,34\,}}} \\ &\times P_{f_3 f_4}^{\mu\nu}(\hat{z}_{3},\hat{k}_{\perp}{};\epsilon) {\boldsymbol{T}}_{\widetilde{\,i\,}} {\boldsymbol{T}}_{\widetilde{\,k\,}} \otimes |{\cal M}_{Q\bar{Q}}^{(0)}{(\widetilde{\,p\,}_1,\widetilde{\,p\,}_2)}|^2\,, \esp \label{eq:Cf3f4Sf3f4}\tag{32}\] where the Altarelli–Parisi kernels are given in eqs. ([eq:Pgg]) and ([eq:Pqq]), while \(\hat{z}_{3}\) and \(\hat{k}_{\perp}^\mu\) are defined in eqs. (22 ) and (23 ). The momenta \(\widehat{\,p\,}_i\), \(\widehat{\,p\,}_k\) and \(\widehat{\,p\,}_{34}\) that appear in the uncontracted eikonal factor above are those obtained by the \(4\to 3\) mapping of eq. (11 ), while the factorized matrix element is evaluated with the same momenta \(\widetilde{\,p\,}_1\) and \(\widetilde{\,p\,}_2\) that enter the definition of the double soft subtraction term.

We remark that this term enters eq. ([eq:A12bbgg]) multiplied with the factor of \({\cal F}^{C}_{34}\). Since this function goes to one in the collinear limit, \({\cal C}_{q_r {\bar q}_s}^{}{\cal S}_{q_r {\bar q}_s}^{(0)}\) correctly regularizes \({\cal S}_{q_r {\bar q}_s}^{(0)}\) in this limit. On the other hand, in the double soft limit \({\cal F}^{C}_{34} \not\to 1\), so \({\cal C}_{g_r g_s}^{}{\cal S}_{g_r g_s}^{(0)}\) must be multiplied with \({\cal F}^{C}_{34}\) to ensure the proper cancellation of this term with \({\cal C}_{g_r g_s}^{(0)}{\cal F}^{C}_{34}\) in the double soft limit.

4.1.0.6 Single soft–double soft subtraction.

The iterated single soft–double soft subtraction term regularizes the double soft subtraction term in the \(p^\mu_s \to 0\) single soft limit, as well as the single soft subtraction term \({\cal S}_{g_s}^{(0)}\) in the double soft limit, \[\esp{g_s}{}{\cal S}_{g_r g_s}^{(0)}(p_1,p_2,p_3,p_4) &\equiv \left[8\pi\frac{\alpha_{\mathrm{s}}\mu^{2\epsilon}}{C(\epsilon)}\right]^2 \bigg\{ \sum_{\widehat{\,i\,},\widehat{\,k\,},\widehat{\,j\,},\widehat{\,l\,} = \widehat{\,1\,},\widehat{\,2\,}} \frac{1}{8}S_{\widehat{\,i\,}\widehat{\,k\,}}(\widehat{\,34\,}) S_{\widehat{\,j\,}\widehat{\,l\,}}(s)\, \{{\boldsymbol{T}}_{\widetilde{\,i\,}} {\boldsymbol{T}}_{\widetilde{\,k\,}} , {\boldsymbol{T}}_{\widetilde{\,j\,}} {\boldsymbol{T}}_{\widetilde{\,l\,}}\} \\ & - \frac{1}{4} C_{\mathrm{A}}\sum_{\widehat{\,i\,},\widehat{\,k\,} = \widehat{\,1\,},\widehat{\,2\,}} \Big[S_{\widehat{\,i\,}\widehat{\,k\,}}(\widehat{\,34\,})\Big(S_{\widehat{\,i\,}\widehat{\,34\,}}(s) + S_{\widehat{\,k\,}\widehat{\,34\,}}(s) - S_{\widehat{\,i\,}\widehat{\,k\,}}(s)\Big) \\ & - S_{\widehat{\,i\,}\widehat{\,i\,}}(\widehat{\,34\,})\Big(S_{\widehat{\,k\,}\widehat{\,34\,}}(s) - S_{\widehat{\,i\,}\widehat{\,k\,}}(s)\Big) \Big] {\boldsymbol{T}}_{\widetilde{\,i\,}} {\boldsymbol{T}}_{\widetilde{\,k\,}} \bigg\} \otimes |{\cal M}_{Q\bar{Q} }^{(0)}{(\widetilde{\,p\,}_1,\widetilde{\,p\,}_2)}|^2\,. \esp \label{eq:SsgSrgsg}\tag{33}\] As always, summation indices can be equal. For the sake of clarity, we emphasize that \[S_{\widehat{\,i\,}\widehat{\,k\,}}(\widehat{\,34\,}) = \frac{(\widehat{\,p\,}_i\cdot \widehat{\,p\,}_k)}{(\widehat{\,p\,}_i\cdot \widehat{\,p\,}_{34}) (\widehat{\,p\,}_k\cdot \widehat{\,p\,}_{34})} \qquadand\qquad S_{\widehat{\,i\,}\widehat{\,i\,}}(\widehat{\,34\,}) = \frac{\widehat{\,p\,}_i^2}{(\widehat{\,p\,}_i\cdot \widehat{\,p\,}_{34})^2}\,.\] Furthermore we have e.g., \[S_{\widehat{\,j\,}\widehat{\,l\,}}(s) = \frac{(\widehat{\,p\,}_j\cdot \widehat{\,p\,}_l)}{(\widehat{\,p\,}_j\cdot p_s)(\widehat{\,p\,}_l\cdot p_s)} \qquadand\qquad S_{\widehat{\,i\,}\widehat{\,34\,}}(s) = \frac{(\widehat{\,p\,}_i\cdot \widehat{\,p\,}_{34})}{(\widehat{\,p\,}_i\cdot p_s)(\widehat{\,p\,}_{34}\cdot p_s)}\,,\] with obvious generalizations for the other terms not displayed explicitly. Here the set of hatted momenta are obtained from the original momenta via the \(4\to 3\) mapping given in eq. (11 ). The tilded momenta entering the factorized matrix element are again equal to those in the double soft subtraction.

Let us remark that this term enters eq. ([eq:A12bbgg]) multiplied with a factor of \({\cal F}^{S}_s\). Since this function goes to one as \(p_s^\mu \to 0\), \({\cal S}_{g_s}^{}{\cal S}_{g_r g_s}^{(0)}\) regularizes \({\cal S}_{g_r g_s}^{(0)}\) correctly in this limit. On the other hand, in the double soft limit, \({\cal F}^S_{s} \not\to 1\), hence \({\cal S}_{g_s}^{}{\cal S}_{g_r g_s}^{(0)}\) must be multiplied by the same factor as \({\cal S}_{g_s}^{(0)}\) in order to ensure the cancellation of these terms in the double soft limit.

4.1.0.7 Soft-collinear–double soft overlap.

The set of subtractions listed so far leads to double subtraction in the soft-collinear limit. In order to avoid this, we introduce the last counterterm in eq. ([eq:A12bbgg]), given by _g_3^S_g_3 g_4^(0)(p_1,p_2,p_3,p_4) & -^2 _, = , S_()
&   C_  \(T\)_ \(T\)_ |M_Q|Q ^(0)(_1,_2)|^2 ,
C_g_3 g_4^S_g_4^S_g_3 g_4^(0)(p_1,p_2,p_3,p_4) & -^2 _, = , S_()
&   C_  \(T\)_ \(T\)_ |M_Q|Q ^(0)(_1,_2)|^2 . Above \(z_{3,\widehat{\,34\,}}\) and \(z_{4,\widehat{\,34\,}}\) are defined in eq. (28 ). As before, hatted momenta are obtained from the \(4\to 3\) mapping of eq. (11 ), while the tilded momenta are given as in the double soft subtraction term, by applying the \(3 \to 2\) mapping of eq. (7 ).

The factor multiplying \({\cal C}_{g_r g_s}^{}{\cal S}_{g_s}^{}{\cal S}_{g_r g_s}^{(0)}\) in eq. ([eq:A12bbgg]) is doubly constrained. In fact, this term must regularize \({\cal S}_{g_s}^{}{\cal S}_{g_r g_s}^{(0)}{\cal F}^{S}_s\) in the collinear limit as well as \({\cal C}_{g_r g_s}^{}{\cal S}_{g_s}^{(0)}{\cal F}^{S}_{s}\) in the double soft limit. Moreover \({\cal F}^{S}_s\not\to 1\) in either of these limits, so we must multiply \({\cal C}_{g_r g_s}^{}{\cal S}_{g_s}^{}{\cal S}_{g_r g_s}^{(0)}\) by \({\cal F}^{S}_s\) in eq. ([eq:A12bbgg]) to achieve the correct pattern of cancellations.

4.1.0.8 Pattern of cancellations.

Finally, the pattern of cancellations in the various limits among the matrix element and subtraction terms we have introduced is schematically illustrated in figure [fig:RR-line] for the case of double gluon emission. The three directions identify the three singular limits, namely single soft (horizontal arrows), single collinear (diagonal double arrows) and double soft (vertical arrows with double arrowheads). The green, red and blue boxes represent \({\cal A}_1\), \({\cal A}_2\) and \({\cal A}_{12}\)-type subtraction terms. The contour of the boxes reflects the multiplicity of the phase space over which the observable is computed: solid, dashed or dotted for \(J_4\), \(J_3\) and \(J_2\) respectively. Last, the magenta vertical arrow connects the two subtraction terms that are multiplied by \({\cal F}^{C}_{34}\), while the cyan arrows connect counterterms that are multiplied by \({\cal F}^{S}_{r}\). Since for light quark-antiquark pair emission only the single colliear and double soft limits require regularization, in that case only the four leftmost terms in figure [fig:RR-line] are present.

4.2 Regularized real-virtual contribution↩︎

The real-virtual contribution to the differential decay rate involves one-loop corrections to the process \(X(P) \to Q(p_1) + \bar{Q}(p_2) + g(p_3)\) and takes the form \[{\mathrm{d}}\Gamma^{{\rm RV}} = \frac{1}{F} {\mathrm{d}}\phi_{3}(p_1,p_2,p_3;P) 2\Re\langle{\cal M}_{Q\bar{Q} g}^{(0)}|{\cal M}_{Q\bar{Q} g}^{(1)}\rangle\,,\] while the three other terms appearing in eq. ([eq:dGNNLO3]) can be written as _[1] ^RR,A_ &= _3(p_1,p_2,p_3;P) \(I\)_1(p_1,p_2,p_3;) |M_Q|Q g^(0)|^2 ,
^RV,A_ &= _3(p_1,p_2,p_3;P) A_1  ,
(_[1] ^RR,A_)^_1 &= _3(p_1,p_2,p_3;P) A_1 . Let us emphasize that throughout this subsection, \(p_1\), \(p_2\) and \(p_3\) denote the momenta of the heavy quark, heavy antiquark and gluon in the three-particle real emission phase space.

Starting with \(\int_{[1]} {\mathrm{d}}\Gamma^{{\rm RR,A}_{\scriptscriptstyle 1}}\), we note that due to the presence of the factors \({\cal F}^C_{34}\) and \({\cal F}^S_{r}\), the integrated counterterms can be computed straightforwardly by a direct evaluation of their corresponding parametric integral representations. To perform the integrations and manipulate the output, we used the PolyLogTools package of ref. [61]. In the case of the soft subtraction, the angular integrals that appeared were evaluated using the results of ref. [62]. After gathering all contributions, we find that the insertion operator can be written as \[\esp_{1}(p_1,p_2,p_3;\epsilon) = \frac{\alpha_{\mathrm{s}}}{2\pi} \left(\frac{\mu^2}{P^2}\right)^\epsilon \bigg[& C_{\mathrm{A}}{{\mathrm C}}_{g}^{}(y_{3P};\epsilon) +(C_{\mathrm{A}}- 2C_{\mathrm{F}}) {{\mathrm S}}_{mm}^{(1,2)}(y_{1P},y_{2P},w;\epsilon) \\& -C_{\mathrm{A}}\Big( {{\mathrm S}}_{m0}^{(1,3)}(y_{13},y_{3P};\epsilon) +{{\mathrm S}}_{m0}^{(2,3)}(y_{23},y_{3P};\epsilon) \Big) \\& + C_{\mathrm{F}}\Big( {{\mathrm S}}_{m}^{(1,1)}(y_{13},y_{3P};\epsilon) +{{\mathrm S}}_{m}^{(2,2)}(y_{23},y_{3P};\epsilon) \Big) \bigg]\,, \esp \label{eq:I10-RR}\tag{34}\] where the variable \(w\) is defined by \[w = \sqrt{1-\frac{4m_Q^2}{s_{12}+2m_Q^2}}\,.\] In order to write eq. (34 ), we used the colour algebra relations \({\boldsymbol{T}}_1^2={\boldsymbol{T}}_2^2 = C_{\mathrm{F}}\), \({\boldsymbol{T}}_3^2=C_{\mathrm{A}}\) together with \[{\boldsymbol{T}}_1 {\boldsymbol{T}}_2 = \frac{C_{\mathrm{A}}-2C_{\mathrm{F}}}{2} \qquadand\qquad {\boldsymbol{T}}_1 {\boldsymbol{T}}_3 = {\boldsymbol{T}}_2 {\boldsymbol{T}}_3 = -\frac{C_{\mathrm{A}}}{2}\,.\] The various functions that enter eq. (34 ) are as follows. First, the integrated single collinear and single soft-collinear subtraction terms are assembled into the function \({{\mathrm C}}_{g}^{}\), \[{{\mathrm C}}_{g}^{}(x;\epsilon) = \frac{1}{2} [{{\mathrm C}}_{34}^{}]_{gg}(x;\epsilon) +n_l [{{\mathrm C}}_{34}^{}]_{q{\bar q}}(x;\epsilon) -[{{\mathrm C}}_{34}^{}{{\mathrm S}}_{r}^{}](x;\epsilon)\] with _34^]_gg(x;) &= + + - + 12
& - (x) + 4 ^2(x) + 4 Li_2(-2) - 4 _2 + \(\mathrm O(\epsilon^{1})\) ,
_q|q(x;) &= - - - + (x) + \(\mathrm O(\epsilon^{1})\) ,
(x;) &= - + () - 2 ^2() + \(\mathrm O(\epsilon^{1})\) . The single soft subtraction involves a double summation over hard momenta, so the integrated soft counterterm also takes the form of a sum, where the contributions correspond to the integrated eikonal function involving two different massive hard momenta, \({{\mathrm S}}_{mm}^{(1,2)}\), one massive and one massless hard momentum \({{\mathrm S}}_{m0}^{(i,r)}\) and finally a single massive hard momentum, \({{\mathrm S}}_{m}^{(i,i)}\), & S_mm^(1,2)(x_1,x_2,w;) = -() +
& { 4Li_2( x_1 (w(2-x_1-x_2)+x_1-x_2)) -2Li_2( 4 w x_1)
&_2( x_1 (w(2-x_1-x_2)+x_2-x_1)) +2Li_2( 4 w x_1)
& + 2() () +^2(1-w) -^2(1+w)
& +() + (1 ) } + \(\mathrm O(\epsilon^{1})\) ,
& S_m0^(i,r)(y_ir,y_rP;) = + () + ^2() + \(\mathrm O(\epsilon^{1})\) ,
& S_m^(i,i)(y_ir,y_rP;) = 2  S_m0^(i,r)(y_ir,y_rP;) . With these definitions, it is straightforward to show that the sum \[{\mathrm{d}}\Gamma^{{\rm RV}} + \int_{[1]} {\mathrm{d}}\Gamma^{{\rm RR,A}_{\scriptscriptstyle 1}} \label{eq:RV43RRA1}\tag{35}\] is free of explicit poles in \(\epsilon\). However, both terms develop non-integrable singular behaviour as the gluon becomes soft. We deal with these divergences by introducing appropriate subtraction terms. Since only the \(p^\mu_3\to 0\) soft gluon limit requires regularization, the structure of the approximate matrix elements are very simple, _1 &= S_g_3^(1) ,
A_[1]&= S_g_3^(I) .

4.2.0.1 Real-virtual single soft subtraction.

Starting with the real-virtual contribution, we consider the general expression for the soft current at one-loop for massive amplitudes computed in [63], [64]. Collecting terms in this formula that do not automatically vanish by colour conservation and using the \(3 \to 2\) momentum mapping \(\{p_1,p_2,p_3\} \to \{\widehat{\,p\,}_1,\widehat{\,p\,}_2\}\) of eq. (7 ), our choice for the counterterm is given by \[\esp{g_r}{(1)}(p_1,p_2,p_3) &\equiv -8\pi\frac{\alpha_{\mathrm{s}}\mu^{2\epsilon}}{C(\epsilon)} \bigg\{ \sum_{\widehat{\,i\,},\widehat{\,k\,}=\widehat{\,1\,},\widehat{\,2\,}} \frac{1}{2}S_{\widehat{\,i\,}\widehat{\,k\,}}(r)\, {\boldsymbol{T}}_{\widehat{\,i\,}} {\boldsymbol{T}}_{\widehat{\,k\,}} \otimes 2\Re\langle{\cal M}_{Q\bar{Q}}^{(0)}(\widehat{\,p\,}_1,\widehat{\,p\,}_2)|{\cal M}_{Q\bar{Q}}^{(1)}(\widehat{\,p\,}_1,\widehat{\,p\,}_2)\rangle \\& +2C_{\mathrm{A}} \sum_{\substack{\widehat{\,i\,},\widehat{\,k\,}=\widehat{\,1\,},\widehat{\,2\,} \\ \widehat{\,i\,} \ne \widehat{\,k\,}}} \left[\frac{1}{2}S_{\widehat{\,i\,}\widehat{\,k\,}}(r) - \frac{1}{2}S_{\widehat{\,i\,}\widehat{\,i\,}}(r)\right] R_{\widehat{\,i\,}\widehat{\,k\,}} {\boldsymbol{T}}_{\widehat{\,i\,}} {\boldsymbol{T}}_{\widehat{\,k\,}} \otimes |{\cal M}_{Q\bar{Q}}^{(0)}{(\widehat{\,p\,}_1,\widehat{\,p\,}_2)}|^2 \\& - \frac{\alpha_{\mathrm{s}}}{2 \pi} \frac{1}{C(\epsilon)} \frac{1}{2 \epsilon} \bigg[\beta_0 + \frac{4}{3} T_{\mathrm{R}}\bigg(\frac{\mu_R^2}{m_Q^2}\bigg)^\epsilon\bigg] \sum_{\widehat{\,i\,},\widehat{\,k\,}=\widehat{\,1\,},\widehat{\,2\,}} \frac{1}{2}S_{\widehat{\,i\,}\widehat{\,k\,}}(r)\, {\boldsymbol{T}}_{\widehat{\,i\,}} {\boldsymbol{T}}_{\widehat{\,k\,}} \otimes |{\cal M}_{Q\bar{Q}}^{(0)}{(\widehat{\,p\,}_1,\widehat{\,p\,}_2)}|^2 \bigg\} \,, \esp \label{eq:S3g-1-loop}\tag{36}\] where the definition of the eikonal factor in eq. (6 ) is recalled here for convenience, \[S_{\widehat{\,i\,}\widehat{\,k\,}}(r) = \frac{2 s_{\widehat{\,i\,}\widehat{\,k\,}}}{s_{\widehat{\,i\,}r}s_{\widehat{\,k\,}r}} = \frac{(\widehat{\,p\,}_i\cdot \widehat{\,p\,}_k)}{(\widehat{\,p\,}_i\cdot p_r)(\widehat{\,p\,}_k\cdot p_r)} \qquadand\qquad S_{\widehat{\,i\,}\widehat{\,i\,}}(r) = \frac{2 s_{\widehat{\,i\,}\widehat{\,i\,}}}{s_{\widehat{\,i\,}r}^2} = \frac{\widehat{\,p\,}_i^2}{(\widehat{\,p\,}_i\cdot p_r)^2}\,.\] Note that the contribution on the last line of eq. (36 ) contains the terms that arise form the renormalization of the one-loop soft current. The one-loop function \(R_{\widehat{\,i\,}\widehat{\,k\,}}\) can be written in the following form, \[R_{\widehat{\,i\,}\widehat{\,k\,}} = \frac{\alpha_{\mathrm{s}}}{2\pi} \left(\frac{1}{2}S_{\widehat{\,i\,}\widehat{\,k\,}}(r) \mu^2\right)^\epsilon \left[-\frac{1}{2\epsilon^2} -\frac{1}{2}\sum_{n=-1}^{1} \epsilon^n R_{\widehat{\,i\,}\widehat{\,k\,}}^{(n)}\right]\,, \label{eq:R95ik-hat}\tag{37}\] where we have adapted the prefactor to our conventions. The functional forms of the \(R_{\widehat{\,i\,}\widehat{\,k\,}}^{(n)}\) coefficients are taken from [64], _^(-1) &= ^(_+) - (^() + ^())  ,
R_^(0) &= + ^2(_+) + _2 ( - )  ,
R_^(1) &= { (1 - (d_+d_)) - 7 _2 ( ^() + ^() )
& + } + { ( ^() + ^() )
& - ( ^3() + ^3() ) + 2 Li_3(1-) + Li_3() + 12 _2 ()
& - ^(_+) } + ^3(_+) - ( + ) _3  . Due to the different choice of prefactors, the forms given above are not exactly equal to those in ref. [64]. In particular, \(R_{\widehat{\,i\,}\widehat{\,k\,}}^{(1)}\) differs from the expression presented in [64] by terms proportional to \(\zeta_{2}R_{\widehat{\,i\,}\widehat{\,k\,}}^{(-1)}\). However, we note that the order \(ep\) coefficient \(R_{\widehat{\,i\,}\widehat{\,k\,}}^{(1)}\) is only relevant to compute the integrated subtraction term, but otherwise does not enter the regularized real-virtual contribution that is actually integrated numerically in four dimensions. The variables in the equations above are defined as [64] \[\begin{gather} \alpha_{\widehat{\,i\,}}\equiv \frac{m_i^2 \, s_{\widehat{\,k\,}r}}{s_{\widehat{\,i\,}\widehat{\,k\,}} s_{\widehat{\,i\,}r}} \,, \quad \alpha_{\widehat{\,k\,}}\equiv \frac{m_k^2 \, s_{\widehat{\,i\,}r}}{s_{\widehat{\,i\,}\widehat{\,k\,}} s_{\widehat{\,k\,}r}} \,, \quad d_{\widehat{\,i\,}}\equiv 1 - 2 \alpha_{\widehat{\,i\,}}\,, \quad d_{\widehat{\,k\,}}\equiv 1 - 2 \alpha_{\widehat{\,k\,}}\,, \\ \widehat{\,v\,} \equiv \sqrt{1 - 4 \alpha_{\widehat{\,i\,}}\alpha_{\widehat{\,k\,}}} \; , \quad \widehat{\,v\,}_{\pm} \equiv \frac{1 \pm \widehat{\,v\,}}{2} \,, \quad \widehat{\,x\,} \equiv \frac{\widehat{\,v\,}_-}{\widehat{\,v\,}_+} \,. \end{gather} \label{eq:1-loop-vars-hat}\tag{38}\] Notice, that similarly to the tree-level single soft subtraction term in eq. (5 ), the eikonal factors and variables are computed using the hard momenta which appear in the factorized matrix elements in eq. (36 ). In our specific case, this leads to simplifications, since \(\widehat{\,v\,}\) reduces to a function of just the (fixed) heavy quark mass and \(P^2\). In particular, we find \[\widehat{\,v\,} = \frac{1 - y^2}{1 + y^2}\] in terms of the variable \(y\) of eq. (10 ).

4.2.0.2 Single soft subtraction to the integrated single unresolved counterterm.

Finally, the subtraction for the integrated single unresolved counterterm, eq. ([eq:A1-IxM0]) can be easily defined for our purpose as follows: \[{\cal S}_{g}^{(I\otimes 0)}(p_1,p_2,p_3) = -8\pi\frac{\alpha_{\mathrm{s}}\mu^{2\epsilon}}{C(\epsilon)} \sum_{\widehat{\,i\,},\widehat{\,k\,}=\widehat{\,1\,},\widehat{\,2\,}} \frac{1}{2}S_{\widehat{\,i\,}\widehat{\,k\,}}(r)\, {\boldsymbol{T}}_{\widehat{\,i\,}} {\boldsymbol{T}}_{\widehat{\,k\,}} \otimes {\boldsymbol{I}}_{1,S}(p_1,p_2,p_3;\epsilon) |{\cal M}_{Q\bar{Q}}^{(0)}{(\widehat{\,p\,}_1,\widehat{\,p\,}_2)}|^2\,. \label{eq:S3g-Ix0}\tag{39}\] Similarly to the real-virtual single soft subtraction term, the momenta entering the factorized matrix element in eq. (39 ) above are obtained from the \(3\to 2\) momentum mapping of eq. (7 ). In the soft limit, the insertion operator reads \[\esp_{1,S}(p_1,p_2,p_3;\epsilon) = \frac{\alpha_{\mathrm{s}}}{2\pi} \left(\frac{\mu^2}{P^2}\right)^\epsilon \bigg[& C_{\mathrm{A}}{{\mathrm C}}_{g}^{}(y_{3P};\epsilon) +(C_{\mathrm{A}}- 2C_{\mathrm{F}}) {{\mathrm S}}_{mm}^{(1,2)}(y_{1P},y_{2P},w_S;\epsilon) \\& -C_{\mathrm{A}}\Big( {{\mathrm S}}_{m0}^{(1,3)}(y_{13},y_{3P};\epsilon) +{{\mathrm S}}_{m0}^{(2,3)}(y_{23},y_{3P};\epsilon) \Big) \\& + C_{\mathrm{F}}\Big( {{\mathrm S}}_{m}^{(1,1)}(y_{13},y_{3P};\epsilon) +{{\mathrm S}}_{m}^{(2,2)}(y_{23},y_{3P};\epsilon) \Big) \bigg]\,, \esp\] and the only difference between \({\boldsymbol{I}}_{1,S}\) and \({\boldsymbol{I}}_{1}\) is that \({{\mathrm S}}_{mm}^{(1,2)}\) must be evaluated with the variable \(w\) computed in the soft limit, that we denote as \(w_S\). Expressing \(w_S\) with the variable \(y\) of eq. (10 ), we find \[w_S = \frac{1 - y}{1 + y}\,.\]

4.2.0.3 Pattern of cancellations.

To finish this subsection, we illustrate the structure of cancellations among the various terms in figure 2. The soft limit of the real-virtual matrix element is regulated in \(d\) dimensions by the \({\cal S}_{3}^{(1)}\) subtraction defined in eq. (36 ). We note that the difference \(\left[2\Re\langle{\cal M}_{Q\bar{Q} g}^{(0)}|{\cal M}_{Q\bar{Q} g}^{(1)}\rangle - {\cal S}_{3}^{(1)}\right]\), although free of non-integrable kinematical singularities is not finite in \(\epsilon\). Similarly, \({\boldsymbol{I}}_{1} \otimes |{\cal M}_{Q\bar{Q} g}^{(0)}{}|^2\) is regularized in the soft limit by the subtraction term \({\cal S}_{g}^{(I\otimes 0)}\) in eq. (39 ), but the difference \(\left[{\boldsymbol{I}}_{1} \otimes |{\cal M}_{Q\bar{Q} g}^{(0)}{}|^2 - {\cal S}_{g}^{(I\otimes 0)}\right]\) still contains poles in \(\epsilon\). However, since the \(\epsilon\) poles of \(2\Re\langle{\cal M}_{Q\bar{Q} g}^{(0)}|{\cal M}_{Q\bar{Q} g}^{(1)}\rangle\) and \({\boldsymbol{I}}_{1} \otimes |{\cal M}_{Q\bar{Q} g}^{(0)}{}|^2\) explicitly cancel (see eq. (35 )), we must have that the explicit poles in \(\left[{\cal S}_{3}^{(1)} + {\cal S}_{g}^{(I\otimes 0)}\right]\) also cancel. This can be easily verified directly, using the explicit expressions presented above. In figure 2, the cancellation of explicit \(\epsilon\)-poles is represented by vertical arrows with full arrowheads, while the regularization of kinematic singularities in the single soft limit is denoted by horizontal arrows. The contour of the boxes again reflects the multiplicity of the phase space over which the observable is computed: dashed or dotted for \(J_3\) and \(J_2\) respectively.

4.3 Regularized double virtual contribution↩︎

Finally, the regularized double virtual contribution is the sum of the two-loop corrections to the \(X(P)\to Q(p_1) + \bar{Q}(p_2)\) process and the four integrated counterterms that we have not yet discussed, see eq. ([eq:dGNNLO2]). Note that throughout this subsection, \(p_1\) and \(p_2\) denote the momenta of the heavy quark \(Q\) and the heavy antiquark \(\bar{Q}\).

In order to integrate the remaining subtraction terms, we follow a dual strategy. First, the double soft subtraction terms can be reduced to known four-particle massive phase space integrals [58], [59] via integration-by-parts (IBP) identities. The IBP reduction is rather straightforward and so the integrated double soft subtraction can be obtained easily. As for the rest of the necessary integrated subtraction terms, we performed a direct evaluation of their various parametric representations, similarly to the case of the single unresolved subtraction terms discussed in 4.2.

The collection of all integrated counterterms in eq. ([eq:dGNNLO2]) can be written in the following form \[\bsp & \int_2 {\mathrm{d}}\Gamma^{{\rm RR,A}_{\scriptscriptstyle 2}} - \int_2 {\mathrm{d}}\Gamma^{{\rm RR,A}_{\scriptscriptstyle 12}} + \int_1 {\mathrm{d}}\Gamma^{{\rm RV,A}_{\scriptscriptstyle 1}} + \int_1 \left(\int_1 {\mathrm{d}}\Gamma^{{\rm RR,A}_{\scriptscriptstyle 1}}\right)^{\!{\rm A}_1} = \\ &\qquad\qquad= \bigg\{{\boldsymbol{I}}_{2}(p_1,p_2;\epsilon) - \frac{\alpha_{\mathrm{s}}}{2 \pi} \frac{1}{C(\epsilon)}\frac{1}{2 \epsilon} \bigg[\beta_0 + \frac{4}{3} T_{\mathrm{R}}\bigg(\frac{\mu_R^2}{m_Q^2}\bigg)^\epsilon\bigg] {\boldsymbol{I}}_{1}(p_1,p_2;\epsilon) \bigg\} \otimes |{\cal M}_{Q\bar{Q}}^{(0)}{}|^2 \\&\qquad\qquad + {\boldsymbol{I}}_{1}(p_1,p_2;\epsilon) \otimes 2\Re\langle{\cal M}_{Q\bar{Q}}^{(0)}|{\cal M}_{Q\bar{Q}}^{(1)}\rangle\,. \esp \label{eq:IVV-sum}\tag{40}\] Notice that the term proportional to \({\boldsymbol{I}}_{1}(p_1,p_2;\epsilon) \otimes |{\cal M}_{Q\bar{Q}}^{(0)}{}|^2\) on the second line corresponds precisely to the renormalization counteterm of the one-loop soft current in eq. (36 ). We find it convenient to keep this term explicit for an easy conversion to the case of multiple heavy quarks.

The \({\boldsymbol{I}}_{1}(p_1,p_2;\epsilon)\) insertion operator has been given explicitly in eq. (9 ), while \({\boldsymbol{I}}_{2}(p_1,p_2;\epsilon)\) takes the following form, _2(p_1,p_2;) = ^2 {& C_^2 ( + + b_0) + C_C_( + + c_0)
&+ C_T_n_l ( + + d_0) + \(\mathrm O(\epsilon^{1})\) } , where the coefficients of the Laurent expansion are functions of the \(y\) variable given in eq. (10 ). The \(b\), \(c\) and \(d\) coefficients are as follows (as before we use compact notation \(G_{a_1,\ldots,a_n} = G_{a_1,\ldots,a_n}(y)\)), _-2 &= 2 +G_0 +G_0,0
b_-1 &= 8 -16 G_1 +G_0 +(2 G_-1,0-4 G_0,1-G_1,0-2 _2) -
&
b_0 &= -156 -64 (G_1+) +128 G_1,1 -G_0,1 -
& -G_1,0,0 -G_0,1,0 -
&
c_-2 &= +G_0
c_-1 &= - G_1 +
& +
c_0 &= -G_1 + - 24 +G_1,1 + +
d_-2 &= - -G_0
d_-1 &= - + G_1 + G_0 - (4 G_-1,0-G_0,0-4 G_0,1+2 G_1,0-4 _2)
d_0 &= - + + G_1 - G_1,1 + +G_0,0 . Similarly to \({\boldsymbol{I}}_{1}\), in eq. ([eq:I20]) we have expanded a factor of \(1/C(\epsilon)^2\) coming from the definition of the renormalized coupling, thereby canceling terms of \(\gamma_E\) and \(\ln(4\pi)\). Switching to the standard \(\overline{\mathrm{MS}}\) convention would imply \[b_0^{\overline{\mathrm{MS}}} = b_0 + \zeta_{2} b_{-2}\,, \qquad c_0^{\overline{\mathrm{MS}}} = c_0 + \zeta_{2} c_{-2} \qquadand\qquad d_0^{\overline{\mathrm{MS}}} = d_0 + \zeta_{2} d_{-2}\,,\] with the rest of the expansion coefficients unchanged.

Let us make two comments about the \({\boldsymbol{I}}_2\) operator. First, \({\boldsymbol{I}}_2\) corresponds to the sum of all integrated counterterms which involve the tree-level factorized matrix element \(|{\cal M}_{Q\bar{Q}}^{(0)}{}|^2\), except for the renormalization term of the one-loop soft current which, as noted above, we keep explicit. Second, although the Laurent expansion of \({\boldsymbol{I}}_2\) starts at \(1/\epsilon^2\), individual contributions to this operator involve qubic poles. These poles come from on the one hand double real configurations where the two gluons are collinear and both soft. Similarly, the double poles of the real-virtual contribution are proportional to the tree level three-parton matrix element that develops an extra pole upon integration over the soft region of phase space. Given that the double virtual matrix element is free of triple poles, these must cancel upon combining all integrated subtraction terms.

The above formulae, together with the subtraction terms given in the previous sections, complete the full set of ingredients of our subtraction scheme. The implementation of the entire procedure in a numeric code is then straightforward.

5 Example: Higgs boson decay to massive bottom quarks↩︎

As stated in the Introduction, the construction given above can be applied to compute fully differential NNLO QCD corrections to any process with a colourless initial state decaying into a massive quark-antiquark pair at leading order. As an illustrative example, in this section we report on such a computation for a Standard Model Higgs boson decaying into a pair of massive bottom quarks.

We recall that the necessary two-loop currents have been computed in ref. [41]–[44] and ref. [45]. We have verified the correctness of our implementations of these formulae by checking the exact agreement among them, after accounting for the different conventions. The three parton one-loop and four parton tree level matrix elements were obtained with a straightforward, direct Feynman-diagram calculation and cross-checked with GoSam [65], [66].

To validate our construction, we start by examining the total decay width, corresponding to \(J=1\), evaluated in the on-shell renormalization scheme for the heavy quark, \[\Gamma^{\mathrm{}}_{b{\bar b}} = \Gamma^{\mathrm{LO}}_{b{\bar b}} \left[1 + \frac{\alpha_{\mathrm{s}}}{\pi} \gamma^{(1)}_{b{\bar b}} + \left(\frac{\alpha_{\mathrm{s}}}{\pi}\right)^2 \gamma^{(2)}_{b{\bar b}} + \mathrm O(\alpha_{\mathrm{s}}^{3}) \right]\,.\] The NLO correction given by \(\gamma^{(1)}_{b{\bar b}}\) has been known analytically for a long time [67], [68], while \(\gamma^{(2)}_{b{\bar b}}\) has been computed as a series expansion in \((m_b^2/m_H^2)\) up to the fourth power [69]. It has also been calculated exactly for physical values of the bottom quark and Higgs boson masses recently [39], [40] and we find perfect agreement with these results.

In order to investigate the validity of the approximate formula for the NNLO correction to the total decay width for values of heavy quark masses approaching the kinematic threshold, in figure 3 we compare it to the exact computation.

The upper panel shows the value of \(\gamma^{(2)}_{b{\bar b}}\) as a function of the heavy quark mass \(m_b\), with the Higgs boson mass fixed to its physical value of \(m_H=125.09\,GeV\). In order to better appreciate the level of agreement, in the lower panels we present the ratio of the exact result to the approximate one. We observe that up to around \(38\,GeV\) (near the threshold for the production of four heavy quarks), the agreement is well within 1%. The reason for the discontinuity observed in the ratio for \(38\,GeV < m_b < 39\,GeV\) is simply due to the fact that the exact and approximate results vanish for slightly different values of the heavy quark mass. Between \(40\,GeV\) and \(46\,GeV\), the difference between the two results is still below 1%. For larger values of the heavy quark mass approaching the threshold, an all-order expansion in \((m_b^2/m_H^2)\) would be needed, that is indeed provided by the exact result.

We illustrate the computation of a differential quantity by clustering the partons in the final state into jets with the Durham algorithm [70] with the resolution variable set to \(y_{\mathrm{cut}}=0.1\). In figure 4, we present the differential decay rate in the \(\overline{\mathrm{MS}}\) scheme with respect to the rapidity of the most energetic jet. As can be seen on the figure, this variable has a non-singular distribution already at LO and so genuine NNLO corrections contribute bin by bin. The results presented in figure 4 were obtained with the following setup. The Higgs boson mass was set to \(m_H = 125.09\,GeV\), while the on-shell bottom quark mass was \(m_b=4.78\,GeV\), which corresponds to \(\overline{m}_b(m_H) \simeq 2.79\,GeV\) in the \(\overline{\mathrm{MS}}\) scheme using two-loop running. The strong coupling at the relevant renormalization scale was evolved using three-loop running starting form \(\alpha_{\mathrm{s}}(M_Z) = 0.118\).

We recall that the relation between results computed in the on-shell and \(\overline{\mathrm{MS}}\) schemes (denoted with a bar) is given as follows, \[\overline{\Gamma}^{\mathrm{}}_{b{\bar b}}[J] = \overline{\Gamma}^{\mathrm{LO}}_{b{\bar b}}[J] + \overline{\Gamma}^{\mathrm{NLO}}_{b{\bar b}}[J] + \overline{\Gamma}^{\mathrm{NNLO}}_{b{\bar b}}[J] + \mathrm O(\alpha_{\mathrm{s}}^{3})\,,\] where _b|b[J] &= ^_b|b[J] ,
^_b|b[J] &= {^_b|b[J] + r_1 ^_b|b[J]} ,
^_b|b[J] &= {^_b|b[J] + r_1 ^_b|b[J] + r_2 ()^2 ^_b|b[J] } . Furthermore [71]–[74] \[r_1 = -2d_1 \qquadand\qquad r_2 = 3d_1^2 - 2d_2\] with _1 &= -C_(1 + L) ,
d_2 &= C_^2 ( - _3 + 3 _2 - _2 + L + L^2)
& + C_C_(- + _3 + _2 - _2 - L - L^2)
& + C_T_n_l ( + _2 + L + L^2) + C_T_( - _2 + L + L^2) , where \(L = \ln(\mu_R^2/m_b^2)\). The relation between the Yukawa couplings in the two schemes is given by \[y_b^2 = \overline{y}_b^2 \bigg[1 + r_1 \frac{\alpha_{\mathrm{s}}(\mu_R)}{\pi} + r_2 \left(\frac{\alpha_{\mathrm{s}}(\mu_R)}{\pi}\right)^2 + \mathrm O(\alpha_{\mathrm{s}}^{3})\bigg]\,.\] Although in the \(\overline{\mathrm{MS}}\) scheme the running mass at the scale around the Higgs boson mass is significantly reduced with respect to the on-shell value, as is customarily done, we prefer to keep the on-shell mass in the definition of the kinematics for the outgoing heavy quark momenta in order to mimic the effects related to hadronization that will produce mesons with masses close to that value.

To exhibit the reduced theoretical uncertainty due to the higher order contributions, we vary the renormalization scale around \(m_H\) by a factor of two in both directions. With the inclusion of the NNLO corrections, we observe a nice convergence of the perturbative expansion and the corresponding reduction in the leftover theoretical uncertainty parametrized by scale variation.

6 Conclusions↩︎

In this paper, we have presented a completely local subtraction scheme for computing fully differential NNLO corrections to the production of a heavy quark-antiquark pair from a colourless initial state. Following the CoLoRFulNNLO method, the construction of our subtraction terms starts from the known singular limits of tree-level and one-loop massive matrix elements supplemented by momentum mappings that enforce exact phase space factorization. However, we furthermore employ a global strategy simplifying simultaneously the computation of the integrated counterterms for single and iterated single unresolved emission. This strategy is quite general and can be applied in principle also for more generic processes. We have implemented further simplifications for the specific case of heavy quark-antiquark pair production by including certain subleading contributions to the general formulae in the double soft limit and the single soft limit for the one-loop heavy quark current. As a result, we were able to obtain a very compact analytic result for the sum of all integrated subtraction terms with a number of terms comparable to that of the two-loop virtual amplitude.

Finally, we have shown the application of our method for the case of a Standard Model Higgs boson decaying to a heavy quark-antiquark pair. First, we have compared our results for the NNLO correction to the inclusive decay rate to the approximate formula of [69], based on a series expansion in \((m_b^2/m_H^2)\). Varying the heavy quark mass, we find excellent agreement up to values of the heavy quark mass where higher order effects in the mass expansion can no longer be neglected. Furthermore, as an illustrative example of a differential calculation, we have presented the leading jet rapidity distribution at NNLO accuracy.

We conclude by remarking that the present paper contains in full detail all formulae that are needed to reproduce the results discussed above, and to extend the computation to other heavy quark-antiquark pair production processes from colourless initial states.

Acknowledgments↩︎

We thank Roberto Bonciani, Luca Buonocore, Vittorio Del Duca and Pierpaolo Mastrolia for helpful discussions. Also we thank Werner Bernreuther, Christian Bogner and Oliver Dekkers for kindly providing us the expressions for the master integrals computed in [59]. This work was supported by grant K 125105 of the National Research, Development and Innovation Fund in Hungary. The research of FT is supported by INFN.

References↩︎