Associated production of a top pair and a Higgs boson beyond NLO

We consider soft gluon emission corrections to the production of a top-antitop pair in association with a Higgs boson at hadron colliders. In particular, we present a soft-gluon resummation formula for this production process and gather all elements needed to evaluate it at next-to-next-to-leading logarithmic order. We employ these results to obtain approximate next-to-next-to-leading order (NNLO) formulas, and implement them in a bespoke parton-level Monte Carlo program which can be used to calculate the total cross section along with arbitrary differential distributions. We use this tool to study the phenomenological impact of the approximate NNLO corrections, finding that they increase the total cross section and the differential distributions which we evaluated in this work.


Introduction
The search for events in which a Higgs boson is produced in association with a topantitop quark pair (ttH production) will be one of the experimental goals of Run 2 of the Large Hadron Collider (LHC). While the Standard Model cross section for this process is quite small (∼ 0.6 pb at a center of mass energy of 14 TeV), it provides important information about the coupling between top quarks and Higgs bosons and, consequently, its measurement can place severe constraints on Beyond the Standard Model scenarios. It is therefore important to have precise theoretical predictions for this process within the Standard Model.
The leading-order (LO) cross section for ttH production scales as O(α 2 s α). The status of higher-order perturbative calculations is as follows. Next-to-leading-order (NLO) QCD corrections to this process were evaluated by two different groups in 2001-2003 [1][2][3][4][5][6]. The calculation of these corrections was repeated a few years ago with tools for the automated calculation of NLO corrections [7,8]; in both papers the NLO corrections were interfaced with parton showers and hadronization effects were studied. The weak corrections of O(α 2 s α 2 ) were evaluated in [9] (where also the QED corrections were considered) and [10]. An additional study of the strong and electroweak corrections to the associated production of a top-quark pair and a massive boson (Z, W or Higgs) was recently carried out in [11]. The NLO QCD and electroweak corrections have been included in the POWHEG BOX framework [12]. Recently, a study of ttH production going beyond stable top quarks was presented in [13], where differential cross sections have been computed, including the decay of the top quarks as well as offshell effects. Finally, the soft gluon emission corrections to the total ttH cross section in the production threshold limit, where the partonic center-of-mass energy approaches 2m t + m H , were evaluated up next-to-leading logarithmic (NLL) accuracy, and they were matched to NLO calculations [14].
In this paper we add to the above literature by approximating the NNLO QCD corrections to the total and differential ttH cross sections using soft-gluon resummation. In contrast to [14], we study soft-gluon corrections in the limit where the partonic center-of-mass energy approaches the invariant mass of the ttH final state, which in turn can be arbitrarily large. This limit is well-suited for the calculation of differential cross sections in addition to the total one. It is the exact analogue of the so-called pair-invariant mass (PIM) threshold limit used to study top-quark pair production at NNLL and approximate NNLO in [15], and we will use this nomenclature throughout the paper. We obtain the approximate NNLO corrections from the perturbative information contained in a soft-gluon resummation formula valid to next-to-next-toleading logarithmic (NNLL) accuracy. The derivation of this formula is based on the soft-collinear effective theory (SCET) methods 1 used to study differential top-quark pair production (tt) cross sections at NNLL in [15,17,18] (see [19][20][21][22] for SCET-based studies of the total tt cross section). In fact, since both tt and ttH production contain four colored partons, the study of soft-gluon corrections to these processes is conceptually identical and differs only because of the underlying kinematics. In particular, the soft-gluon resummation formula for both processes contains three essential ingredients, all of which are matrices in the color space needed to describe four-parton scattering: 1) a hard function, related to virtual corrections; 2) a soft function, related to real emission corrections in the soft limit; and 3) a soft anomalous dimension, which governs the structure of certain all-order soft-gluon corrections through the renormalization group (RG).
Of these three ingredients, both the NLO soft function [15,23] and NLO soft anomalous dimension [24,25] needed for NNLL resummation in processes involving two massless and two massive partons were calculated in such a way that they are valid for arbitrary kinematics (i.e. they do not use momentum conservation particular to twoto-two kinematics) and can thus be adapted directly to ttH production or indeed any tt production process in association with an additional uncolored particle. We perform such an adaptation here, and find agreement with results obtained previously for ttW production in [23]. The main technical challenge to obtaining results at NNLL accuracy is thus the calculation of the hard function to NLO, which unlike the soft function and soft anomalous dimension is process dependent. We carry out this calculation for the first time here, using a modified version of the one-loop providers GoSam [26,27], Openloops [28] and MadLoop [29]. Our result thus adds to the growing literature on hard functions for 2 → 3 processes, i.e. those obtained for ttW [23] production using MadLoop and related calculations for massless 2 → 3 scattering presented in [30]. Our procedure can be easily modified to include other 2 → 3 processes.
Our results are formally valid at NNLL for differential distributions in regions of phase space where the PIM soft limit is respected, which is guaranteed to be the case only when the partonic center-of-mass energy approaches the collider threshold energy. However, due to the mechanism of dynamical threshold enhancement [15,31], it is often the case that also observables sensitive to other regions of phase space receive their dominant contributions from soft-gluon corrections derived in the PIM threshold limit. Obvious examples would be the cross section differential in the invariant mass of the ttH final state at values far away from the machine threshold, or the total cross section obtained by integrating this distribution. Moreover, given that results in the PIM threshold limit are fully differential in the Mandelstam variables characterizing the Born process, we can equally well use them to estimate the NNLO corrections to any differential distribution which is non-vanishing at Born level.
We take advantage of this fact in the present work by implementing our results in an in-house parton level Monte Carlo, which can be used to calculate arbitrary ttH differential distributions along with the total cross section. To illustrate its use, we study approximate NNLO corrections to the p T of the Higgs, the p T of the top quark, the invariant mass of the tt pair, and the rapidities of the top quark or Higgs boson, in addition to the total cross section and differential cross section with respect to the ttH final state. By matching our NNLO approximation in the PIM threshold limit with the complete NLO calculation from MadGraph5_aMC@NLO [32], we obtain the currently most complete result for QCD corrections to differential ttH cross sections. Such a procedure is very much in the spirit of [33], and as in that work could be extended to include the effects of top-quark decays by retaining information on the spins of the final state particles.
The paper is organized as follows: in Section 2 we review the factorization properties of the partonic cross section in the soft emission limit. Furthermore, we discuss the evaluation of the various components which contribute to the approximate NNLO formulas derived in this work. In Section 3 we illustrate the structure of the approximate NNLO formulas obtained by considering the soft limit of the partonic cross section. Section 4 contains numerical calculations of the total ttH production cross section and of some differential distributions for the LHC operating at center of mass energy of 13 TeV. The calculations include the approximate NNLO formulas discussed in this work as well as the full set of NLO QCD corrections. The residual perturbative uncertainty affecting these results is discussed. Finally, we present our conclusions in Section 5.
2 Soft-gluon resummation for ttH hadroproduction We consider the partonic processes where the incoming partons i, j ∈ {q,q, g} and X is a partonic final state. Furthermore, we define the Mandelstam invariantŝ The invariant mass of the ttH final state, is of particular relevance to our work, since it enters in the definition of the soft parameter z z = M 2 s . (2.4) The PIM threshold limit (or, more simply, the soft limit) mentioned in the introduction is defined as the limit where z → 1, such that the unobserved final state X consists of soft partons only. Note that, in contrast to the production threshold limit, where the partonic center-of-mass energy approaches 2m t + m H , the PIM threshold limit does not impose constraints on the velocity of massive particles in the final state. It is therefore well suited for the study of differential cross sections. The starting point for soft-gluon resummation is the factorization of the partonic cross section in the soft limit. One then obtains the hadronic cross section for the collision process involving nucleons N 1 and N 2 at center-of-mass energy √ s by the usual convolution integral with parton distribution functions (PDFs). The form of the factorization of QCD corrections in the soft limit in the ttH case is identical to the tt one, so we can simply quote the result for the cross section in the soft limit by adapting that obtained for tt production using SCET methods in [15]. We write the result for the total cross section as (2.6) The content and notation of (2.5) is as follows. First, the object Tr[H ij S ij ] is proportional to the spin and color averaged squared matrix element for ttH + X s production through two initial-state partons with flavors i and j, where X s is an unobserved final state consisting of any number of soft gluons. The (matrix valued) hard functions H ij are related to color decomposed virtual corrections to the underlying 2 → 3 scattering process, and the (matrix valued) soft functions S ij are related to color-decomposed real emission corrections in the soft limit. To leading order in the soft limit, these soft real emission corrections receive contributions from initial-state partons with flavor indices ij ∈ {qq,qq, gg}; throughout this work we will refer to the channels involving quarks with the generic term "quark annihilation" channel, and the one involving gluons as the "gluon fusion" channel. Channels involving initial-state partons such as qg andqg are subleading in the soft limit, and shall be referred to generically as the "qg" channel. While the hard functions are simple functions of their arguments, the soft functions depend on singular (logarithmic) plus distributions of the form as well as the Dirac delta function δ(1 − z). Second, the parton luminosity function is given by where f i N is the parton distribution function for parton with flavor i in nucleon N . Finally, we write the phase-space integral for the ttH final state in the soft limit (which is identical to the Born-level phase space except that the total energy available is reduced from √ŝ to M due to soft gluon emissions) as where κ is the Kállen function (2.10) The differential of the solid angle of the Higgs boson direction in the laboratory frame is indicated by dΩ H = d cos θ H dφ H , while Ω * t is the solid angle of the top quark in the tt rest frame. The vectorsp 1 andp 2 are reduced momenta defined in such a way that (p 1 +p 2 ) 2 = M 2 .
In order to calculate binned differential cross sections using Monte-Carlo techniques we need explicit parameterizations of the four-momentap 1 . . . p 5 in terms of the integration variables in (2.9). The vectorsp 1 andp 2 in the partonic center-of-mass frame are written asp where we took the z axis to be in the direction of the incoming proton N 1 . The top and antitop vectors in the tt rest frame can be written as One can boost the top and antitop momenta in (2.12) to the partonic center-of-mass frame by using that the relative velocity between the two frames points along the direction of flight of the tt-pair in the partonic rest frame and has magnitude k * 12 E * 12 , where k * 12 and E * 12 are the magnitudes of the three-momentum and energy of the incoming parton pair in the the tt rest frame, respectively: (2.14) The four momentum of the Higgs boson p 5 can be easily written in the partonic centerof-mass frame, in which the Higgs boson recoils against the tt pair: . (2.16) For the rapidity distributions we also need the momenta of the top quarks and Higgs boson in the laboratory frame. In order to implement the required boost we use that the relative velocity between the partonic center-of-mass frame and the laboratory frame is parallel to p 1 and has magnitude ( , where x 1 is the integration variable in the definition of the luminosity (2.8), and x 2 = τ (zx 1 ), with τ and z the integration variables in (2.5).
We should emphasize that the factors of √ z in (2.5) arise by isolating and keeping exact dependence on the parton center-of-mass energy √ŝ = M √ z during two steps in the derivation of the factorized differential cross section. Both are related to the identification of the Fourier transform of the position-space soft function defined in terms of Wilson loops with the momentum-space object quoted in (2.5), and can be understood by examining Eq. (55) of [15]. The first of these is based on the observation that the Fourier transform of the position-space soft function depends on the total energy of the radiated soft partons in the center-of-mass frame, which is equal to √ŝ (1 − z). This explains the form of the first argument of the soft functions S ij , and keeping this √ z dependence is the ttH production equivalent of the PIM SCET scheme defined in [15,17] for tt production. A second factor of √ŝ appears as an overall prefactor between the position-space and Fourier-transformed functions and explains the factor of 1 √ z in the first line of (2.5). 2 Since factorization in the soft limit is valid as z → 1, we could equally well set both of these factors of √ z to unity, but we prefer to keep them as written since they appear "naturally" during the derivation and potentially account for numerically important power corrections away from the soft limit. We study numerical ambiguities due to the prescription used for these terms in Section 4.
The final perturbative ingredient needed for soft-gluon resummation is the soft anomalous dimension Γ H . We define it through the RG equation for the hard function, which reads (suppressing for the moment the dependence on the channel ij) . The hard function, soft function, and soft anomalous dimension in a given production channel all have perturbative expansions in α s . In order to perform soft-gluon resummation at NNLL, one needs their perturbative expansions to NLO. We end this section by explaining how we have extracted or calculated each of these NLO functions.
The results for the soft function and soft anomalous dimension to this order can be read off from results in the literature. The main step in the calculation of the NLO soft function is obtaining the phase-space type integrals I ij , defined as where v i is the velocity vector of the parton carrying momentum p i . These have been calculated in [23], and we have performed an independent calculation and found complete agreement with those results. Rather than collecting the explicit results here, we refer the reader to the list in Eq. (33) of [23]. Most of the notation from that equation matches ours directly, and we furthermore identify θ 3 and β 3 with and similarly for β 4 and θ 4 after obvious replacements. The position space (or, after trivial substitutions, Laplace space) soft function itself is then formed by calculating the I ij for all possible attachments of gluons to partons ij and associating to each attachment a color matrix particular to the partonic production channel. The momentum space function in (2.5) is obtained through an integral transform; all details related to the color matrices and integral transforms can be found in [15], and we shall not reproduce them here. The soft anomalous dimensions in the qq channel and gg channel were calculated to NLO in [24,25]. In our notation the results are and The perturbative expansions of all objects appearing in the soft anomalous dimensions above can be found, for instance, in the Appendix of [15]. Finally, we must determine the hard function at NLO. The definition of and procedure for calculating the hard functions in the quark annihilation and gluon fusion channels is in exact analogy with [15,34]. In a nutshell, the hard function is obtained by projecting out QCD amplitudes onto a particular color basis. The Higgs boson does not carry color charge, therefore the color bases employed for the quark annihilation and gluon fusion channels are chosen to be exactly the same as in [15]. Since the calculation described here requires the hard function up to NLO, we need to evaluate one-loop QCD amplitudes for both the quark annihilation and gluon fusion partonic processes. After UV renormalization, the one-loop amplitudes are still affected by IR divergences, which appear as poles in the limit in which the dimensional regulator ε vanishes. In order to obtain the finite amplitudes needed to build the hard functions, which are finite, one needs to subtract these residual poles. This is done by means of appropriate IR subtraction counterterms [24,25], again following the same procedure employed in [15].
For most 2 → 2 processes with a limited number of mass scales one-loop corrections can be easily evaluated analytically; this fact allowed some of us to evaluate the NLO hard functions for top-quark pair production analytically in [34]. The evaluation of the 2 → 3 amplitudes needed here is considerably more involved. However, in the last decade a number of tools for the automated numerical evaluation of multileg one-loop amplitudes became available. Most of these tools are publicly available and many rely on reduction techniques operating at the integrand level, such as the Ossola-Papadopoulos-Pittau method [35]. For this reason we decided to carry out the calculation of the NLO hard function with three of these tools: GoSam [26,27], MadLoop [29] and Openloops [28]. All of these tools required a certain level of customization in order to make the calculation of the hard function possible. 3 This approach can be easily adapted to the calculation of NLO hard functions for other processes of interest at the LHC. The calculation was tested by checking that the coefficients of the residual IR poles, evaluated numerically in several points of the phase space by means of the automated codes listed above, were correctly canceled by the appropriate IR subtraction counterterm. The finite hard functions obtained with each one of the three automated codes were then compared in a number of phase space points. Numerical agreement to more than eight digits was found in all cases.
The GoSam and Openloops codes were interfaced with an in-house Monte Carlo program which was written in order to evaluate the total cross section and the differential distributions presented in Section 4. (The program can also be easily interfaced with MadLoop if one prefers to use this particular one-loop provider.) Running time can become an issue in a Monte Carlo code, where one needs to evaluate the hard function at millions of phase-space points. The computer time needed to evaluate the hard functions in the two partonic channels is similar if one uses either MadLoop, GoSam or Openloops, provided that the reduction is carried out with CutTools [36] in Openloops. GoSam employs Ninja [37,38] as the default reduction tool, although GoSam can also be configured in such a way that this particular step of the calculation of the one loop amplitudes is done by means of Golem95 [39] or Samurai [40]. The calculation of the hard function is considerably faster if Openloops is run in combination with the (private) Fortran library Collier [41][42][43][44]. 4 The approximate NNLO predictions for the differential distributions which can be found in Section 4 were obtained by employing Openloops (in combination with Collier) and/or GoSam as providers for the hard functions.

Approximate NNLO formulas
By combining the information encoded in the NLO hard function and soft function with the solution of the RG equations that they satisfy, it is possible to resum logarithms of the ratio between the hard scale µ h (which characterizes the hard function) and the soft scale µ s (which is characteristic of the soft emission) up to NNLL accuracy. In particular, the (differential) hard-scattering kernel

1)
would like to acknowledge the very useful exchanges we had with Nicholas Greiner, Giovanni Ossola, Valentin Hirschi, and Philipp Maierhofer. 4 We are grateful to the authors of Collier for allowing us to use a binary version of their code.
can be rewritten in resummed form The definitions of the anomalous dimensions and evolution matrices, as well as the Laplace transformed soft functions found in (3.2) are the same as in [15,23]. Here we simply stress that the hard functions and soft functions are evaluated at their characteristic scale, (µ h and µ s , respectively) and are therefore free from large logarithmic corrections. As such, they can be safely evaluated up to a given order in perturbation theory. Large corrections depending on the ratio of the hard and soft scale are resummed in the evolution factors U. While the all-order hard scattering kernels do not depend on the hard and soft scales but only on the factorization scale at which the PDFs are evaluated, all practical implementations of (3.2) show a residual dependence on the hard and soft scales due to the truncation of the perturbative expansions of the hard functions, soft functions and anomalous dimensions. For this reason, when implementing resummed formulas, one must choose the hard and soft scales judiciously and carefully estimate the related theoretical uncertainty. The fixed-order expansion of the cross section and the resummation of soft emission effects are two complementary approaches to the precise determination of physical observables. For this reason, one typically wants to match resummed and fixed-order calculations in order to account for all of the known effects when obtaining phenomenological predictions. However, there are situations in which the perturbative expansion in α s is still justified, but soft gluon emission effects provide the bulk of the corrections at a given perturbative order. In those cases, one can use the resummed hard scattering kernels in order to obtain approximate formulas which include all of the terms proportional to plus distributions up to a given power of α s in fixed-order perturbation theory. To be specific, one can write where we have set µ f = µ r = µ, with µ r the renormalization scale. 5 The NNLO term in (3.3) has the following structure Note that it is possible to keep these two scales separate using the RG equations for α s .
where the P n distributions are defined as In (3.3,3.4) we dropped all arguments with the exception of µ and z. The approximate NNLO formulas for the partonic cross sections which we obtain in this work include the complete set of functions D i , some of the scale dependent terms in the function C 0 as well as partial information on the function R(z) which is non singular in the z → 1 limit. In particular, here we follow exactly the same procedure employed in [17,45]. That is, the terms included in R(z) arise from the transformation of logarithms in Laplace space back to momentum space. A complete list of those transformations for PIM kinematics can be found for example in Eq. (33) of [45]. As pointed out in [17], the C 0 term is ambiguous; in fact, in order to completely determine the coefficients multiplying the delta functions in the NNLO hard-scattering kernels, one would need to know the complete NNLO hard and soft matrices. Only the scale-dependent part of C 0 can be exactly determined, and one needs to specify which contributions are included there. One contribution to C 0 comes from the conversion of powers of Laplace-space logarithms according to Eq. (33) of [45]. Since these formula are exact, they are not a source of ambiguity for C 0 and those terms are included. Further contributions to C 0 arise from i) the product of the one-loop hard function with the one-loop soft function in Laplace space, ii) the product of the tree-level hard function with the twoloop soft function in Laplace space, and iii) the product of the two-loop hard function with the tree-level soft function in Laplace space. The contribution in i) is known exactly and therefore included while the term in ii) is unknown and dropped. One can reconstruct the scale dependent part of the contribution iii). However, it was observed in [15,17,45] that by including these extra µ-dependent terms one runs the risk of artificially reducing the scale dependence, rendering it an ineffective means of estimating theoretical uncertainties. Therefore, here again we follow [17,45] and drop completely the contributions of the two-loop hard function. The information obtained from approximate NNLO formulas can be added to the complete NLO calculation of a given observable in order to obtain what we refer to as approximate NNLO predictions for a physical quantity. The matching of the approximate NNLO calculation to complete NLO calculations is straightforward; for example, for the total cross section one finds where the subtraction of the last term avoids double counting of NLO terms proportional to plus distributions and delta functions. It must be observed that all of the terms on the r.h.s. of (3.6) must be evaluated with NNLO PDFs. To avoid lengthy superscripts, in the following we indicate matched NLO + approx. NNLO calculations with the symbol "nNLO". In contrast to resummed calculations, nNLO calculations show a residual dependence on the factorization scale only. As usual, the residual dependence of the observable on the factorization scale can be exploited in order to study and estimate the theoretical uncertainty affecting physical predictions. The use of approximate formulas offers an additional advantage: the numerical evaluations of the total cross section and distributions to approximate NNLO accuracy require shorter running times than the evaluation of the corresponding resummed formulas. For this reason, in this work we present predictions based upon approximate NNLO formulas.

Numerical analysis
In this section we present results obtained from the numerical evaluation of the nNLO formulas and discuss their implications. We cover the total cross section in Section 4.1 and differential distributions in Section 4.2.
A central issue is that the soft limit z → 1 is only guaranteed to provide accurate predictions for observables whereŝ → s, with √ s the collider energy; an example would be the case where M → √ s. More realistic observables such as the total cross section or differential distributions at their peaks are also sensitive to regions of phase space far away from z → 1. Thus, in order for corrections in the soft limit to be dominant also in those cases, the mechanism of dynamical threshold enhancement [15,31] must occur. This simply means that the parton luminosities appearing in (2.5) should drop off quickly enough away from the integration region where z → 1, that an expansion under the integrand of the partonic cross section in the soft limit is justified.
In order to address this issue we begin both of the following subsections with a comparison of approximate NLO calculations, valid in the soft limit, with the full NLO calculation. Approximate NLO calculations are obtained by re-expanding the NNLL resummed partonic cross section to NLO; consequently they reproduce completely all of the terms singular in the z → 1 limit in the NLO partonic cross section, but they miss terms which are subleading in the soft limit. We verify in all cases that the soft approximation works quite well at NLO. This obviously does not immediately imply that the same holds at higher orders, but is an important sanity check nonetheless. After these initial studies at NLO we then present the main results of this paper, namely numerical results from the NNLO approximations. We will see that these NNLO corrections tend to enhance both the total cross section and differential distributions to the top of the NLO uncertainty band, and also greatly decrease the uncertainties associated with scale variations. In fact, the residual uncertainties due to scale variation alone at approximate NNLO are so small that it is rather doubtful that they reflect the true theoretical uncertainty associated both with even higher-order soft gluon corrections and with terms subleading in the soft limit. In this section we address this issue and discuss a way to obtain a more conservative estimate of the theoretical uncertainty affecting approximate NLO and nNLO calculations.
The study that follows is meant to be illustrative rather than exhaustive. Therefore, we consider only one LHC energy, namely √ s = 13 TeV, and do not apply any cuts on the momenta of the final state particles. We carry out all of our calculations with MSTW 2008 PDFs [46], along with the additional input parameters shown in Table 1. Throughout the analysis we need exact NLO results for the total cross section and differential distributions. All of the numbers at NLO accuracy reported below are obtained from the code MadGraph5_aMC@NLO [32], which for convenience we indicate with MG5 in the rest of this work.
Finally, we conclude the introduction to this section by pointing out that in our analysis we keep the renormalization and factorization scales equal. However, we explicitly calculated the NLO cross section varying independently the renormalization and factorization scales in the range . In particular, we always find that the smallest NLO cross section is obtained by setting the renormalization and factorization scales equal to 2µ 0 , while the largest cross section is obtained by setting the two scales equal to µ 0 2. Therefore, setting the two scales equal to each other does not underestimate the theoretical uncertainty at NLO, compared to individual variations. For this reason, we feel justified in setting the renormalization and factorization scales equal also in the nNLO analysis, which greatly reduces the amount of running time required to obtain nNLO predictions. Separate variations, but unlikely to increase the final error bands we advocate in Section 4. value (µ 0 ) for the factorization scale (µ) employed in [1,2], namely

Total cross section
In addition to the complete NLO calculation, we show the NLO cross section without the contribution of the quark-gluon channel, which opens up at NLO. This channel is formally subleading in the soft emission limit and is therefore absent in approximate NLO calculations. However, it is important to keep in mind that the quark-gluon channel is accounted for by the matching procedure in nNLO calculations consider later on. Consequently, physical quantities evaluated to nNLO include the same quarkgluon channel contributions included in NLO calculations. By looking at Table 2, we observe that the approximate NLO calculation, based exclusively on the soft emission limit, captures 97.4 % of the full NLO calculation without the contribution of the quark-gluon channel. Very similar results are found for √ s = 7 TeV and √ s = 14 TeV. The residual theoretical uncertainty is estimated by evaluating the cross section also at µ = 235 2 = 117.5 GeV and at µ = 2 × 235 = 470 GeV. By looking at Table 2 one can see that the effect of the quark-gluon channel, which is not included in the last two columns on the right, is quite large on the scale variation. Furthermore, while the complete NLO cross section is a monotonically decreasing function of µ in the range [117 − 470] GeV, the NLO cross section has a maximum close to µ = 235 GeV if the quark-gluon channel is excluded. This fact explains the +0.0 in the scale variation in the third column of Table 2. This behavior is reproduced by the approximate NLO calculation (rightmost column of Table 2). A similar situation was found in the study of top-quark pair production [17]. This kind of behavior is even more pronounced if  NNLO PDFs are employed, as can be seen from Figure 1. In view of the steep decrease of all the curves in Figure 1 for values of the ratio µ (235 GeV) smaller than one, it is reasonable to choose a value for the central scale µ 0 larger than 235 GeV. 6 As an example, we choose µ 0 = 620 GeV (µ 0 (235 GeV) ∼ 2.64), which is a value close to the maximum of the distribution differential with respect to the total final state invariant mass M . We have checked that the location of this maximum is not very sensitive to the LHC energy. Table 3 shows that also when one chooses µ 0 = 620 GeV the approximate NLO calculation reproduces to a very good extent the NLO corrections if one excludes the contribution of the quark-gluon channel from the latter.
The total cross section at LO, NLO and nNLO calculated at µ 0 = 620 GeV can be found in Table 4. If one accounts for the approximate NNLO corrections, the central value of the cross section increases by about 8 % with respect to the NLO calculation, while the scale uncertainty is reduced by a factor of 5. For completeness, in Table 4 we report the LO, NLO and nNLO total cross section obtained by setting µ 0 = 235 GeV. While we do not advocate the use of µ 0 = 235 GeV for the reasons 6 The steep decrease of the cross section for small values of the factorization scale is an unphysical effect, in fact by choosing a factorization scale of the order 10 − 20 GeV one can even obtain negative values for the NLO total cross section. This effect can be cured either by incorporating resummation effects or by choosing a dynamic scale. Our goal in this work is to validate a method for the calculation of the approximate NNLO cross section, therefore we do not analyze this aspect further.    Table 5. Total cross section at the LHC with √ s = 13 TeV with an estimate of the error associated to the scale variation and to the formally subleading terms, as explained in the text. Each order is evaluated with the MSTW2008 PDFs at the corresponding perturbative order. discussed above, we observe that the range of values for the nNLO cross section obtained with µ 0 = 235 GeV, namely [483.1 − 495.6] fb, is reasonably close to the one obtained by setting µ 0 = 620 GeV, which is [472.7 − 490.1] fb. As already pointed out above for the NLO cross section, the nNLO cross section has a maximum close to µ = 235 GeV, and this explains the +0.0 in the scale variation in the third row of Table 4.
So far, we have estimated uncertainties associated with unknown corrections beyond nNLO through scale variations. The motivation for this is that such scale variations induce changes in the result which are beyond the accuracy of the nNLO calculation, that is, both beyond NNLO and also subleading soft terms even at NNLO (since the scale dependence of the approximate NNLO corrections is not exact). Indeed, this method is commonly accepted for standard fixed-order calculations. However, one might question if that method is sufficient here, given that it produces the small un-  certainty estimate observed above. The major difference compared to full fixed-order calculations is that the nNLO calculation misses subleading terms in the soft limit already at NNLO, so it is interesting to study more conservative ways of estimating their size. The most relevant of these subleading terms are next-to-leading power logarithms of the form These logarithms are singular but integrable in the threshold region. In principle an analysis of these next-to-leading-power logarithms is possible within SCET [47]. However, to date, only partial studies of these terms were completed, for the Abelian part  of the Drell-Yan process and without employing the SCET framework, see for example [48,49].
In our case we can easily evaluate the cross section using a factor of 1 z (which was the choice made in [15] and [33] for the case of tt production) rather than 1 √ z in the overall prefactor in (2.5). When expanded around the limit z → 1, each of these two choices of prefactor produces next-to-leading-power logarithms of the form (4.2), but with coefficients which differ by a factor of two. Therefore, the numerical difference between results evaluated with these two choices of prefactor gives additional insight into the generic size of such subleading terms. For both of these choices, one can then evaluate the cross section at µ = µ 0 , 2µ 0 , and µ 0 2, as usual. In this way, one obtains  six different values for the cross section and one can choose the interval between the smallest and largest value as an estimate of the residual perturbative uncertainty.
When this procedure is followed in the evaluation of the total cross section at approximate NLO one obtains the prediction found in the second column of Table 5. The central value of the approximate NLO cross section is determined by calculating the average of the maximum and minimum among the six values of the cross section. For the choice µ 0 = 620 GeV, the central value and the uncertainty interval obtained in this way are quite close to the complete NLO result shown in the first column of Table 5. While this can be somewhat accidental, it shows that, at least for the scale choice µ 0 = 620 GeV, the terms subleading in the soft limit are numerically of the same size of the quark-gluon channel contributions, which is neglected in the soft limit. The last column in Table 5 shows the nNLO total cross section calculated by estimating the residual perturbative uncertainty as it was done in the third column for the approximate NLO case. We stress once more that nNLO results are obtained by  matching the NNLO corrections in the soft limit to the complete NLO results. As such, they include the same quark-gluon channel contribution included in the NLO result. In this case the nNLO total cross section is larger than the NLO one by 5 % and the residual perturbative uncertainty is roughly half the one found at NLO. It is important to keep in mind that we do not want to attribute any special value to this way of estimating the effects of the subleading terms. The procedure is simply motivated by two goals: i) to show that scale variation alone can lead to an underestimate of the residual perturbative uncertainty affecting approximate formulas, ii) to take advantage of the fact that this procedure combined with the choice µ 0 = 620 GeV allows us to obtain approximate NLO uncertainty bands which mimic nicely the scale uncertainty bands of the complete NLO calculations, as shown in Table 5.
For the smaller choice of the reference scale, µ 0 = 235 GeV, the contribution of the quark-gluon channel to the scale uncertainty is dominant, and the method outlined above does not lead to approximate NLO predictions that mimic satisfactorily the complete NLO uncertainty band, as shown in the last row of Table 5. However, the last column of Table 5 shows that the predictions for the nNLO total cross section obtained choosing µ 0 = 235 GeV or µ 0 = 620 GeV and by subsequently estimating the residual perturbative uncertainty with the more conservative method described above are similar. In addition, as it can be seen from the last column of Table 5, the nNLO cross section range obtained starting from µ 0 = 620 GeV is larger than the one obtained starting from µ 0 = 235 GeV; we also observe that the former nearly contains the latter.

Differential distributions
An advantage of our approach is that it can be used to calculate any arbitrary differential cross section. We do this by employing standard Monte Carlo methods. In particular, during the evaluation of the approximate NNLO corrections to the total cross section in (2.5), we use the phase-space and four-momenta parameterizations described in Section 2 in order to create binned distributions. 7 In order to illustrate this approach we consider six differential distributions: • distribution differential in the invariant mass M of the ttH final state; • distribution differential in the invariant mass M tt of the tt pair; • distribution differential in the transverse momentum of the Higgs boson, p H T ; • distribution differential in the transverse momentum of the top quark, p t T ; • distribution differential in the Higgs boson rapidity, y H ; • distribution differential in the top quark rapidity, y t .
All of the distributions are evaluated in the laboratory frame. We have validated the results from our Monte Carlo based method by explicitly changing variables and calculating the first three of the distributions above by standard numerical integration in bins; the agreement between the two methods gives us confidence of the ability of our Monte Carlo implementation to calculate arbitrary distributions which are differential with respect to variables depending on the momenta of the massive particles in the final state. As with the total cross section, we begin with a comparison between the full and approximate NLO results. Figures 2 and 3 show this comparison for differential distributions in M and p H T respectively, evaluated at the default scale choice µ 0 = 620 GeV.    In addition to the full NLO results, we have also shown NLO results with the quarkgluon channel omitted. As seen from the bottom panels of the figures, the approximate NLO results recover much more than 90% of the NLO result across all bins, if one excludes from the latter the contributions of the qg channel. Even if one compares the approximate NLO distributions with the complete NLO calculation one finds that the approximate result never differs from the complete one by more than 10%. As shown in Figures 4 and 5, approximate NLO distributions satisfactorily reproduce the NLO calculations (without the quark gluon channel) also in the case in which one employs the traditional scale choice µ 0 = 235 GeV. In Figure 6 we show a comparison between the NLO result with the quark-gluon channel excluded and the approximate NLO results, this time for all six distributions listed above and including bands from scale variation as described in the caption. We see that in all cases the agreement is quite good also at values of the scale different from our default choice µ 0 = 620 GeV. Figure 7 shows the approximate NLO distributions in the case in which the uncertainty band is estimated by keeping into account the numerical effect of terms subleading in the soft limit, according to the method described in Section 4.1. The figure also shows the full (i.e. including the quark-gluon channel contribution) NLO distributions including their scale variation, evaluated with MG5. By looking at Figure 7 one can see that, as expected, the approximate NLO bands are larger than the ones in Figure 6. However, similarly to the case of the total cross section, one also observes that these larger bands at approximate NLO reproduce quite well the scale uncertainty of the complete NLO distributions. This hints to the fact that the uncertainty bands of the nNLO distributions evaluated in this way could satisfactorily mimic the scale uncertainty bands of the (unknown) full NNLO distributions.
In Figure 8 we show nNLO differential distributions along with the complete NLO results and including uncertainties from scale variation. In all cases the effects of the higher-order corrections contained in the nNLO distributions are quite similar -they enhance the results in the individual bins to the upper portion of the NLO uncertainty band, and greatly reduce the width of the bands obtained by scale variation. A more conservative estimate of the residual perturbative uncertainty affecting our predictions is shown in Figure 9. In this case the uncertainty bands are obtained by following the same procedure already employed in the calculations in Table 5 and in Figure 7. The same features observed in Figure 8 are found also in Figure 9, namely the nNLO band is located in the upper portion of the NLO uncertainty band in all the distributions which we considered. However, in Figure 9 the nNLO bands are larger than in Figure 8 and about half as large as the NLO scale variation bands, which are also shown in Figure 9.
For completeness, in Figure 10 we show the invariant mass and Higgs transverse         obtained by choosing µ 0 = 235 GeV with the ones obtained by choosing µ 0 = 620 GeV. This is done in the upper panels of Figure 11 for the invariant mass and the Higgs transverse momentum distribution. One can see that the thin nNLO scale variation bands obtained by choosing µ 0 = 235 GeV and µ 0 = 620 GeV have a large overlap, which means that the predictions for the nNLO differential distributions show little sensitivity to the choice of µ 0 . The NLO scale variation band for µ 0 = 620 GeV is shown in the background for reference. However, we stress that, as discussed above, the nNLO bands obtained by scale variation alone are likely to underestimate the residual perturbative uncertainty affecting our result. Consequently, we regard the results at µ 0 = 620 GeV with the conservative estimate of the residual perturbative uncertainty shown in Figure 9 as our best estimates for the nNLO differential distributions. In the lower panels of Figure 11 we repeat the analysis shown in the upper panels, but we show the larger uncertainty bands obtained by considering the effects of two different sets of subleading corrections, as explained above. Also in this case we find that the band corresponding to the choice µ 0 = 620 GeV has a significant overlap with the band corresponding to the choice µ 0 = 235 GeV. This fact again indicates that the nNLO predictions have little sensitivity to the choice of µ 0 . Finally, for reference, we compare the nNLO differential distributions with larger bands and µ 0 = 235 GeV to the corresponding complete NLO differential distributions in Figure 12.
We would like to conclude the discussion of the results presented here by emphasizing that the power of our approach is that it can be used to calculate arbitrary differential distributions at nNLO accuracy, a fact that was demonstrated in this section. Furthermore, cuts on the momenta of the final-state particles can easily be applied, allowing for a more direct comparison with experimental results.

Conclusions
We have studied soft-gluon corrections to the total and differential cross sections for ttH hadroproduction. The starting point was a factorization formula for the differential cross section in the PIM threshold limit z → 1, which we derived by adapting results from tt production. We then collected the perturbative ingredients needed to use this formula to perform soft-gluon resummation to NNLL order: namely the hard functions, soft functions, and soft anomalous dimensions in the quark annihilation and gluon fusion channels to NLO. While the soft functions and soft anomalous dimensions could be obtained quite easily from results in the literature, our calculation of the hard function to NLO is new. We performed and cross-checked this calculation by customizing the automated one-loop amplitude providers GoSam, MadLoop, and Openloops to extract the color-decomposed amplitudes which define the hard function. The customized programs could easily be applied to calculate the NLO hard functions for any other 2 → 3 process involving four colored partons, and thus provide an essential building block for NNLL resummations for such processes.
As a first application of our formalism to phenomenology, we studied the soft-gluon corrections in the form of approximate NNLO formulas, which are a fixed-order truncation of the resummed results. In particular, we implemented the NNLO corrections obtained from the soft-gluon resummation formalism into a bespoke parton-level Monte Carlo program, which can be used to calculate the total cross section along with arbitrary differential distributions depending on the momenta of the massive final-state particles. We illustrated the functionality of this tool in Section 4 by studying numerically the soft emission NNLO corrections to six different differential distributions at the LHC with collider energy 13 TeV, and matched these corrections with the exact NLO distributions from the event generator MadGraph5_aMC@NLO in order to obtain the above-mentioned approximate NNLO results, which we labeled nNLO. For the choice of the factorization scale to µ = µ 0 = 620 GeV, which is roughly at the peak of the distribution in the invariant mass of the ttH final state, we observed that approximate NNLO corrections enhance the (differential) cross sections compared to NLO. The same is true for results evaluated with µ = µ 0 2 and µ = 2µ 0 , although the nNLO results are much less sensitive to the choice of factorization scale than the NLO ones. This is best seen through an examination of the total cross sections at different perturbative orders listed in Table 4, and through the plots of binned distributions in Figure 8. Although kinematic cuts were not applied in our calculations, they can be implemented in the Monte-Carlo code in a straightforward way.
The current paper is a significant step forward in the study of higher-order QCD corrections in ttH production, but it is important to emphasize that there are still open issues. The overarching question is to what extent the NNLO corrections generated from the NNLL soft-gluon resummation formula approximate the true, as yet unknown NNLO corrections. We showed in Section 4 that the NLO corrections in the soft limit approximate quite well the full NLO results, and while this speaks in favor of estimating NNLO corrections using the soft limit there is no guarantee that the level of agreement seen at NLO persists at higher orders. For this reason, a conservative use of the nNLO results calculated here would require uncertainty estimates beyond the very small dependence on the factorization scale which we observed.
We addressed this is issue in Section 4. We obtained a more conservative estimate of the residual perturbative uncertainty affecting the nNLO predictions by evaluating the total and differential cross sections by including two different sets of terms which are formally subleading in the soft limit, which correspond to choose two different forms for the z dependence of the integrand in (2.5). For each of these two choices we evaluated the scale variation in the usual way. This provided us with six evaluations for each physical quantity. For the total cross section and for each bin in each distribution, the uncertainty was then determined by looking at the interval between the largest and smallest value obtained in the six calculations. This procedure, when applied to approximate NLO calculations, leads to predictions which are very close to the complete NLO predictions. When applied to nNLO this procedure leads to predictions which are slightly larger than the NLO ones both in the case of the total cross section and in the case of the six differential distributions that we considered in this work. The nNLO uncertainty interval for the total cross section and the bands for the differential distribution obtained with this method are roughly half as large as the ones obtained by evaluating these quantities at NLO. Especially for the differential distributions, it would be very useful to gain further insight into the soft-gluon corrections by implementing the true NNLL resummation instead of its NNLO approximation. We plan to return to this computationally expensive task in future work. Finally, our code has the potential to be upgraded to include the decays of the top quarks and Higgs boson, so that it could be used to evaluate observables with cuts on momenta of the detected particles.