NNLL resummation for the associated production of a top pair and a Higgs boson at the LHC

We study the resummation of soft gluon emission corrections to the production of a top-antitop pair in association with a Higgs boson at the Large Hadron Collider. Starting from a soft-gluon resummation formula derived in previous work, we develop a bespoke parton-level 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 next-to-next-to-leading 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.


Introduction
The associated production of a top-quark 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 Higgs-boson production mode in the currently ongoing Run II. The Standard Model (SM) cross section for this process at a center-of-mass energy of 13 TeV is quite small, of the order of 0.5 pb.
Differences between the measured cross section and the corresponding SM predictions could indicate the presence of new physics which modifies the top-quark 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 α 2 s α, where α s and α denote the strong coupling constant and the electromagnetic fine structure constant, respectively. The next-to-leading 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 off-shell 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 ttH 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 next-to-next-to-leading 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 ttH cross section in the production threshold limit were evaluated up to next-to-leading logarithmic (NLL) accuracy; the production threshold is defined as the kinematic region in which the partonic center-of-mass energy approaches 2m t + m H , which is the minimal energy of the final state. In [16], on the other hand, we applied Soft-Collinear Effective Theory (SCET) methods 1 in order to study the impact of soft-gluon corrections to the associated production of a top pair and a Higgs boson in the partonic threshold limit 2 , i.e. in the limit where the partonic center-of-mass energy approaches the invariant mass M of the ttH final state. The mass M is bounded from above only by the hadronic center-of-mass energy. In [16] a resummation formula for the soft emission corrections was derived and all of the elements necessary for the evaluation of that formula to next-to-next-to-leading logarithmic (NNLL) accuracy were evaluated. By using these results, a study of the approximate NNLO corrections originating from soft gluon emission in the partonic threshold limit was carried out. In particular, an in-house parton level Monte Carlo program was developed and employed to evaluate the total cross section and several differential distributions. However, a direct numerical evaluation of the soft gluon emission corrections to NNLL was not performed in [16]. Recently, results for the total cross section and invariant mass distribution at NLL accuracy in the partonic threshold limit were presented in [18].
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 ttW production [19]. In that work the resummation was carried out up to NNLL accuracy in Mellin space. An in-house 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 center-of-mass energy of 8 and 13 TeV. (The NNLL resummation in the partonic threshold limit for ttW 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 top-quark pair and a Higgs boson in Mellin space. We developed an in-house 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 center-of-mass energy of 13 TeV. Finally, Section 4 contains our conclusions.

Outline of the Calculation
The associated production of a top quark pair and a Higgs boson receives contributions from the partonic process where ij ∈ {qq,qq, gg} at lowest order in QCD, and X indicates the unobserved partonic final-state radiation. Two Mandelstam invariants play a crucial role in our discussion: The soft or partonic threshold limit is defined as the kinematic situation in which 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 In (2.4), s indicates the square of the hadronic center-of-mass energy and We use the symbol {p} to indicate the set of external momenta p 1 , · · · , p 5 . The trace Tr [H ij S ij ] is proportional to the spin and color averaged squared matrix element for ttH + X s production in the process initiated by the two partons i and j, where X s indicates the unobserved soft gluons in the final state. The hard functions H ij , which are matrices in color space, are obtained from the color decomposed virtual corrections to the 2 → 3 tree-level process. The soft functions S ij (which are also matrices in color space) are related to color-decomposed real emission corrections in the soft limit; they depend on plus distributions of the form as well as on the Dirac delta function of argument (1 − z). The parton luminosity functions f f ij are defined as the convolutions of the parton distribution functions (PDFs) for the partons i and j in the protons N 1 and N 2 : In the soft limit the indices ij ∈ {qq,qq, gg}, as at LO. The hard and soft functions are two-by-two matrices for qq-initiated (quark annihilation) processes, and three-by-three matrices for gg-initiated (gluon fusion) processes. Contributions from other production channels such asqg and qg are subleading in the soft limit. We shall refer to such processes collectively as the "quark-gluon" or the "qg" channel in what follows. The hard functions satisfy renormalization group equations governed by the soft anomalous dimension matrices Γ ij H , 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 fixed-order perturbation theory up to a given order in α s . 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 [16,[22][23][24]. In particular, the NLO hard functions were evaluated by customizing two of the one-loop 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 cross-check our results.
The resummation formula for the associated production of a ttH final state in Mellin space is similar to the one which was derived for the production of a ttW final state in [19] and reads where we introduced the Mellin transform of the luminosity functions f f ij , and (2.9) Since the soft limit z → 1 corresponds to the limit N → ∞ in Mellin space, we neglected terms suppressed by powers of 1/N in (2.8). Furthermore, in (2.9) we employed the notationN = N e γ E . The function s ij is the Mellin transform of the soft function S ij found in (2.4). The hard and soft functions in (2.8) can be evaluated in fixed order perturbation theory at scales at which they are free from large logarithms. We indicate these scales with µ h and µ s , respectively. Subsequently, by solving the renormalization group (RG) equations for the hard and soft functions one can evolve the hard scattering kernels in (2.9) to the factorization scale µ f . One obtains Large logarithmic corrections depending on the ratio of the scales µ h and µ s are resummed in the channel-dependent matrix-valued evolution factors U. The expression for the evolution factors is which is formally identical to the expression found for the corresponding quantity in carrying out the resummation for ttW production. For the definition of the various RG factors appearing in (2.11) we refer the reader to [19]. However, while for ttW production one needs to consider the evolution factor in the quark-annihilation channel only, for ttH production one also needs to evaluate the appropriate anomalous dimensions and evolution factor for the gluon fusion channel. The functions U in (2.11) depend on α s evaluated at three different scales: µ h , µ s and µ f . In practice, it is convenient to rewrite the evolution factors in terms of α s (µ h ) only. This can be done by employing the running of α s at three loops [38]. By doing this, logarithms such as ln(µ h /µ s ) appear explicitly in the formula for the evolution matrix, which becomes [19] The leading logarithmic (LL) function g 1 , the NLL function g 2 , and the NNLL function g 3 can be obtained starting from (2.11). One can see that the l.h.s of (2.10) is independent of µ h and µ s 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 µ h and µ s in any numerical evaluation of (2.11) or (2.12). The hard and soft functions are free from large logarithms if one chooses µ h ∼ M and µ s ∼ M/N . It is well known that one then faces the presence of a branch cut for large values of N in the hard scattering kernel, whose existence is related to the Landau pole in α s . In this work, we choose the integration path in the complex N 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]. Table 1. Input parameters employed throughout the calculation.

Numerical Results
In this section we present predictions for the total cross section and differential distributions for the ttH 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 state-of-the-art 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 (2.10) captures the leading terms in the Mellin-space soft limit (N → ∞) of the NLO cross section, namely the single and double powers of ln N as well as N -independent terms. Even though the N -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 ttW production in [19]. In [16], similar comparisons were also performed for ttH production, but with two differences with respect to the current work: the renormalization and factorization scales were fixed (independent of M ) instead of dynamic (correlated with M ), and the leading terms were represented in momentum space instead of Mellin space.
The NLO approximation mentioned above is easily obtained by setting µ s = µ h = µ f in the NNLL resummation formula (2.10). For this reason, the matched NLO+NNLL cross section is given by 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 re-summed perturbation theory, we will also calculate NLO+NLL results, defined as 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 µ h and µ s , in addition to that on µ f . The numerical dependence on these scales is formally of NNLL order (and is indeed canceled through µ s and µ h 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 all-orders corrections are approximated by their NNLO truncation. To this end, we consider "approximate NNLO" calculations based on the NNLL resummation formula (2.10). Approximate NNLO calculations include all powers of ln N and part of the constant terms from a complete NNLO calculation, but neglect terms which vanish as N → ∞. 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 ttW production. A detailed description of which constant pieces are included in that NNLO approximation can be found in Section 4 of [19] 3 . We match these NNLO corrections, obtained in the soft limit, with the NLO ones in the usual way: 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 N 3 LL order. We define this approximation through the matching equation 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 (3.3), which depends only on µ f by construction, the constant pieces of the NNLO expansion of the NNLL result in (3.4) contain explicit dependence on µ h and µ s , in addition to that on µ f . This scale dependence is formally of N 3 LL order, and can be used to estimate the size of such corrections to the NNLL results. Moreover, the NNLO approximation (3.4) differs from the NLO+NNLL result through terms of N 3 LO 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.

Scale choices
Numerical evaluations of the resummed formulas have a residual dependence on the choice of the hard and soft scales µ h and µ s . This feature arises from the fact that the various factors in (2.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 µ h,0 = M and µ s,0 = M/N [19,43,44]. 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 fixed-order and resummed results have a residual dependence on the factorization scale µ f . The factorization scale should be chosen in such a way that logarithms of the ratio µ f /M 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 M . Figure 1 shows the dependence of the total cross section calculated within various perturbative approximations on the choice of the ratio µ f /M at the LHC with √ s = 13 TeV. One can observe that the NLO, NLO+NLL and NLO+NNLL curves intersect each other in the vicinity of µ f /M = 0.5, while the three curves have a very different behavior for small values of µ f . In addition, Figure 1 shows that beyond-NLO corrections are quite significant for µ f /M 0.5, 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 µ f,0 = M/2 and µ f,0 = M . The choice µ f,0 = M/2 may be advantageous because the lower-order perturbative results are larger at lower µ f , 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 M 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 [µ i,0 /2, 2µ i,0 ] (i ∈ {s, f, h}). The scale uncertainty above the central value of an observable O (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 for i = f, h, s. In (3.5) κ i = µ i /µ i,0 andŌ is the value of O evaluated by setting all scales to their default values (κ i = 1 for i = f, h, s). The scale uncertainty below the central value can be obtained in the same way by combining in quadrature the quantities ∆O − i , defined as in (3.5) 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.

Total cross section
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 center-of-mass energy of 13 TeV. The results obtained are summarized in Table 2, where we set µ f,0 = M/2, in Table 3, where we set µ f,0 = M , 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  the approximate NLO results include only the leading-power contributions from the gluon fusion and quark-annihilation channels in the soft limit, the difference between these results and the NLO corrections without the qg 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 qg channel reveals that this channel contributes significantly to the scale uncertainty, in particular when one chooses µ f,0 = M/2. 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.
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  within the range predicted by the uncertainty bands of the lower-order ones. For µ f,0 = M the convergence is still reasonable but not quite as good, mainly because the NLO and NLO+NLL results are noticeably smaller than at µ f,0 = M/2. Interestingly, the NLO+NLL result has a smaller scale uncertainty than the NLO+NNLL one for µ f,0 = M , a fact which looks rather accidental considering its wider variation over a larger range of µ f , as seen in Figure 1. However, one should remember that the scales µ h and µ s 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 µ f,0 = M/2 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 µ f,0 = M , 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 higher-order 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, µ f,0 = M/2 and µ f,0 = M . 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.

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, M ; • the distribution differential with respect to the invariant mass of the top-quark pair, M tt ; • the distribution differential with respect to the transverse momentum of the Higgs boson, p H T ; • the distribution differential with respect to the transverse momentum of the top quark, p t T .  We first set the default value of the factorization scale to µ f,0 = M/2. 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 bin-by-bin ratio of the two distributions, found at the bottom of each panel, shows that the NLO  . Differential distributions at approximate NLO (blue band) compared to the NLO distributions without the quark-gluon channel contribution (red band). All settings are as in Figure 3.
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 qg 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 qg channel agree quite well and the size of the respective uncertainty bands is very similar. We remind the reader that the contribution of the qg-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 M and M tt 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 M and M tt distributions.  We conclude our discussion of the results obtained with the choice µ f,0 = M/2 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 µ i = µ i,0 for i = s, f, h. 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 N 3 LO 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 all-orders resummation. The difference between the dashed red band and the dashed black band is due to constant NNLO corrections, which are of N 3 LL 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-p H T bins). We want at this point to study results for a different choice of the default factorization scale, namely µ f,0 = M . As discussed for the case of the total cross section in Section 3.  Figure 9). Figure 8 shows that at NLO+NLL the overlap between the distributions evaluated at µ f,0 = M and µ f,0 = M/2 is not particularly good, with the band at µ f,0 = M/2 slightly larger than the one at µ f,0 = M in all bins. Figure 9 shows instead that the NLO+NNLL distributions at µ f,0 = M and µ f,0 = M/2 have a large overlap in all bins. The scale uncertainty at NLO+NNLL with µ f,0 = M is larger than the scale uncertainty at µ f,0 = M/2 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.

Conclusions
In this paper we evaluated the resummation of the soft emission corrections to the associated production of a top-quark pair and a Higgs boson at the LHC in the partonic threshold limit z → 1. 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 in-house 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) top-quark pair production cross section and in [19] for the calculation of the cross section for the associated production of a top-quark pair and a W boson.
In the previous sections we presented predictions for the total cross section for this production process at the LHC operating at a center-of-mass energy of 13 TeV. In addition, we showed results for four different differential distributions depending on the four-momenta of the massive particles in the final state: the differential distributions in the invariant mass of the ttH particles, in the invariant mass of the tt 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 µ f . In particular, for the two choices which we analyzed in detail, namely µ f,0 = M/2 and µ f,0 = M , 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 µ f,0 = M than they are at µ f,0 = M/2. However, by comparing NLO+NNLL predictions obtained by setting µ f,0 = M/2 with NLO+NNLL predictions evaluated with µ f,0 = M , 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 µ f,0 .
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 Z 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.