${t {{\bar{t}}}H}$ production at NNLO: the flavour off-diagonal channels

Abstract

We consider QCD radiative corrections to the associated production of a heavy-quark pair ($Q{{\bar{Q}}}$) with a generic colourless system F at hadron colliders. We discuss the resummation formalism for the production of the $Q{{\bar{Q}}}F$ system at small values of its total transverse momentum $q_T$. 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 $q_T$ subtraction formalism to fixed-order calculations for this class of processes. We use the $q_T$ 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.

1 Introduction

The observation of Higgs boson production in association with a top quark–antiquark ($t {{\bar{t}}}$) 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 top–antitop quark pair and a Higgs boson ($t {{\bar{t}}}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 resummation of soft-gluon contributions close to the partonic kinematical threshold was considered in Refs. [14,15,16,17,18,19]. Full off-shell calculations with decaying top quarks were presented at NLO QCD [20] and NLO QCD+EW [21]. The uncertainty of current theoretical predictions for $t {{\bar{t}}}H$ cross sections is at the ${{\mathcal {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 $t {{\bar{t}}}H$ production. At the partonic level, the NNLO calculation of $t {{\bar{t}}}H$ production requires the evaluation of tree-level 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 tree-level and one-loop scattering amplitudes can nowadays be evaluated with automated tools. The two-loop amplitude for $t {{\bar{t}}}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 considering and explicitly computing the $q_T$ distribution of the produced final-state system in the limit $q_T \rightarrow 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 \rightarrow 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.¹ 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 heavy-quark 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-to-back 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 $t {{\bar{t}}}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 $t {{\bar{t}}}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 $q{{\bar{q}}} \rightarrow t {{\bar{t}}}H+X$ and $gg \rightarrow t {\bar{t}}H+X$ (X denotes the unobserved inclusive final state) are not yet known. Their evaluation requires the two-loop amplitudes for the partonic processes $q{{\bar{q}}} \rightarrow t {{\bar{t}}}H$ and $gg \rightarrow t {\bar{t}}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 \rightarrow t {{\bar{t}}}H+X$, with $ab= qg ({{\bar{q}}}g), qq ({{\bar{q}}}{{\bar{q}}}), qq' ({{\bar{q}}}{{\bar{q}}}'), q{{\bar{q}}}' ({{\bar{q}}} q')$ (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 ttH 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.

2 Transverse-momentum resummation and ${\mathbf{q}_{\mathbf{T}}}$ subtraction formalism for ${\mathbf{Q}{{\bar{\mathbf{Q}}}}{} \mathbf{F}}$ production

We consider the associated production of a heavy-quark pair ($Q{\bar{Q}}$) 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 ${\mathbf{q}_{\mathbf{T}}}$ the invariant mass and transverse momentum of the $Q{{\bar{Q}}}F$ system, respectively.

At small values of $q_T$ (i.e., $q_T \ll 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 $Q{{\bar{Q}}}$ 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 $Q{\bar{Q}}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 $\mathbf{b}$ is the Fourier conjugated variable to the transverse-momentum vector ${\mathbf{q}_{\mathbf{T}}}$. The general resummation formula for both $Q{{\bar{Q}}}$ and $Q{{\bar{Q}}}F$ production is given in Eq. (5) of Ref. [38]. The process-dependent contributions to this resummation formula are the LO partonic cross section $\left[ d\sigma _{c{{\bar{c}}}}^{(0)} \right]$ ($c=q,{\bar{q}},g$) and the resummation factor $\bigl ( \mathbf{H} {\varvec{\Delta }} \bigr )$, 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 $\bigl ( {\mathbf{H} {\varvec{\Delta }} } \bigr )$ has an all-order process-independent 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 ${\varvec{\Gamma }}_t$, (ii) the radiative factor $\mathbf{D}$ and (iii) the subtraction operator $\widetilde{\mathbf{I}}_{c{{\bar{c}}}\rightarrow {{\widetilde{F}}}}$.

• (i) The resummation factor $\varvec{\Delta }$ (see Eqs. (15)–(18) in Ref. [38]) depends on the soft anomalous dimension matrix ${\varvec{\Gamma }}_t$, whose perturbative expansion in the QCD coupling $\alpha _{\mathrm {S}}$ reads

  $$\begin{aligned} {\varvec{\Gamma }}_t (\alpha _{\mathrm {S}},\{p_i\})= & {} \frac{\alpha _{\mathrm {S}}}{\pi }\, {\varvec{\Gamma }}_t^{(1)}(\{p_i\})+\left( \frac{\alpha _{\mathrm {S}}}{\pi }\right) ^2 {\varvec{\Gamma }}_t^{(2)}(\{p_i\})\nonumber \\&+{{\mathcal {O}}}(\alpha _{\mathrm {S}}^3). \end{aligned}$$

  (1)

• (ii) The term $\varvec{\Delta }$ also depends on the radiative factor $\mathbf{D}({{\hat{\mathbf{b}}}},\alpha _{\mathrm {S}},\{p_i\})$, which embodies azimuthal correlations of soft origin and, therefore, it depends on the direction ${\hat{\mathbf{b}}}$ of the impact parameter vector $\mathbf{b}$. Its perturbative expansion reads

  $$\begin{aligned} \mathbf{D}({\hat{\mathbf{b}}},\alpha _{\mathrm {S}},\{p_i\})) = 1+\frac{\alpha _{\mathrm {S}}}{\pi }\, \mathbf{D}^{(1)}({\hat{\mathbf{b}}},\{ p_i\})+\mathcal{O}(\alpha _{\mathrm {S}}^2), \end{aligned}$$

  (2)

  with the constraint

  $$\begin{aligned} \langle \mathbf{D}({\hat{\mathbf{b}}},\alpha _{\mathrm {S}},\{p_i\}))\rangle _{\mathrm{av}.} = 1, \end{aligned}$$

  (3)

  where $\langle \cdots \rangle _{\mathrm{av.}}$ denotes the azimuthal average over ${\hat{\mathbf{b}}}$ (i.e., the average over the azimuthal angle ${\phi (\mathbf{b})}$ of the transverse vector $\mathbf{b})$.

• (iii) The subtraction operator $\widetilde{\mathbf{I}}_{c{{\bar{c}}}\rightarrow {{\widetilde{F}}}}$ has the following perturbative expansion,

  $$\begin{aligned}&\widetilde{\mathbf{I}}_{c{{\bar{c}}}\rightarrow {{\widetilde{F}}}}(\alpha _{\mathrm {S}}(M^2), \epsilon ;\{p_i\}) \nonumber \\&\quad = \sum _{n=1}^{\infty } \left( \frac{\alpha _{\mathrm {S}}(\mu _R^2)}{2\pi } \right) ^{\!\!n} \;\widetilde{\mathbf{I}}^{(n)}_{c{{\bar{c}}}\rightarrow {\widetilde{F}}} (\epsilon ,M^2/\mu _R^2;\{p_i\}), \end{aligned}$$

  (4)

  where $\mu _R$ is the renormalisation scale of the QCD coupling $\alpha _{\mathrm {S}}(\mu _R^2)$. Here ${{\widetilde{F}}}$ generically denotes the observed final-state system (i.e., ${{\widetilde{F}}}= Q{{\bar{Q}}}$ for heavy-quark pair production, or ${{\widetilde{F}}}= Q{{\bar{Q}}}F$ for the associated production process) with total invariant mass M. The operator $\widetilde{\mathbf{I}}_{c{{\bar{c}}}\rightarrow {{\widetilde{F}}}}$ embodies IR-divergent contributions that are regularized by the customary procedure of analytic continuation in $d=4-2\epsilon$ space-time dimensions. This subtraction operator contributes to the resummation factor $\mathbf{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 ${\widetilde{{\mathcal {M}}}}_{c{{\bar{c}}} \rightarrow {{\widetilde{F}}}}$ (see Eq. (26) in Ref. [38] and Eq. (15) in the following) of the partonic production process $c{{\bar{c}}} \rightarrow {{\widetilde{F}}}$.

The general transverse-momentum resummation formula for ${\widetilde{F}}= Q{{\bar{Q}}}, Q{{\bar{Q}}}F$ production involves a sole additional ingredient that is process dependent, namely the scattering amplitude ${{\mathcal {M}}}_{c{{\bar{c}}} \rightarrow {{\widetilde{F}}}}$ of the partonic production process $c{{\bar{c}}} \rightarrow {{\widetilde{F}}}$.

The resummation quantities ${\varvec{\Gamma }}_t$, $\mathbf{D}$ and $\widetilde{\mathbf{I}}_{c{{\bar{c}}}\rightarrow {{\widetilde{F}}}}$ have a ‘minimal’ process dependence, which has a soft origin: they depend on the momenta $p_i$ and colour charges $\mathbf{T}_i$ of the colour-charged partons of the process $c{{\bar{c}}} \rightarrow {{\widetilde{F}}}$ (namely, the colliding partons c and ${\bar{c}}$ and the produced heavy quarks and antiquarks). Such dependence is simply denoted by the argument $\{ p_i \}$ in Eqs. (1), (2) and (4). We also recall that ${\varvec{\Gamma }}_t$, $\mathbf{D}$ and $\widetilde{\mathbf{I}}_{c{{\bar{c}}}\rightarrow {{\widetilde{F}}}}$ are actually colour-space operators that act on the colour indices of the corresponding partons. The explicit expressions of the perturbative terms ${\varvec{\Gamma }}_t^{(1)}$, ${\varvec{\Gamma }}_t^{(2)}$, $\mathbf{D}^{(1)}$ and $\widetilde{\mathbf{I}}_{c{{\bar{c}}}\rightarrow {\widetilde{F}}}^{(1)}$ for $Q{{\bar{Q}}}$ production were presented in Ref. [38], and we report the corresponding expressions for $Q{{\bar{Q}}}F$ production in the Appendix of this paper.

According to the $q_T$ subtraction method [27], the formulation of transverse-momentum resummation for $Q{{\bar{Q}}}F$ production allows us to write the (N)NLO partonic cross section $d{\sigma }^{Q{{\bar{Q}}}F}_{(N)NLO}$ as

$$\begin{aligned} d{{{\hat{\sigma }}}}^{Q{{\bar{Q}}}F}_{(N)NLO}= & {} \mathcal{H}^{Q{{\bar{Q}}}F}_{(N)NLO}\otimes d{{{\hat{\sigma }}}}^{Q{{\bar{Q}}}F}_{LO}\nonumber \\&+\left[ d{{{\hat{\sigma }}}}^{Q{{\bar{Q}}}F+\mathrm{jet}}_{(N)LO}- d{\hat{\sigma }}^{Q{{\bar{Q}}}F, \, CT}_{(N)NLO}\right] , \end{aligned}$$

(5)

where $d{\sigma }^{Q{{\bar{Q}}}F+\mathrm{jet}}_{(N)LO}$ is the $Q{\bar{Q}}F$+jet cross section at (N)LO accuracy. To apply Eq. (5) at NLO, the LO cross section $d{\sigma }^{Q{\bar{Q}}F+\mathrm{jet}}_{LO}$ can be directly obtained by integrating the corresponding tree-level scattering amplitudes. To apply Eq. (5) at NNLO, $d{\sigma }^{Q{{\bar{Q}}}F+\mathrm{jet}}_{NLO}$ can be evaluated by using any available NLO method to handle and cancel the corresponding IR divergences, if the relevant tree-level and one-loop QCD amplitudes are available. Therefore, $d{\sigma }^{Q{{\bar{Q}}}F+\mathrm{jet}}_{(N)LO}$ is IR finite, provided $q_T \ne 0$. The square bracket term of Eq. (5) is IR finite in the limit $q_T \rightarrow 0$, but its individual contributions, $d{\sigma }^{Q{{\bar{Q}}}F+\mathrm{jet}}_{(N)LO}$ and $d{\sigma }^{Q{{\bar{Q}}}F, \, CT}_{(N)NLO}$, are separately divergent. The IR-subtraction counterterm $d{\sigma }^{Q{{\bar{Q}} F}, \,CT}_{(N)NLO}$ is obtained from the (N)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{\sigma }^{Q{{\bar{Q}}}F, \,CT}_{(N)NLO}$ 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 ${\varvec{\Gamma }}_t^{(1)}$ and ${\varvec{\Gamma }}_t^{(2)}$ of the soft anomalous dimension matrix in Eq. (1).

The explicit expression of the coefficient ${\varvec{\Gamma }}_t^{(1)}$ for $Q{{\bar{Q}}} F$ production is given in the Appendix. The expression of ${\varvec{\Gamma }}_t^{(2)}$ can be determined (see the Appendix) by exploiting the relation [38] between ${\varvec{\Gamma }}_t$ and the IR singularities of the virtual scattering amplitude ${{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q{{\bar{Q}}}F}$ [45,46,47,48,49].

The IR-finite function ${{\mathcal {H}}}^{Q{{\bar{Q}}}F}$ in Eq. (5) corresponds to the coefficient of the $\delta ^{(2)}({\mathbf{q}_{\mathbf{T}}})$ contribution in the expansion of the resummation formula. It reads [38]

$$\begin{aligned} {{\mathcal {H}}}^{Q{{\bar{Q}}}F}_{c{{\bar{c}}};a_1a_2}=\langle \left[ (\mathbf{H}\, \mathbf{D})C_1C_2\right] _{c{{\bar{c}}};a_1a_2}\rangle _{\mathrm{av.}}, \end{aligned}$$

(6)

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) we have explicitly denoted the parton indices $\{c{{\bar{c}}},a_1,a_2\}$ that are implicit in Eq. (5). The indices c and ${{\bar{c}}}$ correspond to the incoming partons of the LO partonic cross section $d{\hat{\sigma }}^{Q{{\bar{Q}}}F}_{LO}$. The indices $a_1$ and $a_2$ are those of the parton densities $f_{a_1}$ and $f_{a_2}$ of the colliding hadrons. The partonic cross section in Eq. (5) depends on the renormalisation scale $\mu _R$ of $\alpha _{\mathrm {S}}$ and on the factorisation scale $\mu _F$ of the parton densities. In Eq. (6) and in the following (see Eqs. (9) and (12)) the explicit structure of the function ${{\mathcal {H}}}^{Q{{\bar{Q}}}F}$ is presented by setting $\mu _R=\mu _F=M$. The exact dependence on $\mu _R$ and $\mu _F$ can straightforwardly be recovered by using renormalisation group invariance and evolution of the parton densities.

In the quark annihilation channel ($c=q,{{\bar{q}}}$) the functions $C_1$ and $C_2$ do not depend on $\mathbf{b}$, and the symbolic factor $\left[ (\mathbf{H}\, \mathbf{D})C_1C_2\right] _{c{{\bar{c}}};a_1a_2}$ takes the form

$$\begin{aligned} \left[ (\mathbf{H}\, \mathbf{D})C_1C_2\right] _{c{{\bar{c}}};a_1a_2}=\left( \mathbf{H}\, \mathbf{D}\right) _{c{{\bar{c}}}} C_{ca_1}C_{{\bar{c}}a_2}\quad (c=q,{{\bar{q}}}) \end{aligned}$$

(7)

with

$$\begin{aligned} \left( \mathbf{H}{} \mathbf{D}\right) _{c{{\bar{c}}}}=\frac{\langle \widetilde{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q{{\bar{Q}}}F}|\mathbf{D}| \widetilde{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q{{\bar{Q}}}F}\rangle }{\alpha _{\mathrm {S}}^p(M^2)\,|{{\mathcal {M}}}_{c{\bar{c}}\rightarrow Q{{\bar{Q}}}F}^{(0)}(\{p_i\})|^2}\quad (c=q,{{\bar{q}}}). \end{aligned}$$

(8)

The factor $\alpha _{\mathrm {S}}^p|{{\mathcal {M}}}^{(0)}_{c{{\bar{c}}}\rightarrow Q{{\bar{Q}}}F}|^2$ in the denominator is the LO contribution to the squared amplitude $|{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q{{\bar{Q}}}F}|^2$ for the process $c{\bar{c}}\rightarrow Q{{\bar{Q}}}F$ (note that the power p depends on the process, see Eq. (16)). The IR-finite hard-virtual amplitude $\widetilde{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q{{\bar{Q}}}F}$ in Eq. (8) is defined in terms of the all-order renormalised virtual amplitude ${{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}$ through an appropriate subtraction of IR singularities (see Eq. (15)). By using Eq. (3) the contribution of the azimuthal factor $\mathbf{D}$ to Eq. (6) becomes trivial, and we obtain

$$\begin{aligned}&{{\mathcal {H}}}^{Q{{\bar{Q}}}F}_{c{\bar{c}};a_1a_2}=\frac{\langle \widetilde{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q{\bar{Q}}F}|\widetilde{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q{\bar{Q}}F}\rangle }{\alpha _{\mathrm {S}}^p(M^2)\,|{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q{\bar{Q}}F}^{(0)}(\{p_i\})|^2}\, C_{ca_1}\, C_{{{\bar{c}}}a_2}\nonumber \\&(c=q,{\bar{q}}). \end{aligned}$$

(9)

In the gluon fusion channel ($c=g$) the collinear functions $C_1$ and $C_2$ can be decomposed as [50]

$$\begin{aligned} C_{ga}^{\mu \nu }(z,p_1,p_2,{\hat{\mathbf{b}}};\alpha _{\mathrm {S}})= & {} d^{\mu \nu }(p_1,p_2)\,C_{ga}(z;\alpha _{\mathrm {S}})\nonumber \\&+D^{\mu \nu }(p_1,p_2, {\hat{\mathbf{b}}})\,G_{ga}(z;\alpha _{\mathrm {S}}), \end{aligned}$$

(10)

where the tensors $d^{\mu \nu }$ and $D^{\mu \nu }$, which multiply the helicity-conserving and helicity-flip components $C_{ga}$ and $G_{ga}$, read ($b^\mu =(0,\mathbf{b},0)$ with $b_\mu b^\mu =-\mathbf{b}^2$)

$$\begin{aligned} d^{\mu \nu }(p_1,p_2)= & {} -g^{\mu \nu }+\frac{p_1^\mu p_2^{\nu }+p_1^\nu p_2^{\mu }}{p_1\cdot p_2},\nonumber \\ D^{\mu \nu }(p_1,p_2,{\hat{\mathbf{b}}})= & {} d^{\mu \nu }(p_1,p_2)-2\frac{b^\mu b^\nu }{\mathbf{b}^2}. \end{aligned}$$

(11)

Therefore, the function ${{\mathcal {H}}}^{Q{\bar{Q}}F}_{gg;a_1a_2}$ reads

$$\begin{aligned} {{\mathcal {H}}}^{Q{{\bar{Q}}}F}_{gg;a_1a_2}=\langle (\mathbf{H}\, \mathbf{D})_{gg;\mu _1\nu _1,\mu _2\nu _2}\, C_{ga_1}^{\mu _1\nu _1}({\hat{\mathbf{b}}}....)\, C_{ga_2}^{\mu _2\nu _2}({\hat{\mathbf{b}}}....)\rangle _{\mathrm{av.}},\nonumber \\ \end{aligned}$$

(12)

where

$$\begin{aligned}&\!\!\!\left( \mathbf{H} \,\mathbf{D} \right) _{gg;\mu _1 \,\nu _1, \mu _2 \,\nu _2 } \nonumber \\&\!\!\!\quad =\frac{\langle \,\widetilde{{\mathcal {M}}}_{gg \rightarrow Q {\bar{Q}}F}^{\nu _1^\prime \nu _2^\prime } \,| \, \mathbf{D} \, | \,\widetilde{{\mathcal {M}}}_{gg \rightarrow Q {{\bar{Q}}}F}^{\mu _1^\prime \mu _2^\prime } \, \rangle \;d_{\mu _1^\prime \mu _1} \;d_{\nu _1^\prime \nu _1} \;d_{\mu _2^\prime \mu _2} \; d_{\nu _2^\prime \nu _2} }{\alpha _{\mathrm {S}}^p(M^2)\,|{{\mathcal {M}}}_{gg\rightarrow Q {{\bar{Q}}}F}^{(0)}(\{p_i\})|^2}.\nonumber \\ \end{aligned}$$

(13)

The functions $C_{ca}(z;\alpha _{\mathrm {S}})$ ($c=q,{{\bar{q}}},g$) and $G_{ga}(z;\alpha _{\mathrm {S}})$ have perturbative expansions

$$\begin{aligned} C_{ca}(z;\alpha _{\mathrm {S}})= & {} \delta _{ca}\delta (1-z)+\sum _{n=1}^\infty \left( \frac{\alpha _{\mathrm {S}}}{\pi }\right) ^n C^{(n)}_{ca}(z),\nonumber \\ G_{ga}(z;\alpha _{\mathrm {S}})= & {} \sum _{n=1}^\infty \left( \frac{\alpha _{\mathrm {S}}}{\pi }\right) ^n G^{(n)}_{ga}(z) . \end{aligned}$$

(14)

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 $\widetilde{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q{{\bar{Q}}}F}$ in Eq. (9) and (13) is expressed in terms of the all-order renormalised virtual amplitude $\mathcal{M}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}$ as

$$\begin{aligned}&\!\!\!| \,\widetilde{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q {\bar{Q}}F}(\{p_i\}) \rangle \nonumber \\&\!\!\!\quad = \left[ 1 - \widetilde{\mathbf{I}}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}(\alpha _{\mathrm {S}}(M^2), \epsilon ;\{p_i\}) \right] \;| {{\mathcal {M}}}_{c{\bar{c}}\rightarrow Q {{\bar{Q}}}F}(\{p_i\}) \rangle ,\nonumber \\ \end{aligned}$$

(15)

where $\widetilde{\mathbf{I}}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}$ is the subtraction operator whose perturbative expansion is given in Eq. (4). The general expression of the first order coefficient $\widetilde{\mathbf{I}}^{(1)}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}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 ${{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}(\{ p_i \})$ on the right-hand side of Eq. (15) is the renormalised on-shell scattering amplitude and has the perturbative expansion

$$\begin{aligned} {{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}(\{ p_i \})= & {} \alpha _{\mathrm {S}}^{p/2}(\mu _R^2) \,\mu _R^{p\epsilon } \left[ {{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}^{\,(0)}(\{ p_i \}) \right. \nonumber \\&\left. + \sum _{n=1}^{\infty } \left( \frac{\alpha _{\mathrm {S}}(\mu _R^2)}{2\pi }\right) ^{\!\!n} \!\!{{\mathcal {M}}}_{c{\bar{c}}\rightarrow Q {{\bar{Q}}}F}^{\,(n)}(\{ p_i \}; \mu _R) \right] .\nonumber \\ \end{aligned}$$

(16)

The perturbative expansion of $\widetilde{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}$ is completely analogous to that in Eq. (16), with $\widetilde{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}^{\,(0)} = \mathcal{M}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}^{\,(0)}$ and the replacement $\mathcal{M}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}^{\,(n)} \rightarrow \widetilde{{\mathcal {M}}}_{c{\bar{c}}\rightarrow Q {{\bar{Q}}}F}^{\,(n)}$ ($n\ge 1$). At NLO we have

$$\begin{aligned} \widetilde{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}^{\,(1)}=\mathcal{M}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}^{\,(1)}-\widetilde{\mathbf{I}}^{(1)}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F} {{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q {\bar{Q}}F}^{\,(0)}. \end{aligned}$$

(17)

Having discussed the perturbative ingredients entering the function ${{\mathcal {H}}}^{Q{{\bar{Q}}}F}_{c{{\bar{c}}};a_1a_2}$, we can now examine its NLO and NNLO expansions. By inspecting Eqs. (9) and (12) we see that the NLO truncation of $\mathcal{H}^{Q{{\bar{Q}}}F}_{c{{\bar{c}}};a_1a_2}$ receives contributions only from the tree-level and one-loop hard-virtual amplitudes $\mathcal{M}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}^{\,(0)}$ and $\widetilde{\mathcal{M}}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}^{\,(1)}$, and from the first-order helicity-conserving coefficients $C_{ca}^{(1)}(z)$. Indeed, at this perturbative order the azimuthally dependent terms in $\mathbf{D}^{(1)}$ and $C^{\mu \nu }_{ga}$ do not contribute to Eq. (12) because of the azimuthal average. The hard-virtual one-loop amplitude can be computed with available one-loop generators such as OpenLoops [60,61,62] or Recola [63,64,65]. As a consequence, ${{\mathcal {H}}}^{Q{{\bar{Q}}}F}_{NLO}$ is available for the processes of interest, and, in particular, for $t {{\bar{t}}}H$ production.

The second-order coefficients ${{\mathcal {H}}}^{Q{{\bar{Q}}}F}_{NNLO}$ are not available in general. Indeed, they depend on the hard-virtual amplitude $\widetilde{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}^{\,(2)}$, which in turn requires the knowledge of the renormalised two-loop amplitude, and of the subtraction operator $\widetilde{\mathbf{I}}^{(2)}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}$. The computation of the two-loop amplitude for $Q{{\bar{Q}}}F$ production is at the frontier of current techniques, and, moreover, $\widetilde{\mathbf{I}}^{(2)}_{c{\bar{c}}\rightarrow Q {{\bar{Q}}}F}$ is also not known yet. However, the NNLO contribution to ${{\mathcal {H}}}^{Q{{\bar{Q}}}F}_{c{{\bar{c}}};a_1a_2}$ can be completely determined for all the flavour off-diagonal partonic channels $(a_1,a_2)\ne (c,{{\bar{c}}})$. The perturbative ingredients entering the calculation in the quark–antiquark annihilation channel (see Eq. (9)) are the corresponding tree-level (${{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}^{\,(0)}$) and one-loop ($\widetilde{{\mathcal {M}}}_{c{{\bar{c}}}\rightarrow Q {{\bar{Q}}}F}^{\,(1)}$) hard-virtual amplitudes and the first- and second-order helicity-conserving coefficients $C_{ab}^{(1)}(z)$ and $C_{ab}^{(2)}(z)$. In the gluon fusion channel, the azimuthally dependent first-order coefficients of soft ($\mathbf{D}^{(1)}$) and collinear ($G_{ga}^{(1)}$) 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 $\mathbf{D}^{(1)}({\hat{\mathbf{b}}},\{ p_i\})$ is explicitly known (see Appendix). 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\rightarrow Q{{\bar{Q}}}F$.

In summary, the current knowledge of transverse-momentum 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.

3 Results for ${\mathbf{t}{\bar{\mathbf{t}}}H}$ production

Having discussed the content of Eq. (5), we are in a position to apply it to $t {{\bar{t}}}H$ production and to obtain the complete NLO results plus the NNLO corrections in all the flavour off-diagonal partonic channels. Our NLO implementation of the calculation has the main purpose of illustrating the applicability of the $q_T$ subtraction method to $t {{\bar{t}}}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 $t {{\bar{t}}}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 $t {{\bar{t}}}H$ +jet cross section is evaluated by using the Munich code,² 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 $t {{\bar{t}}}H$ production. In both implementations all the required tree-level and one-loop amplitudes are obtained with OpenLoops [60,61,62], including the tree-level spin- and colour-correlated amplitudes required to evaluate the contributions in Eq. (12).

Fig. 1

The $r_{\mathrm{cut}}$ dependence (data points) at $\sqrt{s}=13~\mathrm {TeV}$ (left) and $100~\mathrm {TeV}$ (right) of the NLO total cross section computed by using $q_T$ subtraction. The bands show the extrapolated value at $r_{\mathrm{cut}}\rightarrow 0$ and the NLO results from Madgraph5_aMC@NLO (using FKS subtraction) and Matrix (using dipole subtraction)

In order to numerically evaluate the contribution in the square bracket of Eq. (5), a technical cut-off $r_{\mathrm{cut}}$ is introduced on the dimensionless variable $q_T / M$, where M is the invariant mass of the $t{\bar{t}} H$ system. The final result, which corresponds to the limit $r_{\mathrm{cut}}\rightarrow 0$, is extracted by computing the cross section at fixed values of $r_{\mathrm{cut}}$ in the range $[0.01\%, r_{max}]$. Quadratic least $\chi ^2$ fits are performed for different values of $r_{max}\in [0.5\%, 1\%]$. The extrapolated value is then extracted from the fit with lowest $\chi ^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]).

Table 1 The ttH 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_{\mathrm{cut}}\rightarrow 0$ extrapolation

Fig. 2

The NLO results of Madgraph5_aMC@NLO for the cross section dependence on several kinematic variables at $\sqrt{s}=13~\mathrm {TeV}$. The lower panels show the relative comparison with the corresponding results obtained by using $q_T$ subtraction

We consider pp collisions at the centre-of-mass energies $\sqrt{s}=13~\mathrm {TeV}$ and $\sqrt{s}=100~\mathrm {TeV}$. We use the NNPDF31 [72] parton distribution functions (PDFs) with the QCD running coupling $\alpha _{\mathrm {S}}$ evaluated at each corresponding order (i.e., we use $(n+1)$-loop $\alpha _{\mathrm {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\times 10^{-5}$ GeV$^{-2}$. The renormalisation and factorization scales, $\mu _R$ and $\mu _F$, are fixed at $\mu _R=\mu _F=(2m_t+m_H)/2$. Our predictions for the LO and NLO cross sections and for the NNLO corrections in the flavour off-diagonal channels are presented in Table 1 together with their uncertainties due to the numerical integration and the extrapolation to $r_{\mathrm{cut}}\rightarrow 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 $\sqrt{s}=13$ TeV ($\sqrt{s}=100$ TeV). The flavour off-diagonal $qg+{{\bar{q}}}g$ 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 Madgraph5_aMC@NLO and Matrix results. The quality of the $r_{\mathrm{cut}}\rightarrow 0$ extrapolation can be assessed by studying the behavior of the cross section at fixed values of $r_{\mathrm{cut}}$. In Fig. 1 we investigate this behavior and show also the ($r_{\mathrm{cut}}$ independent) NLO result obtained with Matrix, by using dipole subtraction, and Madgraph5_aMC@NLO, by using FKS subtraction. As expected, the $r_{\mathrm{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_{\mathrm{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 $\sqrt{s}=13~\mathrm {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 transverse-momentum (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 $t {{\bar{t}}}H$ system (bottom right). We find excellent agreement between the two calculations, with bin-wise uncertainties at the percent-level or below. Our results are obtained by using a fixed value of $r_{\mathrm{cut}}$, $r_{\mathrm{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 $\mu _R$ and $\mu _F$.

Fig. 3

The NNLO contribution $\Delta \sigma$ of the qg (top) and $q({{\bar{q}}})q^\prime$ (bottom) partonic channels to the total cross section at $\sqrt{s}=13~\mathrm {TeV}$ (left) and $100~\mathrm {TeV}$ (right). The $r_{\mathrm{cut}}$ dependence of $\Delta \sigma$ is normalized to its extrapolation $\Delta \sigma _{\mathrm{NNLO}}$ at $r_{\mathrm{cut}}\rightarrow 0$

We now move to considering the NNLO contributions to the total cross section. In Table 1 we report our results for the ${{\mathcal {O}}}(\alpha _{\mathrm {S}}^4)$ contributions to the NNLO cross section from the flavour off-diagonal partonic channels $a_1a_2\rightarrow t{{\bar{t}}}H+X$. The contribution from all the channels with $a_1a_2= qg,{{\bar{q}}}g$ is labelled by the subscript qg, and the contribution from all the channels with $a_1a_2=qq, {{\bar{q}}}{{\bar{q}}}, qq', {{\bar{q}}}{{\bar{q}}}', q{{\bar{q}}}', {{\bar{q}}} q'$ $(q\ne q')$ is labelled by the subscript $q({{\bar{q}}})q^\prime$. We see that the NNLO corrections from both contributions are very small, at the few per mille level of the NLO cross section. At $\sqrt{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 ${{\mathcal {H}}}\otimes d{{{\hat{\sigma }}}}$, which is $r_{\mathrm{cut}}$ independent, and the term in the square bracket, which depends on $r_{\mathrm{cut}}$. This cancellation is observed at both $\sqrt{s}=13$ TeV and $\sqrt{s}=100$ TeV, but it is particularly severe at $\sqrt{s}=13$ TeV, downgrading the numerical precision that can be obtained for the relative correction. Similar effects were observed for $t{{\bar{t}}}$ production in Refs. [39, 40].

In Fig. 3 we show the $r_{\mathrm{cut}}$-dependence of the NNLO corrections $\Delta \sigma$ of the flavour off-diagonal channels to the total cross section for $t {{\bar{t}}}H$ production. The result is normalised to our extrapolation $\Delta \sigma _{\mathrm{NNLO}}$ at $r_{\mathrm{cut}}\rightarrow 0$. In the qg channel at $\sqrt{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_{\mathrm{cut}}$ values show relatively large instabilities, in particular again for the qg channel at $\sqrt{s}=13$ in Fig. 3 (top-left). However, these points are not dropped in the fit, which is dominated by the behaviour at $r_{\mathrm{cut}}>0.1\%$. The $r_{\mathrm{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_{\mathrm{cut}}$ behaviour for $t{{\bar{t}}}$ production (see Fig. 1 of Ref. [40]), the behaviour in Fig. 3 for $t{{\bar{t}}}H$ production is qualitatively and quantitatively similar, and we do expect the extrapolation at $r_{\mathrm{cut}} \rightarrow 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 soft contributions and the two-loop amplitudes become available.

4 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 $Q{{\bar{Q}}}$ production in Ref. [38] can be extended to associated $Q{{\bar{Q}}}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 $Q{{\bar{Q}}}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 $t {{\bar{t}}}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.

Appendix

In this Appendix we report the explicit expressions [43] of the resummation coefficients relevant for the production of an arbitrary number of heavy quarks $Q_j$ (with colour charges $\mathbf{T}_j$, pole masses $m_j$ and transverse momenta $\mathbf{p}_{jT}$) accompanied by a colourless system F with total momentum $p_F$. More precisely, we consider the partonic process

$$\begin{aligned}&c(p_1) {{\bar{c}}}(p_2) \rightarrow Q_3(p_3) {{\bar{Q}}}_4(p_4) \dots \nonumber \\&\quad Q_{N-1}(p_{N-1}) {{\bar{Q}}}_{N}(p_N) F(p_F)\, . \end{aligned}$$

(18)

where the final-state system is produced by either massless $q{\bar{q}}$ annihilation ($c=q$) or gg fusion ($c=g$). All the momenta $p_i$ ($i \ge 1$) are in the physical region, i.e. four-momentum conservation reads $p_1+p_2=p_3+...p_N+p_F$, while colour is outgoing, such that $\sum _{i=1}^N \mathbf{T}_i=0$. In particular we have $\mathbf{T}_1^2=\mathbf{T}_2^2= C_c$, where $C_c=C_F$ if $c=q$ and $C_c=C_A$ if $c=g$. This notation corresponds to that used in Ref. [38] for the simpler case of heavy-quark pair production, $c(p_1) {{\bar{c}}}(p_2) \rightarrow Q(p_3) {{\bar{Q}}}(p_4)$. We note that the massless and heavy quarks and antiquarks of the process in Eq. (18) can have different flavours (e.g., we can have $u {{\bar{d}}} \rightarrow t {{\bar{b}}} F$).

The coefficients that are presented below refer to perturbative expansions (see Eqs. (1), (2) and (4)) in which $\alpha _{\mathrm {S}}(\mu ^2)$ denotes the renormalised QCD coupling in the ${\overline{\mathrm{MS}}}$ scheme with decoupling of the heavy quarks $Q_j$ [81]. In the case of $t {{\bar{t}}}H$ production $\alpha (\mu ^2)$ is the QCD coupling in the five-flavour scheme.

The first-order soft anomalous dimension ${\varvec{\Gamma }}_{t}^{(1)}(\{p_i\})$ (see Eq. (1)) reads

$$\begin{aligned}&{\varvec{\Gamma }}_{t}^{(1)}(\{p_i\}) = -\frac{1}{4} \left\{ \sum _{j \ge 3} \mathbf{T}_j^2 \;(1 - i \pi ) + \sum _{\begin{array}{c} i=1,2 \\ j \ge 3 \end{array}} \;\mathbf{T}_i \cdot \mathbf{T}_j \,\ln \frac{(2p_ip_j)^2}{M^2 m_j^2} \right. \nonumber \\&\quad + \left. \sum _{\begin{array}{c} j,k \ge 3 \\ j \ne k \end{array}} \;\mathbf{T}_j \cdot \mathbf{T}_k \left[ \frac{1}{2v_{jk}} \ln \left( \frac{1+v_{jk}}{1-v_{jk}}\right) - i \pi \left( \frac{1}{v_{jk}} + 1 \right) \right] \right\} , \end{aligned}$$

(19)

where

$$\begin{aligned} v_{jk}= \sqrt{1 - \frac{m_j^2 m_k^2}{(p_jp_k)^2}}\,. \end{aligned}$$

(20)

The first-order subtraction operator $\widetilde{\mathbf{I}}^{(1)}$ (see Eq. (4)) can be written as

$$\begin{aligned} \widetilde{\mathbf{I}}^{(1)}&(\epsilon ,M^2/\mu _R^2;\{p_i\})) \nonumber \\ =&- \frac{1}{2} \left( \frac{M^2}{\mu _R^2}\right) ^{-\epsilon } \left\{ \sum _{i=1,2} \left[ \left( \frac{1}{\epsilon ^2} +i\pi \frac{1}{\epsilon } -\frac{\pi ^2}{12}\right) \mathbf{T}_i^2 + \frac{1}{\epsilon } \,\gamma _i \right] \right. \nonumber \\&\left. - \frac{4}{\epsilon } \;{\varvec{\Gamma }}_{t}^{(1)}(\{p_i\}) + \mathbf{F}_t^{(1)}(\{p_i\}) \right\} , \end{aligned}$$

(21)

where the coefficients $\gamma _i$ ($i=q,{{\bar{q}}},g$) originate from collinear radiation and read $\gamma _q=\gamma _{{\bar{q}}}=3C_F/2$ and $\gamma _g= (11C_A-2N_f)/6$, $N_f$ being the number of flavours of massless quarks. The function $\mathbf{F}_t^{(1)}(\{p_i\})$ in Eq. (21) is

$$\begin{aligned} \mathbf{F}_t^{(1)}(\{p_i\})&= \sum _{j \ge 3} \mathbf{T}_j^2 \;\ln \left( \frac{m_j^2 + \mathbf{p}_{jT}^2}{m_j^2}\right) \nonumber \\&\quad - \sum _{\begin{array}{c} i=1,2 \\ j \ge 3 \end{array}} \;\mathbf{T}_i \cdot \mathbf{T}_j \; \mathrm{Li}_2\left( -\frac{\mathbf{p}_{jT}^2}{m_j^2} \right) \nonumber \\&\quad + \sum _{\begin{array}{c} j,k \ge 3 \\ j \ne k \end{array}} \;\mathbf{T}_j \cdot \mathbf{T}_k \;\frac{1}{2 v_{jk}} \,L_{jk}, \end{aligned}$$

(22)

where

$$\begin{aligned} L_{jk}&= \frac{1}{2} \ln \left( \frac{1+v_{jk}}{1-v_{jk}} \right) \, \ln \left[ \frac{(m_j^2+ \mathbf{p}_{jT}^2) (m_k^2+\mathbf{p}_{kT}^2 )}{m_j^2 \,m_k^2} \right] \nonumber \\&\quad - 2 \,\mathrm{Li}_2\left( \frac{2 v_{jk}}{1+v_{jk}}\right) - \frac{1}{4}\ln ^2\left( \frac{1+v_{jk}}{1-v_{jk}} \right) \nonumber \\&\quad + \sum _{i=1,2} \left[ \,\mathrm{Li}_2\left( 1 - \sqrt{\frac{1-v_{jk}}{1+v_{jk}}} \;r_{jk,i} \right) \right. \nonumber \\&\quad \left. + \,\mathrm{Li}_2\left( 1 - \sqrt{\frac{1-v_{jk}}{1+v_{jk}}} \;\frac{1}{r_{jk,i}} \right) + \frac{1}{2} \ln ^2 r_{jk,i} \right] \end{aligned}$$

(23)

with

$$\begin{aligned} r_{jk,i} \equiv \frac{m_k}{m_j} \frac{p_i \cdot p_j}{p_i \cdot p_k}\, . \end{aligned}$$

(24)

The function $\mathrm{Li}_2(z)$ in Eq. (23) is the customary dilogarithm function, $\mathrm{Li}_2(z) = -\!\int _0^1 \frac{dt}{t} \ln (1-zt)$.

The coefficient $\mathbf{D}^{(1)}({\hat{\mathbf{b}}}, \{p_i\})$ (see Eq. (2)) reads

$$\begin{aligned} \mathbf{D}^{(1)}({\hat{\mathbf{b}}}, \{p_i\})&=\; \sum _{j \ge 3} \mathbf{T}_j^2 \left( \frac{{\hat{\mathbf{b}}}\cdot \mathbf{p}_{jT}}{\sqrt{m_j^2 + \left( {\hat{\mathbf{b}}}\cdot \mathbf{p}_{jT}\right) ^2}}\right. \nonumber \\&\quad \times \left. \left[ \mathrm{arcsinh}\left( \frac{{\hat{\mathbf{b}}}\cdot \mathbf{p}_{jT}}{m_j}\right) + \frac{i \pi }{2}\right] - \frac{1}{2} \ln \left( \frac{m_j^2 + \mathbf{p}_{jT}^2}{m_j^2}\right) \right) \nonumber \\&\quad + \sum _{\begin{array}{c} i=1,2 \\ j \ge 3 \end{array}} \mathbf{T}_i \cdot \mathbf{T}_j \left( \mathrm{arcsinh}\left( \frac{{\hat{\mathbf{b}}}\cdot \mathbf{p}_{jT}}{m_j} \right) \right. \nonumber \\&\quad \times \left. \left[ \mathrm{arcsinh}\left( \frac{{\hat{\mathbf{b}}}\cdot \mathbf{p}_{jT}}{m_j}\right) + i \pi \, \right] \! + \frac{1}{2} \mathrm{Li}_2\!\left( -\frac{\mathbf{p}_{jT}^2}{m_j^2}\right) \! \right) \nonumber \\&\quad +\sum _{\begin{array}{c} j,k \ge 3 \\ j \ne k \end{array}} \mathbf{T}_j \cdot \mathbf{T}_k \;\frac{1}{2} \;D_{jk}, \end{aligned}$$

(25)

with the function $D_{jk}$ that is given in terms of the following one-fold integral representation:

$$\begin{aligned} D_{jk}&= \int _0^1 dx \;\frac{p_j\cdot p_k}{w_{jk}^2(x)} \left( \frac{2\,{\hat{\mathbf{b}}}\cdot \mathbf{w}_{jkT}(x)}{\sqrt{w_{jk}^2(x) + \left( {\hat{\mathbf{b}}}\cdot \mathbf{w}_{jkT}(x)\right) ^2}}\right. \nonumber \\&\quad \times \left. \left[ \mathrm{arcsinh}\left( \frac{{\hat{\mathbf{b}}}\cdot \mathbf{w}_{jkT}(x)}{\sqrt{w^2_{jk}(x)}}\right) + \frac{i \pi }{2}\right] \right. \nonumber \\&\quad \left. - \ln \left( 1+\frac{\mathbf{w}_{jkT}^2(x)}{w_{jk}^2(x)}\right) \right) , \end{aligned}$$

(26)

where the four-vector $w_{jk}^\mu (x)$ is defined as

$$\begin{aligned} w_{jk}^\mu (x) = x p_j^\mu + (1-x) p_k^\mu . \end{aligned}$$

(27)

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, $p_1+p_2=p_3+...p_N+p_F$).

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³, $c(p_1) {\bar{c}}(p_2) \rightarrow Q(p_3) {{\bar{Q}}}(p_4)$. To this purpose we can simply use the corresponding constraints $m_3=m_4, \mathbf{T}_1 + \mathbf{T}_2 = - ( \mathbf{T}_3 + \mathbf{T}_4)$ and, importantly, the constraint $\mathbf{p}_{3T} = - \mathbf{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 $\mathbf{p}_{jT}$ ($3 \le j \le N$) in Eq. (18) are independent kinematical variables (since $\mathbf{p}_{FT} \ne 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 $\mathbf{b}$-space resummation coefficient $\mathbf{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 $\phi _{jb} = \phi (\mathbf{p}_{jT}) - \phi (\mathbf{b})$. This dependence produces ensuing azimuthal correlations [38] with respect to the observable angles $\phi _{jq} = \phi (\mathbf{p}_{jT}) - \phi (\mathbf{q}_T)$ of the $\mathbf{q}_T$-differential cross section at $q_T \ne 0$. In the limit $q_T \rightarrow 0$ the azimuthal-correlation cross section is divergent order-by-order in QCD perturbation theory (transverse-momentum resummation eventually makes the cross section finite [76]).

In the case of heavy-quark pair production, $\mathbf{D}^{(1)}$ depends on $\cos ^2 \phi _{3b}= \cos ^2 \phi _{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 \phi _{3q}$ with $k=2,4,6,\ldots$) [76]. This dependence of $\mathbf{D}^{(1)}$ is due to the kinematical relation $\mathbf{p}_{3T} \simeq - \mathbf{p}_{4T}$ at $q_T \rightarrow 0$. In the case of the process in Eq. (18), the transverse momenta $\mathbf{p}_{jT}$ are kinematically independent even if $q_T \rightarrow 0$, and it turns out that $\mathbf{D}^{(1)}$ depends on both $\cos ^2 \phi _{jb}$ and $\cos \phi _{jb}$ (the dependence on $\cos \phi _{jb}$ is due to the contributions that are proportional to ‘$i \pi$’ 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 higher-order) divergent azimuthal correlations for both even and odd harmonics (i.e., harmonics with respect to $\cos ^k \phi _{3q}$ with $k=1,3,5,\dots$). 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+\mathrm{jet}$ and dijet production.

We conclude this Appendix by recalling [38] that the second-order term ${\varvec{\Gamma }}_t^{(2)}$ of the soft anomalous dimension ${\varvec{\Gamma }}_t$ (see Eq. (1)) for the process in Eq. (18) is related to the corresponding term of the anomalous dimension matrix ${\varvec{\Gamma }}(\mu )$ that controls the QCD IR divergences of scattering amplitudes with massive external particles [45,46,47,48,49]. We have

$$\begin{aligned} {\varvec{\Gamma }}_{t}(\alpha _{\mathrm {S}};&\{p_i\})=\frac{1}{2} \;{\varvec{\Gamma }}^{\,\mathrm{sub}.}(\alpha _{\mathrm {S}},\{p_i\}) - \left( \frac{\alpha _{\mathrm {S}}}{\pi }\right) ^2\nonumber \\&\times \frac{1}{4} \left( \; \left[ {\varvec{\Gamma }}_{t}^{(1)}(\{p_i\}) \,, \mathbf{F}_t^{(1)}(\{p_i\}) \right] + \pi \beta _0 \mathbf{F}_t^{(1)}(\{p_i\}) \right) \nonumber \\&+ {{\mathcal {O}}}(\alpha _{\mathrm {S}}^3), \end{aligned}$$

(28)

where $12 \pi \beta _0 = 11 C_A - 2 N_f$, ${\varvec{\Gamma }}_{t}^{(1)}(\{p_i\})$ and $\mathbf{F}_t^{(1)}(\{p_i\}$ are given in Eqs. (19) and (22), and ${\varvec{\Gamma }}^{\,\mathrm{sub.}}$ is given below. The perturbative expansion of the right-hand side of Eq. (28) includes both the first-order and second-order terms ${\varvec{\Gamma }}_{t}^{(1)}$ and ${\varvec{\Gamma }}_{t}^{(2)}$ (obviously, ${\varvec{\Gamma }}^{\,\mathrm{sub.}}= 2 (\alpha _{\mathrm {S}}/\pi ) {\varvec{\Gamma }}_{t}^{(1)} + \mathcal{O}(\alpha _{\mathrm {S}}^2)$), while terms at ${{\mathcal {O}}}(\alpha _{\mathrm {S}}^3)$ and beyond are neglected. The ‘subtracted’ anomalous dimension ${\varvec{\Gamma }}^{\,\mathrm{sub.}}(\alpha _{\mathrm {S}};\{p_i\})$ is

$$\begin{aligned} {\varvec{\Gamma }}^{\,\mathrm{sub.}}(\alpha _{\mathrm {S}},\{p_i\})= & {} {\varvec{\Gamma }}(\mu ) - \left[ \frac{1}{2} (\mathbf{T}_1^2 + \mathbf{T}_2^2) \, \gamma _{\mathrm{cusp}}(\alpha _{\mathrm {S}}) \right. \nonumber \\&\times \left. \left( \ln \frac{M^2}{\mu ^2} - i \pi \right) + 2 \gamma ^c(\alpha _{\mathrm {S}}) \right] , \end{aligned}$$

(29)

where the terms on the right-hand side are written by exactly using the notation of Eq. (5) of Ref. [48]. The term ${\varvec{\Gamma }}(\mu )$ is the anomalous-dimension matrix that controls the IR divergences of the scattering amplitude ${{\mathcal {M}}}_{c{\bar{c}} \rightarrow Q_3 {{\bar{Q}}}_4 \dots 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 ${\varvec{\Gamma }}(\mu )$ for a generic process $c{{\bar{c}}}\rightarrow F$, where the system F is colourless. The expressions of ${\varvec{\Gamma }}(\mu ), \gamma _{\mathrm{cusp}}(\alpha _{\mathrm {S}})$ and $\gamma ^c(\alpha _{\mathrm {S}})$ up to ${{\mathcal {O}}}(\alpha _{\mathrm {S}}^2)$ are explicitly given in Ref. [48].