NNLL resummation for the associated production of a top pair and a Higgs boson at the LHC
Abstract
We study the resummation of soft gluon emission corrections to the production of a topantitop pair in association with a Higgs boson at the Large Hadron Collider. Starting from a softgluon resummation formula derived in previous work, we develop a bespoke partonlevel Monte Carlo program which can be used to calculate the total cross section along with differential distributions. We use this tool to study the phenomenological impact of the resummation to nexttonexttoleading logarithmic (NNLL) accuracy, finding that these corrections increase the total cross section and the differential distributions with respect to NLO calculations of the same observables.
TUMHEP1064/16,
IPPP/16/101
a,]Alessandro Broggio, b,c]Andrea Ferroglia, d]Ben D. Pecjak, e,f,g]and Li Lin Yang \affiliation[a]Physik Department T31, Technische Universität München, James FranckStraße 1, D85748 Garching, Germany \emailAddalessandro.broggio@tum.de \affiliation[b]Physics Department, New York City College of Technology, The City University of New York, 300 Jay Street, Brooklyn, NY 11201 USA \affiliation[c]The Graduate School and University Center, The City University of New York, 365 Fifth Avenue, New York, NY 10016 USA \emailAddaferroglia@citytech.cuny.edu \affiliation[d]Institute for Particle Physics Phenomenology, Ogden Centre for Fundamental Physics, Department of Physics, University of Durham, Science Laboratories, South Rd, Durham DH1 3LE, United Kingdom \emailAddben.pecjak@durham.ac.uk \affiliation[e]School of Physics and State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, China \affiliation[f]Collaborative Innovation Center of Quantum Matter, Beijing, China \affiliation[g]Center for High Energy Physics, Peking University, Beijing 100871, China \emailAddyanglilin@pku.edu.cn
1 Introduction
The associated production of a topquark pair and a Higgs boson can provide direct information on the Yukawa coupling of the Higgs boson to the top quark, which is crucial for verifying the origin of fermion masses and may shed light on the hierarchy problem related to the mass of the Higgs boson. For this reason, experimental collaborations at the Large Hadron Collider (LHC) are actively searching for this Higgsboson production mode in the currently ongoing Run II. The Standard Model (SM) cross section for this process at a centerofmass energy of TeV is quite small, of the order of pb.
Differences between the measured cross section and the corresponding SM predictions could indicate the presence of new physics which modifies the topquark Yukawa coupling. Consequently, a large amount of work has been devoted to the study of this process beyond leading order (LO) in the SM. The LO cross section scales as , where and denote the strong coupling constant and the electromagnetic fine structure constant, respectively. The nexttoleading order (NLO) QCD corrections to this process were first evaluated more than ten years ago [1, 2, 3, 4, 5, 6]. This process also served as a benchmark for validating automated tools for NLO calculations; in [7, 8] the NLO corrections were calculated automatically and interfaced with Monte Carlo event generators, thus including parton shower and hadronization effects. Electroweak corrections to this process were studied in [9, 10, 11]. NLO QCD and electroweak corrections were included in the POWHEG framework in [12]. In [13] the NLO corrections to the associated production of a top pair and a Higgs boson were studied by considering also the decay of the top quark and offshell effects. The cross section for the associated production of a top pair, a Higgs boson and an additional jet at NLO was evaluated in [14].
Perturbative calculations for the production process
are difficult and involved, due to the presence of five external legs, four of which carry color
charges. Consequently, it is not likely that the nexttonexttoleading order (NNLO) QCD
corrections for this process will be computed in the near
future. For this reason, the impact of soft gluon emission corrections beyond
NLO was the subject of recent studies. In [15] the
soft gluon emission corrections to the total cross
section in the production threshold limit were evaluated up to
nexttoleading logarithmic (NLL) accuracy; the production threshold
is defined as the kinematic region in which the partonic
centerofmass energy approaches , which is the minimal
energy of the final state. In [16], on the other
hand, we applied SoftCollinear Effective Theory (SCET) methods
From the technical point of view, the associated production of a top pair and a W boson shows several similarities to the associated production of a top pair and a Higgs boson. However, the former process involves only one partonic production channel in the partonic threshold limit, namely the quark annihilation channel, while the latter also receives large contributions from the gluon fusion channel. For this reason some of us recently studied the resummation of the soft gluon corrections in the partonic threshold limit to production [19]. In that work the resummation was carried out up to NNLL accuracy in Mellin space. An inhouse parton level Monte Carlo program for the numerical evaluation of the resummation formulas was developed and employed to obtain predictions for the total cross section and several differential distributions at the LHC operating at a centerofmass energy of and TeV. (The NNLL resummation in the partonic threshold limit for production in momentum space was studied in [20].)
By building upon the results of [16] and [19], in this paper we study the resummation of soft gluon emission corrections to the associated production of a topquark pair and a Higgs boson in Mellin space. We developed an inhouse parton level Monte Carlo code which allows us to evaluate numerically soft emission corrections to this process up to NNLL accuracy. By matching these results with complete NLO calculations carried out with MadGraph5_aMC@NLO
[21] (which we will indicate with MG5_aMC in the rest of this paper) we obtain predictions for the total cross section and several differential distributions which are valid to NLO+NNLL accuracy. We also compute the observables at NLO+NLL accuracy
and using NNLO approximations of the NLO+NNLL results, and show that these less precise computations
miss important effects.
The paper is organized as follows: In Section 2 we review the salient features of the technique employed to obtain and evaluate the relevant resummation formulas. In Section 3 we present predictions, valid to NLO+NNLL accuracy, for the total cross section and several differential distributions for the associated production of a top pair and a Higgs boson at the LHC operating at a centerofmass energy of TeV. Finally, Section 4 contains our conclusions.
2 Outline of the Calculation
The associated production of a top quark pair and a Higgs boson receives contributions from the partonic process
(1) 
where at lowest order in QCD, and indicates the unobserved partonic finalstate radiation. Two Mandelstam invariants play a crucial role in our discussion:
(2) 
The soft or partonic threshold limit is defined as the kinematic situation in which
(3) 
In this region, the unobserved final state can contain only soft radiation.
The factorization formula for the QCD cross section in the partonic threshold limit was derived in [16] and reads
(4) 
In (4), indicates the square of the hadronic centerofmass energy and
(5) 
We use the symbol to indicate the set of external momenta . The trace is proportional to the spin and color averaged squared matrix element for production in the process initiated by the two partons and , where indicates the unobserved soft gluons in the final state. The hard functions , which are matrices in color space, are obtained from the color decomposed virtual corrections to the treelevel process. The soft functions (which are also matrices in color space) are related to colordecomposed real emission corrections in the soft limit; they depend on plus distributions of the form
(6) 
as well as on the Dirac delta function of argument . The parton luminosity functions are defined as the convolutions of the parton distribution functions (PDFs) for the partons and in the protons and :
(7) 
In the soft limit the indices , as at LO. The hard and soft functions are twobytwo matrices for initiated (quark annihilation) processes, and threebythree matrices for initiated (gluon fusion) processes. Contributions from other production channels such as and are subleading in the soft limit. We shall refer to such processes collectively as the “quarkgluon” or the “” channel in what follows.
The hard functions satisfy renormalization group equations governed by the soft anomalous dimension matrices , which depend on the partonic channel considered. These anomalous dimension matrices, which are needed to carry out the resummation of soft gluon corrections, were derived in [22, 23]. The hard functions, soft functions, and soft anomalous dimensions must be computed in fixedorder perturbation theory up to a given order in . In this work we study the resummation up to NNLL accuracy. For this task we need to evaluate the hard functions, soft functions and soft anomalous dimensions to NLO. All of these elements were already evaluated to the order needed here [22, 23, 24, 16]. In particular, the NLO hard functions were evaluated by customizing two of the oneloop provider programs available on the market, GoSam [25, 26, 27, 28, 29] and Openloops [30]. The numerical evaluation of the hard functions for this work has been performed by using a modified version of Openloops in combination with Collier [31, 32, 33, 34, 35]. GoSam in combination with Ninja [29, 36, 37] was used to crosscheck our results.
The resummation formula for the associated production of a final state in Mellin space is similar to the one which was derived for the production of a final state in [19] and reads
(8) 
where we introduced the Mellin transform of the luminosity functions , and
(9) 
Since the soft limit corresponds to the limit in Mellin space, we neglected terms suppressed by powers of in (8). Furthermore, in (9) we employed the notation . The function is the Mellin transform of the soft function found in (4).
The hard and soft functions in (8) can be evaluated in fixed order perturbation theory at scales at which they are free from large logarithms. We indicate these scales with and , respectively. Subsequently, by solving the renormalization group (RG) equations for the hard and soft functions one can evolve the hard scattering kernels in (9) to the factorization scale . One obtains
(10) 
Large logarithmic corrections depending on the ratio of the scales and are resummed in the channeldependent matrixvalued evolution factors . The expression for the evolution factors is
(11) 
which is formally identical to the expression found for the corresponding quantity in carrying out the resummation for production. For the definition of the various RG factors appearing in (11) we refer the reader to [19]. However, while for production one needs to consider the evolution factor in the quarkannihilation channel only, for production one also needs to evaluate the appropriate anomalous dimensions and evolution factor for the gluon fusion channel.
The functions in (11) depend on evaluated at three different scales: , and . In practice, it is convenient to rewrite the evolution factors in terms of only. This can be done by employing the running of at three loops [38]. By doing this, logarithms such as appear explicitly in the formula for the evolution matrix, which becomes [19]
(12) 
with
(13) 
The leading logarithmic (LL) function , the NLL function , and the NNLL function can be obtained starting from (11). One can see that the l.h.s of (10) is independent of and if the evolution factors and the hard and soft functions are known to all orders in perturbation theory. This is impossible in practice, which introduces a residual dependence on the choice of the scales and in any numerical evaluation of (11) or (12). The hard and soft functions are free from large logarithms if one chooses and . It is well known that one then faces the presence of a branch cut for large values of in the hard scattering kernel, whose existence is related to the Landau pole in . In this work, we choose the integration path in the complex plane when evaluating the inverse Mellin transform according to the Minimal Prescription (MP) introduced in [39]. In the numerics, we need the parton luminosity functions in Mellin space. These can be constructed using techniques described in [40, 41].
3 Numerical Results
GeV  GeV  

GeV  GeV  
GeV  from MMHT 2014 PDFs 
In this section we present predictions for the total cross section and differential distributions for the production process. The main goal of this work is to obtain predictions for physical observables which are valid to NLO+NNLL accuracy. However, we also perform some systematic studies meant to provide insight into the validity of various approximations to this stateoftheart result. In all cases, we use the input parameters listed in Table 1, and MMHT 2014 PDFs [42]. We switch PDF orders as appropriate for a given perturbative approximation according to the scheme given in Table 2, where we also specify the computer code used in each case.
As a preliminary step we check that with our choice of scales and input parameters the NLO expansion of the NNLL resummation formula (which we refer to as “approximate NLO”) provides a satisfactory approximation to the exact NLO calculation. Such an approximation of (10) captures the leading terms in the Mellinspace soft limit () of the NLO cross section, namely the single and double powers of as well as independent terms. Even though the independent terms depend on the Mandelstam variables, we will refer to them as “constant” terms in what follows. Analogous comparisons of approximate NLO and complete NLO calculations were carried out for production in [19]. In [16], similar comparisons were also performed for production, but with two differences with respect to the current work: the renormalization and factorization scales were fixed (independent of ) instead of dynamic (correlated with ), and the leading terms were represented in momentum space instead of Mellin space.
The NLO approximation mentioned above is easily obtained by setting in the NNLL resummation formula (10). For this reason, the matched NLO+NNLL cross section is given by
(14) 
The difference of terms in the square brackets contributes at NNLO and beyond, adding NNLL resummation onto the NLO result. In order to study the convergence of resummed perturbation theory, we will also calculate NLO+NLL results, defined as
(15) 
The difference of terms in the square brackets contributes at NNLO and beyond, adding NLL resummation onto the NLO result. However, in contrast to the approximate NLO result, the constant piece of the NLO expansion of the NLL resummation formula contains explicit dependence on the matching scales and , in addition to that on . The numerical dependence on these scales is formally of NNLL order (and is indeed canceled through and dependence in the NLO hard and soft functions in the NNLL result), and provides an additional handle on estimating the size of NNLL corrections using the NLL resummation formula.
While we are mainly interested in NNLL resummation effects, it is also interesting to study to what extent these allorders corrections are approximated by their NNLO truncation. To this end, we consider “approximate NNLO” calculations based on the NNLL resummation formula (10). Approximate NNLO calculations include all powers of and part of the constant terms from a complete NNLO calculation, but neglect terms which vanish as . Since the constant terms are not fully determined by an NNLL calculation (only their dependence is, through the RG equations), there is some freedom as to how to construct such approximations.
Here we consider two possibilities. The first follows the procedure
used in [19] for the case of
production. A detailed description of which constant pieces are
included in that NNLO approximation can be found in Section 4 of
[19]
(16) 
where we introduced the acronym nNLO to indicate approximate NNLO corrections matched to full NLO calculations. The second NNLO approximation we consider is based on the direct expansion of the NLO+NNLL result to NNLO. This differs from the approximate NNLO result used above by constant terms, which are formally of NLL order. We define this approximation through the matching equation
(17) 
In both cases above, the difference of terms in the square brackets is a pure NNLO correction. Contrary to the approximate NNLO result used in (16), which depends only on by construction, the constant pieces of the NNLO expansion of the NNLL result in (17) contain explicit dependence on and , in addition to that on . This scale dependence is formally of NLL order, and can be used to estimate the size of such corrections to the NNLL results. Moreover, the NNLO approximation (17) differs from the NLO+NNLL result through terms of NLO and higher, so comparing the two results gives a direct measure of how important such terms are numerically. In fact, were an exact NNLO calculation to appear, adding to it these beyond NNLO terms would achieve NNLO+NNLL resummation.
3.1 Scale choices
Numerical evaluations of the resummed formulas have a residual dependence on the choice of the hard and soft scales and . This feature arises from the fact that the various factors in (10) have to be evaluated at a given order in perturbation theory. When the resummation is carried out in Mellin space the standard default choice of these scales is and [43, 44, 19]. This choice is the same one followed in the “direct QCD” resummation method [39, 45, 46], and is the one we shall use here.
Furthermore, both the fixedorder and resummed results have a residual dependence on the factorization scale . The factorization scale should be chosen in such a way that logarithms of the ratio are not large [47]. Since we are working in the partonic threshold limit it is natural to choose a dynamical value for the factorization scale which is correlated with the final state invariant mass . Figure 1 shows the dependence of the total cross section calculated within various perturbative approximations on the choice of the ratio at the LHC with TeV. One can observe that the NLO, NLO+NLL and NLO+NNLL curves intersect each other in the vicinity of , while the three curves have a very different behavior for small values of . In addition, Figure 1 shows that beyondNLO corrections are quite significant for , as the NLO result falls rather steeply away to smaller values in that region, while the other three curves remain reasonably stable.
Because of these considerations, in the following we employ two different default choices for the factorization scale, namely and . The choice may be advantageous because the lowerorder perturbative results are larger at lower , so that the apparent convergence of the perturbative series is improved, but other than this numerical fact there is no obvious reason to exclude the natural hard scale as a default choice so we study this as well. In both cases, the uncertainty associated to the choice of a default value for the scale is estimated by varying each scale in the interval (). The scale uncertainty above the central value of an observable (the total cross section, or the value of a differential cross section in a given bin) is then evaluated by combining in quadrature the quantities
(18) 
for . In (18) and is the value of evaluated by setting all scales to their default values ( for ). The scale uncertainty below the central value can be obtained in the same way by combining in quadrature the quantities , defined as in (18) but with “max” replaced by “min”. We use this procedure to obtain the perturbative uncertainties given in all of the tables and figures that follow.
3.2 Total cross section
order  PDF order  code  [fb] 
LO  LO  MG5_aMC  
app. NLO  NLO  inhouse MC  
NLO no  NLO  MG5_aMC  
NLO  NLO  MG5_aMC  
NLO+NLL  NLO  inhouse MC +MG5_aMC  
NLO+NNLL  NNLO  inhouse MC +MG5_aMC  
nNLO (Mellin)  NNLO  inhouse MC +MG5_aMC  
(NLO+NNLL)  NNLO  inhouse MC +MG5_aMC 
We begin our analysis by considering the total cross section for the associated production of a top pair and a Higgs boson at the LHC operating at a centerofmass energy of 13 TeV. The results obtained are summarized in Table 2, where we set , in Table 3, where we set , and in Figure 2, which presents a visual comparison between the main results at the two different scales.
We first compare the approximate NLO corrections generated from NNLL softgluon resummation (second row of each table), with the full NLO corrections without (third row of each table) and with (fourth row of each table) the channel. Since the approximate NLO results include only the leadingpower contributions from the gluon fusion and quarkannihilation channels in the soft limit, the difference between these results and the NLO corrections without the channel gives a measure of the importance of power corrections away from this limit. The two results are seen to differ by no more than a few percent, even though the NLO corrections are large. This shows that at NLO the power corrections away from the soft limit for these channels are quite small. Comparing the NLO results with and without the channel reveals that this channel contributes significantly to the scale uncertainty, in particular when one chooses . The fact that the leading terms in the soft limit make up the bulk of the NLO correction provides a strong motivation to resum them to all orders. No information is lost by doing this, as both sources of power corrections are taken into account by matching with NLO as discussed above. Since the power corrections are treated in fixed order, the perturbative uncertainties associated with them are estimated through the standard approach of scale variations.
order  PDF order  code  [fb] 
LO  LO  MG5_aMC  
app. NLO  NLO  inhouse MC  
NLO no  NLO  MG5_aMC  
NLO  NLO  MG5_aMC  
NLO+NLL  NLO  inhouse MC +MG5_aMC  
NLO+NNLL  NNLO  inhouse MC +MG5_aMC  
nNLO (Mellin)  NNLO  inhouse MC +MG5_aMC  
(NLO+NNLL)  NNLO  inhouse MC +MG5_aMC 
We next turn to the NLO+NLL and NLO+NNLL cross sections, which are the main results of this section. The exact numbers are given in Tables 2 and 3, and a pictorial representation is given in Figure 2. The results for the default scale choice converge quite nicely. The scale uncertainties get progressively smaller when moving from NLO to NLO+NLL to NLO+NNLL, and the higherorder results are roughly within the range predicted by the uncertainty bands of the lowerorder ones. For the convergence is still reasonable but not quite as good, mainly because the NLO and NLO+NLL results are noticeably smaller than at . Interestingly, the NLO+NLL result has a smaller scale uncertainty than the NLO+NNLL one for , a fact which looks rather accidental considering its wider variation over a larger range of , as seen in Figure 1. However, one should remember that the scales and are kept fixed at their default values in the NLO+NLL and NLO+NNLL curves of Figure 1, while they are varied as explained above in order to obtain the scale uncertainty reported in the tables.
Finally, we discuss the NNLO approximations to the NNLL resummation formula. The results in Table 2 show that for the importance of resummation effects beyond NNLO is rather small, roughly at or below the 5% level after taking scale uncertainties into account. An examination of Table 3 shows that the effects are noticeably larger at , approximately at the 10% level. In either case, Figure 2 shows very clearly that the nNLO results display an artificially small scale dependence compared to the NLO+NNLL results, confirming the cautionary statements made in [16] about the reliability of the nNLO scale dependence in estimating higherorder perturbative corrections.
The results in this section highlight the importance of an NNLL calculation. Taken as a whole, they show that both NLO+NLL and approximate NNLO calculations are a poor proxy for the more complete NLO+NNLL calculation. We have considered two default scale choices, and . However, we should emphasize that in the end the default scale choice is arbitrary, and it would not be unreasonable to combine the envelope of results from the two choices into a single, larger perturbative uncertainty. The NLO+NNLL results quoted at either scale would not change significantly through such a combination.
3.3 Differential distributions
In this section we discuss results for differential distributions. In particular, we consider:

the distribution differential with respect to the invariant mass of the top pair and Higgs boson in the final state, ;

the distribution differential with respect to the invariant mass of the topquark pair, ;

the distribution differential with respect to the transverse momentum of the Higgs boson, ;

the distribution differential with respect to the transverse momentum of the top quark, .
We first set the default value of the factorization scale to . Figure 3 shows the comparison between complete NLO calculations and approximate NLO calculations for all of the distributions listed above. We observe that for all of the distributions the approximate NLO scale uncertainty band (in blue) is included in the NLO scale uncertainty band (bins with the red frame). However, the approximate NLO uncertainty is smaller than the NLO uncertainty in all bins. Furthermore the binbybin ratio of the two distributions, found at the bottom of each panel, shows that the NLO and approximate NLO corrections have somewhat different shapes.
As for the case of the total cross section, it is reasonable to look at how the approximate NLO distributions compare to the NLO calculations when the contribution of the channel is left out from the latter. This comparison can be found in Figure 4. One can see that approximate NLO and NLO distributions without the channel agree quite well and the size of the respective uncertainty bands is very similar. We remind the reader that the contribution of the channel at NLO is included in the NLO+NLL, NLO+NNLL and nNLO predictions discussed below through the matching procedure.
The comparison between the NLO and the NLO+NNLL calculations of the differential distributions can be found in Figure 5. We see that the NLO+NNLL uncertainty band is included in the NLO scale uncertainty band in almost all bins of the distributions considered here. The exception is the bins in the far tail of the and distributions, where the NLO+NNLL band is not completely included in the NLO one, but is higher than the NLO one. In general one can observe that the central value of the NLO+NNLL calculation is slightly larger than the central value of the NLO one in almost all bins of the distributions shown in Figure 5.
Figure 6 shows a comparison between NLO+NLL and NLO+NNLL results. The central value of these two calculations is quite close in all bins. The main effect of the corrections at NLO+NNLL is to shrink slightly the scale uncertainty bands with respect to the NLO+NLL results everywhere with the exception of the bins in the far tail of the and distributions.
We conclude our discussion of the results obtained with the choice by comparing in Figure 7 the NLO+NNLL, nNLO and NLO+NNLL expanded predictions for the various distributions. The figure shows the ratio, separately for each bin, of the distribution to the NLO+NNLL prediction evaluated with for . The blue band refers to NLO+NNLL calculations, the dashed red band to nNLO calculations and the dashed black band to the NNLO expansion of the NLO+NNLL resummation. The dashed black band and the blue band thus differ by NNLL resummation effects of order NLO and higher. Numerically, these effects contribute roughly at the 5% level, and as for the total cross section the NNLO truncation of the NLO+NNLL resummation formula tends to underestimate the uncertainty of the allorders resummation. The difference between the dashed red band and the dashed black band is due to constant NNLO corrections, which are of NLL order. Taking the envelope of the two NNLO approximations (i.e. the black and red bands) gives a more realistic estimate of the scale uncertainty, which is generally contained within the NLO+NNLL result (the exception is the high bins).
We want at this point to study results for a different choice of the default factorization scale, namely . As discussed for the case of the total cross section in Section 3.2, the numerical impact of the soft emission corrections with the choice is significantly larger than the impact of the same corrections with the choice . However, NLO+NNLL predictions obtained with the two choices are in good agreement. For what concerns the differential distributions studied here this can be seen by comparing NLO+NLL calculations carried out with the choice or (Figure 8), and NLO+NNLL calculations with or (Figure 9). Figure 8 shows that at NLO+NLL the overlap between the distributions evaluated at and is not particularly good, with the band at slightly larger than the one at in all bins. Figure 9 shows instead that the NLO+NNLL distributions at and have a large overlap in all bins. The scale uncertainty at NLO+NNLL with is larger than the scale uncertainty at in all bins. The good agreement between the two bands shown in each panel of Figure 9 indicates that NLO+NNLL predictions are quite stable with respect to different (but reasonable) choices of the standard value for the factorization scale.
4 Conclusions
In this paper we evaluated the resummation of the soft emission corrections to the associated production of a topquark pair and a Higgs boson at the LHC in the partonic threshold limit . The calculation is carried out to NNLL accuracy and it is matched to the complete NLO cross section in QCD. The numerical evaluation of observables at NLO+NNLL was carried out by means of an inhouse parton level Monte Carlo code developed for this work, based on the resummation formula derived in [16]. The resummation procedure is however carried out in Mellin space, following the same approach employed in [43, 44] for the calculation of the (boosted) topquark pair production cross section and in [19] for the calculation of the cross section for the associated production of a topquark pair and a boson.
In the previous sections we presented predictions for the total cross section for this production process at the LHC operating at a centerofmass energy of TeV. In addition, we showed results for four different differential distributions depending on the fourmomenta of the massive particles in the final state: the differential distributions in the invariant mass of the particles, in the invariant mass of the pair, in the transverse momentum of the Higgs boson, and in the transverse momentum of the top quark. We found that the relative size of the NNLL corrections with respect to the NLO cross section is rather sensitive to the choice of the factorization scale . In particular, for the two choices which we analyzed in detail, namely and , it was found that the NNLL corrections enhance the total cross section and differential distributions in all bins considered. The NNLL soft emission corrections expressed as a percentage of the NLO observables are larger at than they are at . However, by comparing NLO+NNLL predictions obtained by setting with NLO+NNLL predictions evaluated with , and after accounting for the scale uncertainty affecting both predictions, we find compatible results. This fact shows that the NLO+NNLL predictions are quite stable with respect to the factorization scale choice. Indeed, it would not be unreasonable to combine the envelope of the results at the two different scale choices into a single result with a larger perturbative uncertainty, which for the case of the total cross section would be at about the 20 % level. By taking the envelope of the corresponding NLO results, one finds instead an uncertainty larger than 30 %. We also studied the total cross section and differential distributions at NLO+NLL accuracy and with NNLO approximations of the NLO+NNLL resummation formula, and found that both of these are a poor proxy for the more complete NLO+NNLL results, especially for higher values of .
The parton level Monte Carlo developed for this paper could be extended to include the decays of the top quarks and the Higgs boson following the work done in [48]. This would allow one to impose cuts on the momenta of the detected particles. Furthermore, our code could serve as a template for the calculation of the NNLL soft emission corrections to the associated production of a top pair and a boson at the LHC. The latter is a process of significant phenomenological interest which has already been investigated experimentally at both the Run I and Run II of the LHC. We plan to study the NLO+NNLL cross section for this process in future work.
Acknowledgments
The inhouse Monte Carlo code which we developed and employed to evaluate the (differential) cross sections presented in this paper was run on the computer cluster of the Center for Theoretical Physics at the Physics Department of New York City College of Technology. We thank A. Signer for discussions about the Monte Carlo implementation, P. Maierhofer and S. Pozzorini for their assistance with the program Openloops, and and G. Ossola for useful discussions. The work of A.F. is supported in part by the National Science Foundation under Grant No. PHY1417354. B.P. would like to thank the ESI Vienna for hospitality and support during the completion of this work. The work of L.L.Y. was supported in part by the National Natural Science Foundation of China under Grant No. 11575004.
Footnotes
 See [17] for an introduction to SCET.
 Often this limit is referred to as PIM kinematics. The acronym PIM stands for Pair Invariant Mass and was extensively employed in the context of topquark pair production. While the generalization to our case is trivial, the word “pair” should not be applied to the process under study here, where the final state invariant mass involves 3 particles.
 In [16] such approximate NNLO formulas were obtained starting from the resummation formula in momentum space, and thus differ from Mellin space results through power corrections and constant terms. However, we have checked that the two approaches lead to results which are numerically almost identical.
References
 W. Beenakker, S. Dittmaier, M. Kramer, B. Plumper, M. Spira and P. M. Zerwas, Higgs radiation off top quarks at the Tevatron and the LHC, Phys. Rev. Lett. 87 (2001) 201805, [hepph/0107081].
 W. Beenakker, S. Dittmaier, M. Kramer, B. Plumper, M. Spira and P. M. Zerwas, NLO QCD corrections to t antit H production in hadron collisions, Nucl. Phys. B653 (2003) 151–203, [hepph/0211352].
 L. Reina, S. Dawson and D. Wackeroth, QCD corrections to associated t antit h production at the Tevatron, Phys. Rev. D65 (2002) 053017, [hepph/0109066].
 L. Reina and S. Dawson, Nexttoleading order results for t antit h production at the Tevatron, Phys. Rev. Lett. 87 (2001) 201804, [hepph/0107101].
 S. Dawson, L. H. Orr, L. Reina and D. Wackeroth, Associated top quark Higgs boson production at the LHC, Phys. Rev. D67 (2003) 071503, [hepph/0211438].
 S. Dawson, C. Jackson, L. H. Orr, L. Reina and D. Wackeroth, Associated Higgs production with top quarks at the large hadron collider: NLO QCD corrections, Phys. Rev. D68 (2003) 034022, [hepph/0305087].
 R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, R. Pittau and P. Torrielli, Scalar and pseudoscalar Higgs production in association with a topantitop pair, Phys. Lett. B701 (2011) 427–433, [1104.5613].
 M. V. Garzelli, A. Kardos, C. G. Papadopoulos and Z. Trocsanyi, Standard Model Higgs boson production in association with a top antitop pair at NLO with parton showering, Europhys. Lett. 96 (2011) 11001, [1108.0387].
 Y. Zhang, W.G. Ma, R.Y. Zhang, C. Chen and L. Guo, QCD NLO and EW NLO corrections to production with top quark decays at hadron collider, Phys. Lett. B738 (2014) 1–5, [1407.1110].
 S. Frixione, V. Hirschi, D. Pagani, H. S. Shao and M. Zaro, Weak corrections to Higgs hadroproduction in association with a topquark pair, JHEP 09 (2014) 065, [1407.0823].
 S. Frixione, V. Hirschi, D. Pagani, H. S. Shao and M. Zaro, Electroweak and QCD corrections to toppair hadroproduction in association with heavy bosons, JHEP 06 (2015) 184, [1504.03446].
 H. B. Hartanto, B. Jager, L. Reina and D. Wackeroth, Higgs boson production in association with top quarks in the POWHEG BOX, Phys. Rev. D91 (2015) 094003, [1501.04498].
 A. Denner and R. Feger, NLO QCD corrections to offshell topantitop production with leptonic decays in association with a Higgs boson at the LHC, JHEP 11 (2015) 209, [1506.07448].
 H. van Deurzen, G. Luisoni, P. Mastrolia, E. Mirabella, G. Ossola and T. Peraro, NexttoLeadingOrder QCD Corrections to Higgs Boson Production in Association with a Top Quark Pair and a Jet, Phys. Rev. Lett. 111 (2013) 171801, [1307.8437].
 A. Kulesza, L. Motyka, T. Stebel and V. Theeuwes, Soft gluon resummation for associated production at the LHC, JHEP 03 (2016) 065, [1509.02780].
 A. Broggio, A. Ferroglia, B. D. Pecjak, A. Signer and L. L. Yang, Associated production of a top pair and a Higgs boson beyond NLO, JHEP 03 (2016) 124, [1510.01914].
 T. Becher, A. Broggio and A. Ferroglia, Introduction to SoftCollinear Effective Theory, vol. 896. Springer, 2015, 10.1007/9783319148489.
 A. Kulesza, L. Motyka, T. Stebel and V. Theeuwes, Soft gluon resummation at fixed invariant mass for associated production at the LHC, in 4th Large Hadron Collider Physics Conference (LHCP 2016) Lund, Sweden, June 1318, 2016, 2016. 1609.01619.
 A. Broggio, A. Ferroglia, G. Ossola and B. D. Pecjak, Associated production of a top pair and a W boson at nexttonexttoleading logarithmic accuracy, JHEP 09 (2016) 089, [1607.05303].
 H. T. Li, C. S. Li and S. A. Li, Renormalization group improved predictions for production at hadron colliders, Phys. Rev. D90 (2014) 094009, [1409.1460].
 J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer et al., The automated computation of treelevel and nexttoleading order differential cross sections, and their matching to parton shower simulations, JHEP 07 (2014) 079, [1405.0301].
 A. Ferroglia, M. Neubert, B. D. Pecjak and L. L. Yang, Twoloop divergences of scattering amplitudes with massive partons, Phys. Rev. Lett. 103 (2009) 201601, [0907.4791].
 A. Ferroglia, M. Neubert, B. D. Pecjak and L. L. Yang, Twoloop divergences of massive scattering amplitudes in nonabelian gauge theories, JHEP 11 (2009) 062, [0908.3676].
 V. Ahrens, A. Ferroglia, M. Neubert, B. D. Pecjak and L. L. Yang, RenormalizationGroup Improved Predictions for TopQuark Pair Production at Hadron Colliders, JHEP 09 (2010) 097, [1003.5827].
 G. Cullen, N. Greiner, G. Heinrich, G. Luisoni, P. Mastrolia, G. Ossola et al., Automated OneLoop Calculations with GoSam, Eur. Phys. J. C72 (2012) 1889, [1111.2034].
 G. Cullen et al., GS2.0: a tool for automated oneloop calculations within the Standard Model and beyond, Eur. Phys. J. C74 (2014) 3001, [1404.7096].
 T. Binoth, J. P. Guillet, G. Heinrich, E. Pilon and T. Reiter, Golem95: A Numerical program to calculate oneloop tensor integrals with up to six external legs, Comput. Phys. Commun. 180 (2009) 2317–2330, [0810.0992].
 P. Mastrolia, G. Ossola, T. Reiter and F. Tramontano, Scattering AMplitudes from Unitaritybased Reduction Algorithm at the Integrandlevel, JHEP 08 (2010) 080, [1006.0710].
 T. Peraro, Ninja: Automated Integrand Reduction via Laurent Expansion for OneLoop Amplitudes, Comput. Phys. Commun. 185 (2014) 2771–2797, [1403.1229].
 F. Cascioli, P. Maierhofer and S. Pozzorini, Scattering Amplitudes with Open Loops, Phys. Rev. Lett. 108 (2012) 111601, [1111.5206].
 A. Denner and S. Dittmaier, Reduction of one loop tensor five point integrals, Nucl. Phys. B658 (2003) 175–202, [hepph/0212259].
 A. Denner and S. Dittmaier, Reduction schemes for oneloop tensor integrals, Nucl. Phys. B734 (2006) 62–115, [hepph/0509141].
 A. Denner and S. Dittmaier, Scalar oneloop 4point integrals, Nucl. Phys. B844 (2011) 199–242, [1005.2076].
 A. Denner, S. Dittmaier and L. Hofer, COLLIER  A fortranlibrary for oneloop integrals, PoS LL2014 (2014) 071, [1407.0087].
 A. Denner, S. Dittmaier and L. Hofer, Collier: a fortranbased Complex OneLoop LIbrary in Extended Regularizations, 1604.06792.
 P. Mastrolia, E. Mirabella and T. Peraro, Integrand reduction of oneloop scattering amplitudes through Laurent series expansion, JHEP 06 (2012) 095, [1203.0291].
 H. van Deurzen, G. Luisoni, P. Mastrolia, E. Mirabella, G. Ossola and T. Peraro, Multileg Oneloop Massive Amplitudes from Integrand Reduction via Laurent Expansion, JHEP 03 (2014) 115, [1312.6678].
 S. Moch, J. A. M. Vermaseren and A. Vogt, Higherorder corrections in threshold resummation, Nucl. Phys. B726 (2005) 317–335, [hepph/0506288].
 S. Catani, M. L. Mangano, P. Nason and L. Trentadue, The Resummation of soft gluons in hadronic collisions, Nucl. Phys. B478 (1996) 273–310, [hepph/9604351].
 M. Bonvini, Resummation of soft and hard gluon radiation in perturbative QCD. PhD thesis, Genoa U., 2012. 1212.0480.
 M. Bonvini and S. Marzani, Resummed Higgs cross section at NLL, JHEP 09 (2014) 007, [1405.3654].
 L. A. HarlandLang, A. D. Martin, P. Motylinski and R. S. Thorne, Parton distributions in the LHC era: MMHT 2014 PDFs, Eur. Phys. J. C75 (2015) 204, [1412.3989].
 A. Ferroglia, B. D. Pecjak, D. J. Scott and L. L. Yang, QCD resummations for boosted top production, PoS TOP2015 (2016) 052, [1512.02535].
 B. D. Pecjak, D. J. Scott, X. Wang and L. L. Yang, Resummed differential cross sections for topquark pairs at the LHC, Phys. Rev. Lett. 116 (2016) 202001, [1601.07020].
 M. Bonvini, S. Forte, M. Ghezzi and G. Ridolfi, Threshold Resummation in SCET vs. Perturbative QCD: An Analytic Comparison, Nucl. Phys. B861 (2012) 337–360, [1201.6364].
 M. Bonvini, S. Forte, G. Ridolfi and L. Rottoli, Resummation prescriptions and ambiguities in SCET vs. direct QCD: Higgs production as a case study, JHEP 01 (2015) 046, [1409.0864].
 J. C. Collins, D. E. Soper and G. F. Sterman, Factorization of Hard Processes in QCD, Adv. Ser. Direct. High Energy Phys. 5 (1989) 1–91, [hepph/0409313].
 A. Broggio, A. S. Papanastasiou and A. Signer, Renormalizationgroup improved fully differential cross sections for top pair production, JHEP 10 (2014) 98, [1407.2532].