\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${t {{\bar{t}}}H}$$\end{document}tt¯H production at NNLO: the flavour off-diagonal channels

We consider QCD radiative corrections to the associated production of a heavy-quark pair (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q{{\bar{Q}}}$$\end{document}QQ¯) with a generic colourless system F at hadron colliders. We discuss the resummation formalism for the production of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q{{\bar{Q}}}F$$\end{document}QQ¯F system at small values of its total transverse momentum \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_T$$\end{document}qT. We present the results of the corresponding resummation coefficients at next-to-leading and, partly, next-to-next-to-leading order. The perturbative expansion of the resummation formula leads to the explicit ingredients that can be used to apply the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_T$$\end{document}qT subtraction formalism to fixed-order calculations for this class of processes. We use the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_T$$\end{document}qT subtraction formalism to perform a fully differential perturbative computation for the production of a top-antitop quark pair and a Higgs boson. At next-to-leading order we compare our results with those obtained with established subtraction methods and we find complete agreement. We present, for the first time, the results for the flavour off-diagonal partonic channels at the next-to-next-to-leading order.


Introduction
The observation of Higgs boson production in association with a top quark-antiquark (tt) pair was reported by the ATLAS and CMS collaborations in 2018 [1,2]. This production mode allows for a direct measurement of the top-quark Yukawa coupling.
The first theoretical studies of the production of a topantitop quark pair and a Higgs boson (tt H) were carried out in Refs. [3,4] at leading order (LO) in QCD perturbation theory, and in Refs. [5][6][7][8][9][10] at next-to-leading order (NLO). NLO EW corrections were reported in Refs. [11][12][13]. The resuma e-mail: stefan.kallweit@cern.ch (corresponding author) mation of soft-gluon contributions close to the partonic kinematical threshold was considered in Refs. [14][15][16][17][18][19]. Full offshell calculations with decaying top quarks were presented at NLO QCD [20] and NLO QCD+EW [21]. The uncertainty of current theoretical predictions for tt H cross sections is at the O(10%) level [22]. To match the experimental precision expected at the end of the high-luminosity phase of the LHC, next-to-next-to-leading order (NNLO) predictions in QCD perturbation theory are required. This paper is devoted to the NNLO (and NLO) QCD calculation of tt H production. At the partonic level, the NNLO calculation of tt H production requires the evaluation of treelevel contributions with two additional unresolved partons in the final state, of one-loop contributions with one unresolved parton and of purely virtual contributions. The required treelevel and one-loop scattering amplitudes can nowadays be evaluated with automated tools. The two-loop amplitude for tt H production is not known. Being a five-leg amplitude involving particles with different masses, its computation is at the frontier of current possibilities [23].
Even having all the required amplitudes, their implementation in a complete NNLO calculation at the fully differential (exclusive) level is a highly non-trivial task because of the presence of infrared (IR) divergences at intermediate stages of the calculation. In particular, these divergences do not permit a straightforward implementation of numerical techniques. Various methods have been proposed and used to overcome these difficulties at the NNLO level (see Refs. [23][24][25][26] and references therein).
In this work we will use the transverse-momentum (q T ) subtraction method [27]. The q T subtraction formalism is a method to handle and cancel the IR divergences in QCD computations at NLO, NNLO and beyond. The method uses IR subtraction counterterms that are constructed by consid-ering and explicitly computing the q T distribution of the produced final-state system in the limit q T → 0. If the produced final-state system is composed of non-QCD (colourless) partons (such as vector bosons, Higgs bosons, and so forth), the behaviour of the q T distribution in the limit q T → 0 has a universal (process-independent) structure that is explicitly known up to the NNLO level through the formalism of transverse-momentum resummation [28]. These results on transverse-momentum resummation are sufficient to fully specify the q T subtraction formalism up to NNLO for this entire class of processes. 1 The resummation formalism can, however, be extended to the production of final states containing a heavy-quark pair [36][37][38]. Exploiting such extension, the NNLO computations of top-quark and bottom-quark production were recently completed [39][40][41][42].
In this paper we consider the associated production of a top-quark pair with a Higgs boson. Since the Higgs boson is colourless, the structure of transverse-momentum resummation for this process is closely analogous to that for heavyquark production. The only important difference is that the emission of a Higgs boson off the top-quark pair changes the kinematics, and the top and antitop quarks are not back-toback anymore at Born level. We show that this feature can be controlled through the knowledge of appropriate resummation coefficients. We present the explicit results of the resummation coefficients at NLO and, partly, NNLO for the process of associated production of an arbitrary number of heavy quark-antiquark pairs and a generic colourless system F [43]. This allows us to obtain first results on the application of the q T subtraction method to the NLO and NNLO computations of tt H production in hadron collisions. We exploit the formulation of transverse-momentum resummation in Ref. [38] that includes the complete dependence on the kinematics of the heavy-quark pair. This dependence and, in particular, the complete control on the heavy-quark azimuthal correlations are essential (see Sect. 2) to extract all the NNLO counterterms of the q T subtraction method. Although the structure of transverse-momentum resummation for tt H production is fully worked out up to NNLO, the explicit NNLO results for the hard-virtual factors [38] in the flavour diagonal partonic channels qq → tt H + X and gg → tt H + X (X denotes the unobserved inclusive final state) are not yet known. Their evaluation requires the two-loop amplitudes for the partonic processes qq → tt H and gg → tt H, as well as related soft contributions whose computation is available only for heavy-quark pair production. Therefore, in the NNLO calculation of this paper we present numerical results for all the flavour off-diagonal channels ab → tt H + X , with 1 For this class of processes transverse-momentum resummation has recently been extended [29][30][31][32][33] to next-to-next-to-next-to leading order (N 3 LO), thus allowing first N 3 LO applications [34,35] of the q T subtraction method. ab = qg(qg), qq(qq), qq (qq ), qq (qq ) (q and q denote quarks with different flavours).
The paper is organised as follows. In Sect. 2 we recall the transverse-momentum resummation formalism for the production of a high-mass system containing a heavy-quark pair and discuss the perturbative ingredients needed for the NNLO calculation. In Sect. 3 we present our numerical results for tt H production at NLO and NNLO. In Sect. 4 we summarise our findings. In the Appendix we report the explicit expressions of the resummation coefficients required for the calculation and highlight ensuing dynamical features related to azimuthal correlations and asymmetries.

Transverse-momentum resummation and q T subtraction formalism for QQF production
We consider the associated production of a heavy-quark pair (QQ) and an arbitrary colourless system F in hadron collisions. The system F consists of one or more colourless particles, such as vector bosons, Higgs bosons, and so forth. We denote by M and q T the invariant mass and transverse momentum of the QQ F system, respectively.
At small values of q T (i.e., q T M) the perturbative QCD computation for this class of production processes is affected by large logarithmic terms of the type ln n (M/q T ). These terms can formally be resummed to all perturbative orders.
In the case of QQ production (i.e., no accompanying system F) the transverse-momentum resummation formalism was developed in Ref. [38] at arbitrary logarithmic accuracy, and by including the azimuthal-correlation contributions (see also Refs. [36,37] for the azimuthally averaged case up to next-to-leading logarithmic accuracy). The resummation formalism of Ref. [38] can be extended to QQ F associated production since the system F is formed by colourless particles. The extension is straightforward at the formal level and, at the practical level, it requires the explicit computation of process-dependent resummation factors at the necessary perturbative order. In the following we briefly summarize the key steps and the ingredients that are involved in this extension.
Transverse-momentum resummation is performed through Fourier transformation from impact parameter space, where the impact parameter vector b is the Fourier conjugated variable to the transverse-momentum vector q T . The general resummation formula for both QQ and QQ F production is given in Eq. (5) of Ref. [38]. The process-dependent contributions to this resummation formula are the LO partonic cross section dσ (0) cc (c = q,q, g) and the resummation factor H , while all the other contributions are process independent. These process-independent contributions and the corresponding resummation coefficients are the same terms that control the production of a colourless high-mass system [28] (see below). The factor H has an all-order processindependent structure (see Eqs. (10)- (16) and (26) in Ref. [38]) that is controlled by resummation coefficients that can be explicitly computed at the required perturbative order. These resummation coefficients are (i) the soft anomalous dimension matrix t , (ii) the radiative factor D and (iii) the subtraction operator I cc→ F .
(i) The resummation factor (see Eqs. (15)- (18) in Ref. [38]) depends on the soft anomalous dimension matrix t , whose perturbative expansion in the QCD coupling α S reads (ii) The term also depends on the radiative factor D(b, α S , { p i }), which embodies azimuthal correlations of soft origin and, therefore, it depends on the directionb of the impact parameter vector b. Its perturbative expansion reads with the constraint where · · · av. denotes the azimuthal average overb (i.e., the average over the azimuthal angle φ(b) of the transverse vector b). (iii) The subtraction operator I cc→ F has the following perturbative expansion, where μ R is the renormalisation scale of the QCD coupling α S (μ 2 R ). Here F generically denotes the observed final-state system (i.e., F = QQ for heavy-quark pair production, or F = QQ F for the associated production process) with total invariant mass M. The operator I cc→ F embodies IR-divergent contributions that are regularized by the customary procedure of analytic continuation in d = 4 − 2 space-time dimensions. This subtraction operator contributes to the resummation factor H (see Eqs. (12) and (13) in Ref. [38], and Eqs. (8) and (13) in the following) through the definition of the (IR-finite) hard-virtual amplitude M cc→ F (see Eq. (26) in Ref. [38] and Eq. (15) in the following) of the partonic production process cc → F.
The general transverse-momentum resummation formula for F = QQ, QQ F production involves a sole additional ingredient that is process dependent, namely the scattering amplitude M cc→ F of the partonic production process cc → F.
The resummation quantities t , D and I cc→ F have a 'minimal' process dependence, which has a soft origin: they depend on the momenta p i and colour charges T i of the colour-charged partons of the process cc → F (namely, the colliding partons c andc and the produced heavy quarks and antiquarks). Such dependence is simply denoted by the argu- (2) and (4). We also recall that t , D and I cc→ F are actually colour-space operators that act on the colour indices of the corresponding partons. The explicit expressions of the perturbative terms t , D (1) and I (1) cc→ F for QQ production were presented in Ref. [38], and we report the corresponding expressions for QQ F production in the Appendix of this paper.
According to the q T subtraction method [27], the formulation of transverse-momentum resummation for QQ F production allows us to write the (N)NLO partonic cross section NLO perturbative expansion (see, e.g., Ref. [44]) of the resummation formula of the logarithmically enhanced contributions to the corresponding q T distribution [36][37][38]. The explicit form of dσ QQ F, CT (N )N L O can be completely worked out up to NNLO accuracy. It depends on the resummation coefficients that control transverse-momentum resummation for the production of a colourless final-state system and, additionally, on the first two coefficients (1) t and (2) t of the soft anomalous dimension matrix in Eq. (1).
The explicit expression of the coefficient (1) t for QQ F production is given in the Appendix. The expression of (2) t can be determined (see the Appendix) by exploiting the relation [38] between t and the IR singularities of the virtual scattering amplitude M cc→QQ F [45][46][47][48][49].
The IR-finite function H QQ F in Eq. (5) corresponds to the coefficient of the δ (2) (q T ) contribution in the expansion of the resummation formula. It reads [38] where the perturbative functions C 1 and C 2 are process independent and describe the emission of collinear radiation off the incoming partons. In Eq. (6) 6) and in the following (see Eqs. (9) and (12)) the explicit structure of the function H QQ F is presented by setting μ R = μ F = M. The exact dependence on μ R and μ F can straightforwardly be recovered by using renormalisation group invariance and evolution of the parton densities.
In the quark annihilation channel (c = q,q) the functions C 1 and C 2 do not depend on b, and the symbolic factor [(H D)C 1 C 2 ] cc;a 1 a 2 takes the form cc→QQ F | 2 in the denominator is the LO contribution to the squared amplitude |M cc→QQ F | 2 for the process cc → QQ F (note that the power p depends on the process, see Eq. (16)). The IR-finite hard-virtual amplitude M cc→QQ F in Eq. (8) is defined in terms of the allorder renormalised virtual amplitude M cc→QQ F through an appropriate subtraction of IR singularities (see Eq. (15)). By using Eq. (3) the contribution of the azimuthal factor D to Eq. (6) becomes trivial, and we obtain In the gluon fusion channel (c = g) the collinear functions C 1 and C 2 can be decomposed as [50] C μν where the tensors d μν and D μν , which multiply the helicityconserving and helicity-flip components C ga and G ga , read Therefore, the function H where The functions C ca (z; α S ) (c = q,q, g) and G ga (z; α S ) have perturbative expansions The helicity-conserving coefficients C (n) ca (z) are known up to n = 2 [51][52][53][54], and they are the same that contribute to Higgs boson [27] and vector-boson [55] production. Recently, their computation has been extended to the third order (n = 3) [31][32][33]. The helicity-flip coefficients G (n) ga (z) are known up to n = 2 [56,57].
The hard-virtual amplitude M cc→QQ F in Eq. (9) and (13) is expressed in terms of the all-order renormalised virtual amplitude M cc→QQ F as where I cc→QQ F is the subtraction operator whose perturbative expansion is given in Eq. (4). The general expression of the first order coefficient I (1) cc→QQ F in Eq. (4) is known (see Appendix), while the result for the second-order coefficient is available only in the case of heavy-quark production [40,58,59].
The quantity M cc→QQ F ({ p i }) on the right-hand side of Eq. (15) is the renormalised on-shell scattering amplitude and has the perturbativeexpansion The perturbative expansion of M cc→QQ F is completely analogous to that in Eq. (16) Having discussed the perturbative ingredients entering the function H QQ F cc;a 1 a 2 , we can now examine its NLO and NNLO expansions. By inspecting Eqs. (9) and (12) we see that the NLO truncation of H QQ F cc;a 1 a 2 receives contributions only from the tree-level and one-loop hard-virtual amplitudes M (0) cc→QQ F and M (1) cc→QQ F , and from the first-order helicity-conserving coefficients C (1) ca (z). Indeed, at this perturbative order the azimuthally dependent terms in D (1) (2) cc→QQ F , which in turn requires the knowledge of the renormalised two-loop amplitude, and of the subtraction operator I (2) cc→QQ F . The computation of the two-loop amplitude for QQ F production is at the frontier of current techniques, and, moreover, I (2) cc→QQ F is also not known yet. However, the NNLO contribution to H QQ F cc;a 1 a 2 can be completely determined for all the flavour off-diagonal partonic channels (a 1 , a 2 ) = (c,c). The perturbative ingredients entering the calculation in the quark-antiquark annihilation channel (see Eq. (9)) are the corresponding tree-level (M (0) cc→QQ F ) and one-loop ( M (1) cc→QQ F ) hard-virtual amplitudes and the first-and second-order helicity-conserving coefficients C (1) ab (z) and C (2) ab (z). In the gluon fusion channel, the azimuthally dependent first-order coefficients of soft (D (1) ) and collinear (G (1) ga ) origin are also required. Indeed at NNLO such coefficients produce [38] non-vanishing mixed collinear-collinear and soft-collinear contributions in the expansion of Eq. (12). The general expression of the first-order coefficient D (1) . The evaluation of the ensuing NNLO contributions can be performed by computing the corresponding spin and colour-correlated squared tree-level amplitudes for the process gg → QQ F.
In summary, the current knowledge of transversemomentum resummation for high-mass systems containing a heavy-quark pair and a colour singlet system allows us to use Eq. (5) to obtain the complete NLO corrections for this class of processes plus the NNLO corrections in the flavour off-diagonal partonic channels.

Results for tt H production
Having discussed the content of Eq. (5), we are in a position to apply it to tt H production and to obtain the complete NLO results plus the NNLO corrections in all the flavour offdiagonal partonic channels. Our NLO implementation of the calculation has the main purpose of illustrating the applicability of the q T subtraction method to tt H production and, in particular, of cross-checking the q T subtraction methodology by numerical comparisons with NLO calculations performed by using more established NLO methods. Our NNLO results on tt H production represent a first step (due to the missing flavour diagonal partonic channels) towards the complete NNLO calculation for this production process.
Our results are obtained with two independent computations, which show complete agreement. In the first computation, up to NLO, we use the phase space generation routines from the MCFM program [66], suitably modified for q T subtraction along the lines of the corresponding numerical programs for Higgs boson [27] and vector-boson [55] production. The numerical integration is carried out using the Cuba library [67]. At NNLO accuracy the tt H +jet cross section is evaluated by using the Munich code, 2 which provides a fully automated implementation of the NLO dipole subtraction formalism [68][69][70] and an efficient phase space integration. The remaining flavour off-diagonal contributions at NNLO are evaluated with a dedicated fortran implementation. The second computation is directly implemented within the Matrix framework [71], suitably extended to tt H production. In both implementations all the required treelevel and one-loop amplitudes are obtained with OpenLoops [60-62], including the tree-level spin-and colour-correlated amplitudes required to evaluate the contributions in Eq. (12).
In order to numerically evaluate the contribution in the square bracket of Eq. (5), a technical cut-off r cut is introduced on the dimensionless variable q T /M, where M is the invariant mass of the tt H system. The final result, which corresponds to the limit r cut → 0, is extracted by computing the cross section at fixed values of r cut in the range [0.01%, r max ]. Table 1 The tt H total cross section at LO and NLO, and its NNLO corrections in the flavour off-diagonal partonic channels. The numerical uncertainties at LO and NLO (Madgraph5_aMC@NLO, Matrix) are due to numerical integration, while at NLO (q T subtraction) and NNLO they also include the systematics uncertainty from the r cut → 0 extrapolation Quadratic least χ 2 fits are performed for different values of r max ∈ [0.5%, 1%]. The extrapolated value is then extracted from the fit with lowest χ 2 /degrees-of-freedom, and the uncertainty is estimated by comparing the results obtained by the different fits. This procedure is the same as implemented in matrix [71] and it has been shown to provide a conservative estimate of the systematic uncertainty in the q T subtraction procedure for various processes (see Sec. 7 in Ref. [71]). We consider pp collisions at the centre-of-mass energies √ s = 13 TeV and √ s = 100 TeV. We use the NNPDF31 [72] parton distribution functions (PDFs) with the QCD running coupling α S evaluated at each corresponding order (i.e., we use (n + 1)-loop α S at N n LO, with n = 1, 2). The pole mass of the top quark is m t = 173.3 GeV, the Higgs boson mass m H = 125 GeV, and the Fermi constant G F = 1.16639 × 10 −5 GeV −2 . The renormalisation and factorization scales, μ R and μ F , are fixed at μ R = μ F = (2m t + m H )/2. Our predictions for the LO and NLO cross sections and for the NNLO corrections in the flavour offdiagonal channels are presented in Table 1 together with their uncertainties due to the numerical integration and the extrapolation to r cut → 0, computed as explained above. The NLO cross section computed with q T subtraction is compared with the result obtained with Madgraph5_aMC@NLO [73], which uses FKS subtraction [74,75] and with the corresponding result obtained with Matrix, which implements dipole subtraction [68][69][70].
We start our discussion from the NLO results. The NLO corrections increase the LO result by 27% (31%) at √ s = 13 TeV ( √ s = 100 TeV). The flavour offdiagonal qg +qg channel contributes about 15% (23%) of the total NLO correction. As expected, from Table 1 we observe excellent agreement between the NLO cross section obtained with Madgraph5_aMC@NLO and Matrix. The result obtained with q T subtraction also agrees with Mad-graph5_aMC@NLO and Matrix results. The quality of the r cut → 0 extrapolation can be assessed by studying the behavior of the cross section at fixed values of r cut . In Fig. 1 we investigate this behavior and show also the (r cut independent) NLO result obtained with Matrix, by using dipole subtraction, and Madgraph5_aMC@NLO, by using FKS subtraction. As expected, the r cut dependence is linear [76,77], contrary to what happens in the case of the production of a colourless final-state system (see Sec. 7 of Ref. [71]), where the power-like dependence of the total cross section on r cut is known [78][79][80] to be quadratic (modulo logarithmic enhancements).
In Fig. 2 we present the NLO results for several differential distributions at √ s = 13 TeV and compare them with those obtained by using Madgraph5_aMC@NLO. In particular we consider the transverse-momentum (top left) and rapidity (top right) distributions of the Higgs boson, the transversemomentum (center left) and rapidity (center right) distributions of the top quark, and the invariant-mass distributions of the top-quark pair (bottom left) and of the tt H system (bottom right). We find excellent agreement between the two Fig. 1 The r cut dependence (data points) at √ s = 13 TeV (left) and 100 TeV (right) of the NLO total cross section computed by using q T subtraction. The bands show the extrapolated value at r cut → 0 and the NLO results from Madgraph5_aMC@NLO (using FKS subtraction) and Matrix (using dipole subtraction) Fig. 2 The NLO results of Madgraph5_aMC@NLO for the cross section dependence on several kinematic variables at √ s = 13 TeV. The lower panels show the relative comparison with the corresponding results obtained by using q T subtraction calculations, with bin-wise uncertainties at the percent-level or below. Our results are obtained by using a fixed value of r cut , r cut = 0.1%, which, as suggested by Fig. 1, already provides a good estimate of the NLO cross section. A similar level of agreement is found with the results obtained with dipole subtraction. We checked that the agreement also holds for different values of μ R and μ F .
We now move to considering the NNLO contributions to the total cross section. In Table 1 we report our results for the O(α 4 S ) contributions to the NNLO cross section from the flavour off-diagonal partonic channels a 1 a 2 → tt H + X . The contribution from all the channels with a 1 a 2 = qg,qg is labelled by the subscript qg, and the contribution from all the channels with a 1 a 2 = qq,qq, qq ,qq , qq ,qq (q = q ) is labelled by the subscript q(q)q . We see that the NNLO corrections from both contributions are very small, at the few per mille level of the NLO cross section. At √ s = 13 TeV they contribute with similar size and opposite sign, and, therefore, their overall quantitative effect in this setup is completely negligible. We also see that the numerical uncertainty of the NNLO correction in the qg channel is rather large. This is due to a cancellation between the two terms in Eq. (5): the term H ⊗ dσ , which is r cut independent, and the term in the square bracket, which depends on r cut . This cancellation is observed at both √ s = 13 TeV and √ s = 100 TeV, but it is particularly severe at √ s = 13 TeV, downgrading the numerical precision that can be obtained for the relative correction. Similar effects were observed for tt production in Refs. [39,40].
In Fig. 3 we show the r cut -dependence of the NNLO corrections σ of the flavour off-diagonal channels to the total cross section for tt H production. The result is normalised to our extrapolation σ NNLO at r cut → 0. In the qg channel at √ s = 13 TeV the extrapolation is particularly delicate, due to the cancellation discussed above. For some of the channels, the first few points at low r cut values show relatively large instabilities, in particular again for the qg channel at √ s = 13 in Fig. 3  (top-left). However, these points are not dropped in the fit, which is dominated by the behaviour at r cut > 0.1%. The r cut dependence confirms that our calculation can control the NNLO contributions in the off-diagonal partonic channels at the few percent level. Comparing with the r cut behaviour for tt production (see Fig. 1 of Ref. [40]), the behaviour in Fig. 3 for tt H production is qualitatively and quantitatively similar, and we do expect the extrapolation at r cut → 0 in the flavour diagonal channels to work in a similar way for both production processes. Based on the experience with the NNLO calculations for heavy-quark production [40][41][42], this should be fully sufficient to obtain precise NNLO results by using the q T subtraction method once the presently unknown Fig. 3 The NNLO contribution σ of the qg (top) and q(q)q (bottom) partonic channels to the total cross section at √ s = 13 TeV (left) and 100 TeV (right). The r cut dependence of σ is normalized to its extrapolation σ NNLO at r cut → 0 soft contributions and the two-loop amplitudes become available.

Summary
In this paper we have considered the associated production of the SM Higgs boson with a top-quark pair, and, more generally, processes in which heavy-quark pairs are produced in association with a colourless final-state system F. We have pointed out that the transverse-momentum resummation formalism developed for QQ production in Ref. [38] can be extended to associated QQ F production. This extension, which requires the evaluation of the appropriate resummation coefficients at the necessary perturbative accuracy, is also sufficient to apply the q T subtraction method to this class of processes.
Using the resummation coefficients presented in this paper and the current knowledge of scattering amplitudes, it is possible to apply the q T subtraction formalism to QQ F production up to NLO and to obtain the NNLO corrections in all the flavour off-diagonal partonic channels.
We have implemented for the first time the q T subtraction formalism for tt H production, and we have presented first quantitative results at NLO and NNLO. The calculation is accurate at NLO in QCD, and the NNLO corrections have been computed for the flavour off-diagonal partonic channels. At NLO we have checked the correctness of our implementation by comparing with the results obtained by using tools that are based on established subtraction methods. We found complete agreement for the total cross section and for single-differential distributions. Within the setup that we have considered, we have found that the NNLO contribution of the off-diagonal partonic channels to the total cross section has a very small quantitative effect. The extension of this calculation to the diagonal channels requires further theoretical work to compute the two-loop virtual amplitudes, and the NNLO soft contributions to the resummation coefficients.
Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors' comment: There are no external data associated with the manuscript.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 . where the coefficients γ i (i = q,q, g) originate from collinear radiation and read γ q = γq = 3C F /2 and γ g = (11C A − 2N f )/6, N f being the number of flavours of massless quarks. The function F (1) where with The function Li 2 (z) in Eq. (23) is the customary dilogarithm function, Li 2 (z) = − 1 0 dt t ln(1 − zt). The coefficient D (1) with the function D jk that is given in terms of the following one-fold integral representation: where the four-vector w We note that the expressions of the resummation coefficients in Eqs. (19), (21), (22) and (25) for the process in Eq. (18) have no explicit dependence on the flavour of the (massless and heavy) quarks and on the colourless system F. In particular, there is no dependence on the quantum numbers of F and on its momentum p F (though the dependence on p F enters implicitly through momentum conservation, Using Eqs. (19), (21), (22) and (25), we also note that we can recover the resummation coefficients of Ref. [38] for the simpler process of heavy-quark pair production 3 , c( p 1 )c( p 2 ) → Q( p 3 )Q( p 4 ). To this purpose we can simply use the corresponding constraints m 3 = m 4 , T 1 + T 2 = −(T 3 +T 4 ) and, importantly, the constraint p 3T = −p 4T that follows from momentum conservation (i.e., p F = 0). Indeed, to a large extent, the main difference between heavy-quark pair production and the process in Eq. (18) is due to the fact that the N transverse momenta p j T (3 ≤ j ≤ N ) in Eq. (18) are independent kinematical variables (since p FT = 0). This main kinematical difference leads to technical computational complications (e.g., in the analytic evaluation of the resummation coefficients) and to new dynamical features.
One of these dynamical features regards azimuthal correlations. The b-space resummation coefficient D (1) (which is due [38] to soft-parton radiation from the final state and from initial-state/final-state interferences) depends on the relative azimuthal angles φ jb = φ(p j T ) − φ(b). This dependence produces ensuing azimuthal correlations [38] with respect to the observable angles φ jq = φ(p j T ) − φ(q T ) of the q Tdifferential cross section at q T = 0. In the limit q T → 0 the azimuthal-correlation cross section is divergent orderby-order in QCD perturbation theory (transverse-momentum resummation eventually makes the cross section finite [76]).
In the case of heavy-quark pair production, D (1) depends on cos 2 φ 3b = cos 2 φ 4b [38], and this implies that it produces azimuthal-correlation divergences only in the case of even harmonics (i.e., harmonics with respect to cos k φ 3q with k = 2, 4, 6, . . .) [76]. This dependence of D (1) is due to the kinematical relation p 3T −p 4T at q T → 0. In the case of the process in Eq. (18), the transverse momenta p j T are kinematically independent even if q T → 0, and it turns out that D (1) depends on both cos 2 φ jb and cos φ jb (the dependence on cos φ jb is due to the contributions that are proportional to 'iπ ' in the right-hand side of Eqs. (25) and (26)). This dependence implies that the process in Eq. (18) is affected by lowest-order (and higherorder) divergent azimuthal correlations for both even and odd harmonics (i.e., harmonics with respect to cos k φ 3q with k = 1, 3, 5, . . . ). A quite general discussion of azimuthal correlations at small q T is presented in Ref. [76], where is also pointed out that lowest-order divergent odd harmonics can occur in other processes, such as V + jet and dijet production.
We conclude this Appendix by recalling [38] that the second-order term (2) t of the soft anomalous dimension t (see Eq. (1)) for the process in Eq. (18) is related to the corresponding term of the anomalous dimension matrix (μ) that controls the QCD IR divergences of scattering amplitudes with massive external particles [45][46][47][48][49]. We have t (α S ; {p i }) = (1) t ({ p i }) and F (1) t ({ p i } are given in Eqs. (19) and (22), and sub. is given below. The perturbative expansion of the right-hand side of Eq. (28) includes both the first-order and second-order terms (1) t and (2) t (obviously, sub. = 2(α S /π ) (1) t + O(α 2 S )), while terms at O(α 3 S ) and beyond are neglected. The 'subtracted' anoma- where the terms on the right-hand side are written by exactly using the notation of Eq. (5) of Ref. [48]. The term (μ) is the anomalous-dimension matrix that controls the IR divergences of the scattering amplitude M cc→Q 3Q4 ...F for the process in Eq. (18), while the square-bracket term on the right-hand side of Eq. (29) is the corresponding expression of (μ) for a generic process cc → F, where the system F is colourless. The expressions of (μ), γ cusp (α S ) and γ c (α S ) up to O(α 2 S ) are explicitly given in Ref. [48].