Associated $t \bar{t} H$ production at the LHC: theoretical predictions at NLO+NNLL accuracy

We perform threshold resummation of soft gluon corrections to the total cross section and the invariant mass distribution for the process $pp \to t\bar{t}H$. The resummation is carried out at next-to-next-to-leading-logarithmic (NNLL) accuracy using the direct QCD Mellin space technique in the three-particle invariant mass kinematics. After presenting analytical expressions we discuss the impact of resummation on the numerical predictions for the associated Higgs boson production with top quarks at the LHC. We find that NLO+NNLL resummation leads to predictions for which the central values are remarkably stable with respect to scale variation and for which theoretical uncertainties are reduced in comparison to NLO predictions.


Introduction
The measurement of Higgs boson production rates in the pp → ttH process is of the central importance to the LHC research program. The process has been intensively searched for in Run 1 [1][2][3][4][5] and its measurement is among the highest priorities of the LHC Run 2 physics programme [6][7][8]. The associated production process offers a direct way to probe the strength of the top-Higgs Yukawa coupling without making any assumptions regarding its nature. As the top-Higgs Yukawa coupling is especially sensitive to the underlying physics, ttH production provides a vital test of the Standard Model (SM) and possibly a means to probe the beyond the SM physics indirectly. It is thus highly important that precise and reliable theoretical predictions are available for this process.
For these reasons, a large amount of effort has been invested to improve theoretical description of the ttH production. While the next-to-leading-order (NLO) QCD, i.e. O(α 3 s α) predictions were obtained some time ago [9,10], they have been newly recalculated and matched to parton showers in [11][12][13][14]. In the last years, the mixed QCD-weak corrections [15] and QCD-EW corrections [16,17] of O(α 2 s α 2 ) are also available. Furthermore, the NLO EW and QCD corrections to the hadronic ttH production with off-shell top and antitop quarks have been recently obtained [18,19]. For the most part, the NLO QCD corrections are ∼ 20% at the Run 2 LHC energies, whereas the size of the (electro)weak correction is more than ten times smaller. The scale uncertainty of the NLO QCD corrections is estimated to be ∼ 10% [9,10,20].
In general, if for a given process one expects that a significant part of higher order corrections originates from emission of soft and/or collinear gluons, it is possible to improve the accuracy of theoretical predictions by employing methods of resummation. Relying on principles of factorization between various dynamical modes, they allow an all-order calculation of dominant logarithmic corrections originating from a certain kinematical limit. Supplementing fixed-order results with resummation leads not only to change in the value of the cross section but also offers a better control over the theoretical error, in particular due to cancellations of factorization scale dependence between parton distribution functions (pdfs) and the partonic cross sections. The universality of resummation concepts warrants their applications to scattering processes with arbitrary many partons in the final state [21,22], thus also to a class 2 → 3 processes and in particular the associated ttH production at the LHC.
The first step in this direction was performed by us in [23], where we presented the first calculation of the resummed total cross section for the ttH production at the next-toleading-logarithmic (NLL) accuracy. The calculation relied on application of the traditional Mellin-space resummation formalism in the absolute threshold limit, i.e. in the limit of the partonic energy √ŝ approaching the production threshold M = 2m t + m H ,ŝ → M 2 , where m t is the top quark mass and m H is the Higgs boson mass. In [23], we have achieved an all-order improvement of the theoretical predictions by taking into account a well-defined subclass of higher order corrections. However, due to the suppression of the available 3particle phase-space in the absolute production threshold limit, it is to be expected that the numerical impact of formally large logarithmic corrections resummed in these kinematics will be somewhat diminished and that contributions prevailing numerically might come from regions further away from the absolute threshold scale M . Subsequently we have performed [24] resummation of NLL corrections arising in the limit of √ŝ approaching the invariant mass threshold Q with Q 2 = (p t + pt + p H ) 2 . We have considered cross sections differential in the invariant mass Q, as well as the total cross sections obtained after integration over Q. For a 2 → 2 process, this type of resummation is often referred to as threshold resummation in the pair-invariant mass (PIM) kinematics. Threshold resummation can be also performed in the framework of the soft-collinear effective theory (SCET). The first application of this technique to a 2 → 3 process, more specifically to the process pp → ttW ± , was presented in [25]. The SCET framework was also used to obtain an approximation of the next-to-next-to-leading order (NNLO) ttH cross section and distributions [26], following from the expansion of the NNLL resummation formula. Recently, NNLL results for ttW [27], ttH [28] and ttZ [29] associated production processes appeared, based on expressing the SCET formulas in Mellin space.
In this paper, we continue the work presented in [24] and perform threshold resummation in the invariant mass limit at the NNLL accuracy using the direct QCD [30] Mellinspace approach. Compared to NLL calculations, the anomalous dimensions governing resummation need to be implemented with accuracy higher by one order. In contrast to the absolute threshold limit considered in [23], the soft anomalous dimension is a matrix in the colour space containing non-zero off-diagonal elements, thus requiring an implementation of the diagonalization procedure. We then match our NNLL cross section with the fixed-order cross section at NLO. The invariant mass kinematics also offers an opportunity to perform resummation for the differential distributions in Q.
The paper is structured as follows. In Section 2 we review threshold resummation in Mellin space, stressing the difference between the resummation in the invariant mass and the absolute threshold limits. The numerical results and their discussion is presented in Section 3, where we also compare our results to those in [28]. We summarize our most important findings in Section 4.
2 NNLL resummation in the triple invariant mass kinematics for 2 → 3 processes with two massive coloured particles in the final state The resummation of soft gluon corrections to the differential cross section dσ pp→ttH /dQ 2 is performed in Mellin space, where the Mellin moments are taken w.r.t. the variable ρ = Q 2 /S. At the partonic level, the Mellin moments for the process ij → klB, where i, j denote massless coloured partons, k, l two massive quarks and B a massive colour-singlet particle, are given by withρ = Q 2 /ŝ and {m 2 } denoting all masses entering the calculations. Taking the Mellin transform allows one to systematically treat the logarithmic terms of the form α n s [log m (1 − z)/(1 − z)] + , with m ≤ 2n − 1 and z = Q 2 /ŝ, appearing in the perturbative expansion of the partonic cross section to all orders in α s . In Mellin space these logarithms turn into logarithms of the variable N , and the threshold limit z → 1 corresponds to the limit N → ∞.
The resummed cross section in the N -space has the form [31,32] where the trace is taken over colour space. The appearance of colour dependence in Eq. (2) is inherently related to the fact that soft radiation is coherently sensitive to the colour structure of the hard process from which it is emitted. The matrix H ij→klB indicates the hard-scattering contributions, absorbing the off-shell effects, projected onto the chosen colour basis. The colour matrix S ij→klB represents the soft wide-angle emission. The functions ∆ i and ∆ j sum the logarithmic contributions due to (soft-)collinear radiation from the incoming partons. The radiative factors are thus universal for a specific initial state parton, i.e. they depend neither on the underlying colour structure nor the process. At LO the ttH production receives contributions from the qq and gg channels. We analyze the colour structure of the underlying processes in the s-channel colour bases, {c q I } and {c g I }, with The hard function H ij→klB carries no dependence on N and is given by the perturbative expansion In order to perform resummation at NLL accuracy one needs to know H (0) ij→klB , whereas NNLL accuracy requires the knowledge of the H (1) ij→klB coefficient. The soft function, on the other hand, resums logarithms of N at the rate of one power of the logarithm per power of the strong coupling. These single logarithms due to the soft emission can be confronted with double logarithms due to soft and collinear emissions resummed by the jet factors ∆ i and ∆ j . As a dimensionless function, S ij→klB depends only on the ratio of the scales. At the same time, the dependence on N enters only via Q/N [33], making S ij→klB dependent on the ratio Q/(N µ R ). The soft function is given by a solution of the renormalization group equation [31,34] and has the form whereS ij→klB plays a role of a boundary condition and is obtained by taking S ij→klB at Q 2 /(N 2 µ 2 R ) = 1 withN = N e γ E and γ E denoting the Euler constant. It is a purely eikonal function [31,34,35] and can be calculated perturbativelỹ At the lowest order the colour matrix is given by: Similarly to the hard function, knowledge of S ij→klB is required in order to perform resummation at NLL accuracy and a result for S (1) ij→klB at the NNLL accuracy. The scale of α s inS ij→kl , equal to Q 2 /N 2 , results in an order α 2 s logN term if we expandS ij→kl in α s (µ 2 R ).
The soft function evolution matricesŪ ij→klB , U ij→klB contain logarithmic enhancements due to soft wide-angle emissions [31,36] where P andP denote the path-and reverse path-ordering in the variable q, respectively. The soft anomalous dimension Γ ij→klB is a perturbative function in α s : In order to perform resummation at NLL accuracy we need to know Γ (1) ij→klB , whereas NNLL accuracy requires including Γ (2) ij→klB . The one-loop soft anomalous dimension for the process ij → klB with k, l being heavy quarks can be found e.g. in [23]. The two-loop contributions to the soft anomalous dimension were calculated in [37,38] 3 . In the tripleinvariant mass (TIM) kinematics, the soft anomalous dimension matrix in general contains off-diagonal terms, thus complicating the evaluation of the resummed cross section. At NNLL additional difficulty arises because of non-commutativity of Γ We make use of the method of [31] in order to diagonalize the one-loop soft anomalous dimension matrix. Denoting the diagonalization matrix for Γ where the diagonalized matrix is given by eigenvalues λ (1) and can be also written as Γ D . The other matrices are transformed as: At NLL accuracy, by changing the colour basis to the one in which Γ (1) ij→klB is diagonal, the path ordered exponentials in Eq. (4) reduce to sum over simple exponentials. At NNLL accuracy, to recast the path ordered exponential of the soft anomalous dimension matrix in a form containing simple exponential functions, we make use of a technique detailed in e.g. [41,42] resulting in with the subscript D indicating a diagonal matrix. The matrix K is given by where b 0 and b 1 are the first two coefficients of expansion β QCD in α s : In our calculation we set n f = 5.
In the diagonal basis of the one-loop soft anomalous dimension, up to NNLL accuracy Eq. (2) can be written as In the above equation, the H R andS R are hard and soft function matrices projected onto R colour basis. They are calculated at the NLO accuracy, i.e. including the O(α s ) terms in Eqs. (3) and (5). The LO hard matrix is derived from the Born cross section. The NLO hard matrix contains non-logarithmic contributions which are independent of N . They consist of virtual loop contributions, real terms of collinear origin and the contributions from the evolution matrices U R andŪ R , corresponding to evolution between µ R and Q, expanded up to O(α s ). The colour-decomposed virtual corrections are extracted from the calculations of the NLO cross section in the PowHel framework [13]. Aside from evolution terms, the remaining terms in H (1) R are obtained from the infrared-limit of the real corrections [43] using the method initially proposed in [44,45]. Additionally, we recalculate the one-loop soft functionS (1) [25,42]. The dependence on N in the soft functioñ S R enters only through the argument of α s in Eq. (5).
Substituting the expression for the running coupling, we obtain up to NNLL accuracy for the soft matrix evolution factors in Eq. (15) where: and The U R andŪ R factors in Eqs. (16), (17) correspond to evolution from Q/N up to Q and depend on µ R only through the argument α s . The N -independent evolution from µ R to Q is incorporated into the hard function, as noted earlier.
The other factors contributing to the resummation of logarithms, i.e. the radiative factors for the incoming partons, ∆ i and ∆ j are widely known. The results at NLL accuracy can be found for example in [39,46] and at the NNLL level in [40].
As already noted, at NLL accuracy, by changing the colour basis to R-basis, the path ordered exponentials in Eq. (4) reduce to simple exponentials. Equivalently, the NLL accuracy can be obtained by neglecting terms suppressed by a factor of α s in Eqs. (16), (17) and (18). This results in the soft matrix evolution factors turning into exponential functions for the eigenvalues of the soft anomalous dimension matrix. At NLL, it is also enough to only know the LO contributions to the hard and soft function, which results in the following expression for the resummed cross section in the Mellin space where the color indices I and J are implicitly summed over. The trace of the product of two matrices H R (0) andS (0) R returns the LO cross section.The incoming parton radiative factors ∆ i are now considered only at NLL accuracy.
In order to improve the accuracy of the numerical approximation provided by resummation, it is customary to include terms up to O(α s ) in the expansion of the hard and soft function leading to where R . We will refer to this result as "NLL w C", since the N -independent O(α s ) terms in the hard and soft function are often collected together in one function, known as the hard matching coefficient, C. Although we choose to treat these terms as in Eq. (21), we keep the name "w C" ("w" standing for "with") as a useful shorthand.
The resummation-improved cross sections for the pp → ttH process are obtained through matching the resummed expressions with the full NLO cross sections where "matched" can stand for "NLO+NNLL", "NLO+NLL" or "NLO+NLL w C" and "res" for "NNLL", "NLL" or "NLL w C", correspondingly. The inclusive cross section is obtained by integrating the invariant mass distribution given in Eq. (15) over Q 2 and represents its perturbative expansion truncated at NLO. The moments of the parton distribution functions (pdf) f i/h (x, µ 2 F ) are defined in the standard way The inverse Mellin transform (23) is evaluated numerically using a contour C in the complex-N space according to the "Minimal Prescription" method developed in Ref. [39].
3 Numerical results for the pp → ttH process at NLO+NNLL accuracy In this section we present and discuss our state-of-the-art NLO+NNLL predictions for the ttH production process at the LHC for two collision energies √ S = 13 TeV and √ S = 14 TeV. The results for the total cross section which we present below are obtained by integrating out the invariant mass distribution over invariant mass Q. The distribution in Q undergoes resummation of soft gluon corrections in the threshold limitŝ → Q 2 , i.e. in the invariant mass kinematics. This approach is different from directly resumming corrections to the total cross section in the absolute threshold limitŝ → M 2 , which we performed in [23]. Numerical results involving resummation are obtained using two independently developed in-house computer codes. Apart from NLL and NNLL predictions matched to NLO according to Eq. (22), we also show the NLL predictions supplemented with the O(α s ) non-logarithmic contributions ("NLL w C"), also matched to NLO.
In the phenomenological analysis we use m t = 173 GeV and m H = 125 GeV. The NLO cross section is calculated using the aMC@NLO code [47]. We perform the current analysis employing PDF4LHC15_100 sets [48][49][50][51][52][53] and use the corresponding values of α s . In particular, for the NLO+NLL predictions we use the NLO sets, whereas the NLO+NNLL predictions are calculated with NNLO sets. For the sake of comparison with Broggio et al. [28], we adopt the same choice of pdfs, i.e. MMHT2014 [50].
We present most of our analysis for two choices of the central values of the renormalization and factorization scales: µ 0 = µ F,0 = µ R,0 = Q and µ 0 = µ F,0 = µ R,0 = M/2. The former choice is motivated by invariant mass Q being the natural scale for the invariant mass kinematics used in resummation. The latter choice of the scale is often made in the NLO calculations of the total cross section reported in the literature, see e.g. [20]. By studying results for these two relatively distant scales, we aim to cover a span of scale choices relevant in the problem. The theoretical error due to scale variation is calculated using the so called 7-point method, where the minimum and maximum values obtained with (µ F /µ 0 , µ R /µ 0 ) = (0.5, 0.5), (0.5, 1), (1, 0.5), (1, 1), (1, 2), (2, 1), (2, 2) are considered. For reasons of technical simplicity, the pdf error is calculated for the NLO predictions, however we expect that adding the soft gluon corrections only minimally influences the value of the pdf error.
As discussed in the previous section for the evaluation of the first-order hard function matrix H IJ we need to know one-loop virtual corrections to the process, decomposed into various colour transitions IJ. We extract them numerically using the publicly available PowHel implementation of the ttH process [13]. In particular, we use analytical relations to translate between virtual corrections split into various colour configurations in the colour flow basis used in [13] and our default singlet-octet(s) bases. We cross check the consistency of results obtained in this way by comparing the colour-summed one-loop virtual contributions to Tr H (1)S(0) with the full one-loop virtual correction given by the PowHel package [13], as well as the POWHEG implementation of the ttH process [14] and the standalone MadLoop implementation in aMC@NLO [11].
We begin the discussion of numerical results by analyzing how well the full NLO result for the total cross section is approximated by the expansion of the resummed cross section up to the same accuracy in α s as in NLO. It was first pointed out in [25] in the context of the ttW production and then later in [26] for the ttH process that the qg production channel carries a relatively large numerical significance, especially in relation to the scale uncertainty. This is due to the fact that a non-zero contribution from the qg channel appears first at NLO, i.e. it is subleading w.r.t. contributions from the qq and gg channels. Correspondingly, no resummation is performed for this channel and it enters the matched resummation-improved formula Eq. (22) only through a fixed order contribution at NLO. It is then clear that in order to estimate how much of the NLO result is constituted by the terms accounted for in the resummed expression, Eq. (15), its expansion should not be directly compared with full NLO but with NLO cross section without a contribution from the qg channel. We obtain the latter result from the PowHel package [13]. Its comparison with the expansion of the resummed expression in Eq. (15) up to NLO accuracy in α s as a function of the scale µ = µ F = µ R is shown in Figure 1 for √ S = 14 TeV and two choices of the central scale µ 0 = Q and µ 0 = M/2. While in both cases the expansion of the resummed cross section differs significantly from the full NLO, the NLO result with the qg channel contribution subtracted is much better approximated by the expansion, especially for the dynamical scale choice µ = Q and for the fixed scale choice µ ≥ M/2, for the physically motivated scale choices. Such good agreement lets us conclude that the NNLL resummed formula will indeed take into account a prevailing part of the higher-order contributions from the qq and gg channels to all orders in α s . Our numerical predictions for the total cross sections at √ S = 13 TeV and √ S = 14 TeV are shown in Table 1. We report results obtained with our default scale choice µ 0 = Q as well as the fixed scale µ 0 = M/2. Additionally, we also provide results for the 'inbetween' choice of µ 0 = Q/2. While for these choices of central scale the NLO results vary by 20 %, the variation 5 reduces as the accuracy of resumation increases. In particular, the NLO+NNLL results show a remarkable stability w.r.t. the scale choice. We also observe that the 7-point method scale uncertainty of the results gets reduced with the increasing accuracy. In particular, for all scale choices, the scale uncertainty of NLO+NNLL cross section is reduced compared to the NLO scale uncertainty calculated in the same way. The degree up to which the scale uncertainty is reduced depends on the specific choice of the central scale. For example, for µ 0 = Q/2 the theoretical precision of the NLO+NNLL prediction is improved by about 40% with respect to the NLO result, bringing the scale   Table 1.
The size of the K NNLL factor measuring the impact of the higher-order logarithmic corrections, defined as the ratio of the NLO+NNLL to NLO cross sections, is shown in Table 2. It varies depending on the value of the central scale. The variation is almost entirely driven by the scale dependence of the NLO cross section. For the choice µ 0 = Q the K NNLL -factor can be as high as 1.19.
Given the conspicuous stability of the NLO+NNLL results, see Fig. 2, we are encouraged to combine our results obtained for various scale choices. For this purpose we adopt the method proposed by the Higgs Cross Section Working Group [20]. In this way, we obtain for the ttH cross section at 13 TeV σ NLO+NNLL = 500 +7.5%+3.0% −7.1%−3.0% fb, and at 14 TeV σ NLO+NNLL = 604 +7.6%+2.9% −7.1%−2.9% fb, where the first error is the theoretical uncertainty due to scale variation and the second error is the pdf uncertainty. 1.01 Table 2: Total cross section predictions at NLO+NNLL for pp → ttH at various LHC collision energies and central scale choices. The first error is the theoretical error due to scale variation calculated using the 7-point method and the second is the pdf error.
Our findings are further illustrated in the plots in Fig. 3 and Fig. 4. We show there the scale dependence of ttH total cross sections calculated with the factorization and renormalization scale kept equal, µ = µ F = µ R for two LHC collision energies √ S = 13 TeV and √ S = 14 TeV. As readily expected, apart from quantitative differences there is no visible disparity between the qualitative behaviour of results for the two energies. For the central scale choice of µ 0 = Q, we observe a steady increase in the stability of the cross section value w.r.t. scale variation as the accuracy of resummation improves from NLO+NLL to NLO+NNLL. Our final NLO+NNLL prediction is characterised by a very low scale dependence if µ F = µ R choice is made. Correspondingly, if calculated only along the µ F = µ R direction, the theoretical error on the NLO+NNLL prediction due to scale variation would be at the level of 1%, which is a significant reduction from the 10% variation of the NLO, c.f. Table 1. Results obtained with the scale choice of µ 0 = M/2 behave mostly in a similar way. Only in the very low scale regime, µ 0.2M , the NLO+NNLL cross section shows a stronger scale dependence. For this scale choice, the rise of the matched resummed predictions with the diminishing scale is driven by the fall of the expanded resummed NLL| NLO results, cf. Fig.1, and a therefore is a consequence of the relatively large scale dependence of NLO contributions stemming from the qg channel.
We further investigate the dependence on the scale but showing separately the renormalization and factorization scale dependence while keeping the other scale fixed. Fig. 5 shows the dependence on µ R and Fig. 6 on µ F for the √ S = 14 TeV. We conclude that the weak scale dependence present when the scales are varied simultaneously is a result of the opposite behaviour of the total cross section under µ F and µ R variations. The effect is similar to the cancellations between renormalization and factorization scale dependencies for threshold resummation in the absolute threshold limit which we observed in [23]. The typical decrease of the cross section with increasing µ R originates from running of α s . The behaviour under variation of the factorization scale, on the other hand, is related to the effect of scaling violation of pdfs at probed values of x. In this context, it is interesting to observe that the NLO+NLL predictions in Fig. 6 show very little µ F dependence around the central scale, in agreement with expectation of the factorization scale dependence in the resummed exponential and in the pdfs cancelling each other, here up to NLL. The relatively strong dependence on µ F of the NLO+NLL (w C) predictions can be then easily understood: the resummed expression will take into account higher-order scale-dependent terms which involve higher-order terms of both logarithmic (in N ) and non-logarithmic origin. The latter terms do not have their equivalent in the pdf evolution since the pdfs do not carry any process-specific information. Correspondingly, the µ F dependence does not cancel and can lead to strong effects if the non-logarithmic terms are numerically significant.
Given apparent cancellations between µ R and µ F scale dependence, we believe that the 7-point method of estimating the scale error, allowing for an independent variation of µ R only (for µ F = µ 0 ), is better suited here as an estimate of the theory error than the often used variation of µ = µ F = µ R . Another reason for our preference of this conservative estimate is presence of the hard and soft functions in the resummation formula, Eq. (15), which involve virtual corrections and are known only up to the order α s . Due to suppression of the LO phase space, they provide a relatively significant part of the NLO+NNLL corrections to the total cross sections, cf. Table 1. It is then justified to suppose that a similar situation might take place also at higher logarithmic orders and that the value of the yet unknown two-loop virtual corrections which feed into the secondorder coefficients in Eq. (15) can have a non-negligible impact on the predictions. With the 7-point method error estimate, we expect that this effect is included within the size of the error.
Our observation of stability of the predictions w.r.t. scale variation is also confirmed at the differential level. In Fig. 7 we show the differential distribution in the invariant mass Q of the ttH system produced at √ S = 14 TeV. While the NLO distributions calculated with µ 0 = Q and µ 0 = M/2 differ visibly, the NLO+NNLL distributions for these scale choices are very close in shape and value. The stability of the NLO+NNLL distribution w.r.t. the scale choice is demonstrated explicitly in Fig. 8. Correspondingly, the ratios of the NLO+NNLL to NLO distributions differ. In particular for the choice of µ 0 = Q the NNLL differential K-factor grows with the invariant mass and can be higher than 1.2 at large Q. The scale error for the invariant mass distribution is also calculated using the 7-point method. The error bands are slightly narrower for the NLO+NNLL distributions than at NLO. If the scale errors were calculated by variation of µ = µ F = µ R by factors of 0.5 and 2, the NLO+NNLL error bands would be considerably narrower.
We complete this part of the discussion by comparing resummed results obtained using the invariant mass kinematics with those obtained earlier by us in the absolute mass threshold limit [23]. At 13 (14) TeV, our most accurate prediction in these kinematics, i.e. the NLO+NLL cross section including the first-order hard-matching coefficient, evaluated with PDF4LHC15_100 pdf sets, amounted to σ NLO+NLL w C = 530 +7.8% −5.5% (641 +7.9% −5.5% ). The absolute mass threshold approach allows only for a fixed scale choice, which is taken to be µ 0 = µ F = µ R = M/2. Comparing this result with our NLO+NLL predictions for the same scale choice in the invariant mass kinematics, cf. Table 1, we see that the results calculated using the two resummation methods agree within errors.
In the remaining part of this section we comment on the relation of our results to the results of Broggio et al. [28]. That work relies on a resummation formula derived in the SCET framework in [26], though for the purpose of numerical calculations the Mellin space is adopted. In order to facilitate a comparison with results of [28] we recalculate our results as a function of the scale µ = µ F = µ R using MMHT2014 pdfs as in [28]. The outcome is presented in Fig. 9, where we show the NLO cross section and the matched resummed cross sections at various accuracy as a function of µ = µ F = µ R for the range of scales same as in Fig. 1 of [28]. Comparing the two figures, we find a qualitatively similar behaviour of the NLO+NNLL cross sections as a function of the scale. Likewise, we obtain bigger NNLL corrections for the µ 0 = µ F = µ R = Q scale choice than for µ 0 = µ F = µ R = Q/2. However, our NLO+NNLL results appear to be more stable wrt. the scale variation, leading to very little difference between the predictions for µ 0 = Q/2 and µ 0 = Q (cf. also Table 1). Fig. 9 additionally illustrates another feature of our results, namely that for physically relevant values of µ 0 0.3 Q the scale dependence diminishes as the accuracy of the predictions increases, independently on the choice of the central scale µ 0 .
However, it has to be noted that the scale choices made to obtain results reported in this paper and [28] are not equivalent. While our resummed expressions depend on µ F and µ R , the formulas used in [28] contain dependence on the hard and soft scales µ h and µ s , as well as µ F . The µ s scale in [28] is chosen in such a way as to mimic the scale of soft radiation in the Mellin-space framework, i.e. µ s = Q/N . Furthermore, for a given µ F the resummed central results of [28] are obtained with a fixed hard scale µ h = Q, while the exact and approximate NLO results are evaluated keeping all other scales equal to the factorization scale. There is one choice of factorization scale for which the scale setting procedure of [28] corresponds to simultaneous variation of µ = µ F = µ R , that is µ F = Q. For this choice we obtain σ NLO+NNLL = 501.7 +38.6 −34.6 fb, to be compared with 514.3 +42.9 −39.5 fb reported in [28], i.e. the central results of the two calculations agree within 2.5%. (The scale errors given together with the central values are expected to vary due to the different methods used for calculating them.) At NLO+NLL accuracy we do not find an agreement with [28]. We conclude that the differences in properties of the NLO+NNLL cross sections reported here and in [28] are likely related to handling of scale setting in the two resummation approaches.

Summary
In this work, we have investigated the impact of the soft gluon emission effects on the total cross section for the process pp → ttH at the LHC. The resummation of soft gluon emission has been performed using the Mellin-moment resummation technique at the NLO+NNLL accuracy in the three particle invariant mass kinematics. We have considered the differential distribution in the invariant mass as well as the total cross section, obtained by integrating the distribution. Our NLO+NNLL predictions are very stable with respect to a choice of the central scale µ 0 for the invariant mass distribution, and consequently also for the total cross section. As this is not the case for the NLO predictions, the NNLL corrections vary in size, depending on the choice of the scale. In general, for the energies and scale choices considered they provide a non-negative modification of the cross section, which for the scale choice of µ 0 = Q can be even higher than 20% at larger values of Q 2 .
We estimate the theoretical error due to scale variation by using the 7-point method, allowing for independent variation of renormalization and factorization scales. The overall size of the theoretical scale error becomes gradually smaller as the accuracy of resummation increases, albeit the reduction is relatively modest. The reduction would have been much more significant if the scale error had been estimated by simultaneous variation of renormalization and factorization scales, i.e. of µ = µ F = µ R . However, as it seems that the reduction in this case is a result of cancellations between factorization and renormalization scale dependencies, we choose a more conservative 7-point approach for estimating the error.
The stability of NLO+NNLL results w.r.t. the scale choice allows us to derive our best prediction for the pp → ttH total cross section at 13 TeV σ NLO+NNLL = 500 +7.5%+3.0% −7.1%−3.0% fb, and at 14 TeV σ NLO+NNLL = 604 +7.6%+2.9% −7.1%−2.9% fb, where the first error is the theoretical uncertainty due to scale variation and the second error is the pdf uncertainty. We note that the predictions are very close in their central value to the corresponding NLO predictions obtained for the scale choice µ 0 = M/2 and are compatible with them within errors, vindicating the appropriatness of this commonly made choice. However, in comparison with the NLO predictions obtained in this way, our NLO+NNLL predictions are characterized by the overall smaller size of the theory error related to scale variation. For an equivalent scale choice setup, our NLO+NNLL results for the ttH production process at the LHC agree with the results previously obtained by Broggio et al. [28].