Resummation of soft and Coulomb corrections for $t\bar{t}h$ production at the LHC

In this paper, a combined resummation of soft and Coulomb corrections is performed for the associated production of the Higgs boson with a top quark pair at the LHC. We illustrate the similarities and critical differences between this process and the $t\bar{t}$ production process. We show that up to the next-to-leading power, the total cross section for $t\bar{t}h$ production admits a similar factorization formula in the threshold limit as that for $t\bar{t}$ production. This fact, however, is not expected to hold at higher powers. Based on the factorization formula, we perform the resummation at the improved next-to-leading logarithmic accuracy, and match to the next-to-leading order result. This allows us to give NLL$'$+NLO predictions for the total cross sections at the LHC. We find that the resummation effects enhance the NLO cross sections by about 6\%, and significantly reduce the scale dependence of the theoretical predictions.


Introduction
After the discovery of the Higgs boson at the Large Hadron Collider (LHC) in 2012 [1,2], a main task of particle physics is to investigate its properties. One particularly important property is the Yukawa coupling between the top quark and the Higgs boson, which is crucial to understand the origin of the large top quark mass. Precise knowledge of the top quark Yukawa coupling will also help us to constrain new physics effects in other couplings such as the Higgs boson self-coupling [3][4][5][6]. At the LHC, the top-Higgs Yukawa coupling can be probed by measuring the cross section for Higgs boson production in association with a top quark pair (tth production). The observation of this process has been established very recently by the ATLAS [7] and CMS [8] collaborations. The measured cross sections are in agreement with theoretical predictions, although there are still large experimental uncertainties due to limited statistics. In the future, more precise measurements will be carried out at the upgraded LHC and the High Luminosity phase of the LHC (HL-LHC) [9]. For this reason, it is necessary to improve the theoretical understanding of this process.
At the leading order (LO), tth production can be initialized by qq or gg from the colliding protons, which has been studied many years ago [10][11][12]. The next-to-leading order (NLO) quantum chromodynamics (QCD) corrections were calculated in [13][14][15][16][17][18], while the NLO electroweak contributions were calculated in [19][20][21]. Given the complicated structure of this scattering process, it will be very challenging to perform an exact calculation of the next-to-next-to-leading order corrections. Therefore, a lot of efforts have been devoted to approximations of the higher order QCD corrections in various kinematic limits [22][23][24][25]. The derivation of an approximation often involves a factorization formula, which can be used to resum a class of large corrections to all orders in perturbation theory. For example, Refs. [23][24][25] have investigated the limitŝ → M 2 tth , where √ŝ is the partonic center-of-mass energy and M tth is the invariant mass of the tth system. In this limit, large logarithmic corrections are present at each order in perturbation theory. These corrections are resummed to all orders in the strong coupling α s up to the next-to-next-to-leading logarithmic (NNLL) accuracy [24,25]. In this paper, we consider a different kinematic limit s → 2m t + m h , where m t and m h are the masses of the top quark and the Higgs boson, respectively.
The limit under consideration is a bit different from the limit taken in [24,25]. There are power-like corrections arising from exchanges of Coulomb gluons, besides logarithmic corrections coming from soft gluon emissions. This limit has been studied in [22], where the soft gluon contributions are resummed, but the Coulomb gluon contributions are only incorporated at fixed-order. In this paper, we will derive a factorization formula which can resum simultaneously both kinds of higher-order corrections. The framework presented in this paper closely resembles that in tt production [26][27][28]. However, it should be emphasized that tth production is more complicated than tt production. In particular, in the threshold limit, the tt pair has a non-vanishing transverse momentum given by the recoil against the extra Higgs boson. This will lead to more involved interplay between the soft and Coulomb gluons, as will be clear later. The derivation of the factorization formula therefore requires new analyses other than those in [26,27]. It is also much more difficult to calculate the hard function describing the contributions from hard gluons with typical momentum scale around √ŝ . The paper is organized as follows. In Sec. 2, we use the analytic form of the LO partonic cross sections to analyze their behavior in the threshold limit. In Sec. 3, we present the derivation of the factorization and resummation formulas using effective field theory methods. In Sec. 4, we show the numeric results based on our resummation formula. We conclude in Sec. 5.

Analyses of leading order results
The total cross section for inclusive tth production at hadron colliders can be expressed as [29] σ(s, m t , m h ) = i,j The above factorization formula involves the partonic cross sectionσ ij and the effective parton luminosity function ff ij . The definition of the latter is where f i/N is the non-perturbative parton distribution function (PDF) of the parton i in the hadron N . It is universal and can be extracted from experimental data. The partonic cross sectionσ ij can be calculated in perturbative QCD. In this work, we are interested in its behavior near the threshold limit τ τ min , or √ŝ 2m t + m h . Here √ŝ is the partonic center-of-mass energy defined byŝ ≡ τ s. To analyze this region, we define β = 1 − (2m t + m h ) 2 /ŝ. The parameter β goes to zero in the threshold limit, and represents the typical momenta of final state particles. To see this, consider the energy conservation condition in the partonic center-of-mass frame where p t(t) and p h are the 3-momenta of the (anti-)top quark and the Higgs boson, respectively, and E X is the total energy of other emitted particles in the final state. In the limit β → 0, the 3-momenta of the (anti-)top quark and the Higgs boson becomes much smaller than their rest mass. The right side of the above equation can then be expanded in the small momenta and we obtain the following relation: We therefore have the power-counting √ŝ This will be important for establishing the effective field theory description later in the next section. For the moment, we are going to investigate the behavior of the partonic cross sections in the limit β → 0. The Born-level partonic cross section can be written aŝ (2.7) To the first order in β, one can approximate the δ-functions in the above formula as   and the squared-amplitudes can be expanded according to the counting in eq. (2.6). The integrals over the 3-momenta can then be carried out, and we arrive at approximate expressions of the partonic cross sections in the following form where v is the vacuum expectation value of the Higgs field and α s is the strong coupling constant.
A crucial feature of eq. (2.9) is that the leading order partonic cross sections are proportional to β 4 in the threshold limit. This should be contrasted to the case of a similar process, tt production, where the threshold behavior is given by [30], (2.10) That is, the Born partonic cross sections are linear in β in the threshold limit. The different behaviors imply that the threshold region is less important in the case of tth production than for tt production. In order to study these different behaviors more precisely, we will numerically compare the approximate results to the exact ones in the following. For our numerical computations, we set m t = 173.5 GeV, m h = 125.09 GeV and v = 246.22 GeV [31]. The strong coupling constant is evolved from the initial condition α s (m Z ) = 0.1181 to the renormalization scale µ r = m t + m h /2, where m Z is the mass of the Z boson. In Fig. 1, we numerically compare the approximate Born partonic cross sections for tth production with the exact ones, i.e,σ (0)Exact qq(gg)→tth . The approximate results are obtained from Eq. (2.9). To calculateσ (0)Exact qq(gg)→tth , we first employ the programs FeynArts [32] and FeynCalc [33,34] to generate the transition amplitudes and then use the Cuba library [35,36] to perform the phase space integration. Our numerical results are checked against the automatic program MadGraph5 aMC@NLO (MG5) [37]. As shown in Fig. 1, the approximate and exact partonic cross sections approach each other as β becomes small.  They are both highly suppressed in the threshold limit, as can be expected. When β grows larger, the approximate results start to overestimate the exact ones, indicating that there are important negative power corrections to the approximate formula. In Fig. 2, we show a similar comparison for the tt production process. The exact results are obtained from [38], while the approximate results are computed using Eq. (2.10). The first impression is that the small-β region is less suppressed compared to the tth case, which is clear from the β 1 vs. β 4 behaviors. This means that the reliability of small-β resummation can be quite different in these two processes, a fact not often mentioned in the literatures. Beside this, here we also observe that the approximate result in the qq channel overestimates the exact one as β grows. Interestingly, for the gg → tt subprocess, the approximate result underestimates the exact one for most values of β, contrary to the gg → tth case shown in Fig. 1. It is therefore possible that, for the sum of the two channels, the approximate result stays closer to the exact one than in the case of individual channels. This should be investigated at the hadron level as we are going to do below.
The partonic cross sections need to be convoluted with the parton distribution functions (PDFs) to arrive at the hadronic cross sections. It is therefore interesting to see how the approximate and exact results compare at the hadron level. We show in Fig. 3 the hadronic differential cross sections dσ (0) /dβ for tth production. From the left plot, we see again that the approximate results in both the qq and the gg channels overshoot a lot over the exact ones for large β. It is also clear that the small-β region is highly suppressed due to the β 4 behavior. In the right plot, we show the ratio between the approximate results and the exact ones as a function of β. One can see that in the small-β region, the approximate results are in good agreements with the exact ones. As β goes above ∼ 0.3, the approximation quickly fails. On the other hand, we show in Fig. 4 the results for tt production. Again, we find that the small-β region is more important in this case compared to tth production due to the β 1 behavior. Besides, one can also see the accidental cancellation between the power corrections in the qq and gg channel mentioned in the last paragraph. The above two facts lead to the observation that the small-β expansion        provides a reasonable approximation to the hadronic total cross section for tt production [26-28, 38, 40, 41]. This is clearly not the case for tth production. While the above analyses show that the small-β region for tth production is not as important as that for tt production, it is still interesting to study the small-β behavior of the cross section at higher orders in QCD. First of all, theoretically, tth production is similar but slightly different from tt production. The tt pair is recoiled against the Higgs boson and therefore has a non-vanishing transverse momentum already at the lowest order in QCD. The interplay between the soft gluons and the Coulomb gluons can therefore be a bit different from the case of tt production. This poses a question of how to properly factorize these contributions to all orders in perturbation theory, which was not addressed in the literature. Secondly, this process is closely related to the e + e − → tth process at future electron-positron colliders, which receives important QED corrections, especially in the threshold region [42]. The investigation of soft and Coulomb gluons in the pp → tth process can therefore be applied straightforwardly to the soft and Coulomb photons in the e + e − → tth process. Finally, while the small-β limit is not significant for the total cross section, it could be important if one specifically wants to study certain kinematic configurations sensitive to the threshold region by, e.g., vetoing additional jets.
In this work, we are going to study the threshold region by applying a cut on the β variable. While this is not a physical cut (since β cannot be measured), it simplifies the theoretical considerations. We define the following quantity [38] (2.11) Apparently, for sufficiently small β cut , σ cut should be well approximated by the leading power expression of the threshold expansion. It is also obvious that as β cut → 1, σ cut approaches the total cross section defined in Eq. (2.1). The constrained cross section σ cut is the main object we are going to study in the rest of the paper. Before entering the technical details, we first investigate its leading order behavior against the variation of β cut . We again choose to work with the 13 TeV LHC. We take the exact results and the approximate results for the leading order partonic cross sectionsσ At higher orders in perturbation theory, exchanges of Coulomb gluons and soft gluons lead to threshold-enhanced terms such as 1/β n and ln n β. For small β cut , these terms represent the dominant contributions to the constrained cross section σ cut . The rest of the paper will be devoted to deriving a factorization formula for the partonic cross section in the threshold limit, and resumming these enhanced terms to all orders in perturbation theory.

Higher order QCD corrections in the threshold limit
Beyond the Born level, the cross sections receive contributions from exchange of virtual gluons and emission of real gluons. We will investigate these contributions as a power expansion in β in the threshold limit β → 0. In this limit, there will be ln β-enhanced terms and 1/β-enhanced terms at higher orders in α s . Schematically, we are going to consider corrections of the form The collection of these terms are referred to as the improved next-to-leading logarithmic (NLL ′ ) corrections. Note that due to the presence of two kinds of terms, one needs to insist on a consistent logarithmic counting for both of them, which we take as λ ∼ α s ∼ β ∼ 1/ ln β. Using this counting, it can be seen that the NLL ′ corrections include terms up to order λ 1 . It is also clear from this counting that one needs to include formally O(β 1 ) nextto-leading power (NLP) terms besides the O(β 0 ) leading power (LP) ones in the power expansion. This greatly complicates the analysis of factorization, as will be clear below.
The behavior of higher order corrections in the threshold limit can be studied using the method of regions [43,44]. We work in the partonic center-of-mass frame where the momenta of the two incoming partons are given by where n andn are two light-like vectors satisfying n 2 =n 2 = 0 and n ·n = 2. For a given momentum k, we perform the light-cone decomposition as with k + =n · k and k − = n · k. We identify the following momentum regions relevant to our problem: hard : These serve as the basis for constructing the effective field theoretic description of the process, and for deriving the factorization formula for the cross sections. At this point, it should be noted that there is a subtle difference between tth production here and tt production discussed in [26-28, 41, 45]. In tt production, the 3-momentum of the tt pair is of the ultrasoft scale √ŝ β 2 . This means that the tt rest frame is formally equivalent to the partonic center-of-mass frame. On the other hand, in tth production, the tt pair is recoiled by the Higgs boson and has a 3-momentum of the potential scale √ŝ β. Therefore, an ultrasoft mode in the partonic center-of-mass frame will become a potential mode in the tt rest frame. The impact of this difference on the factorization and resummation will be discussed in this section.

Effective field theories
In order to derive the factorization and resummation formulas in the threshold limit, it is useful to employ the language of effective field theories (EFTs). According to the momentum regions in Eq. (3.4), the relevant EFTs are the soft-collinear effective theory (SCET) and the non-relativistic QCD (NRQCD).
SCET [46][47][48][49][50] describes the interactions among collinear, anticollinear and ultrasoft modes. At leading power and next-to-leading power in β, the effective Lagrangians are given by where ξ n and q us denote the collinear and ultrasoft quark fields; A n and A (us) represent the collinear (ultrasoft) gluon fields, with F µν n(us) their field strength tensors; W n is the collinear Wilson line.
To describe the interactions among the potential, soft and ultrasoft modes, we employ the potential non-relativistic QCD (pNRQCD) [51][52][53][54]. The leading power and next-toleading power effective Lagrangians can be written as [53][54][55] where ψ and χ are Pauli spinor fields annihilating the top quark and creating the antitop quark, respectively; E i us = F i0 us are the chromoelectric components of the ultrasoft field strength tensor. The coefficient a 1 was calculated in [56,57] and is given by a 1 = 31C A /9 − 10n f /9. The one-loop coefficient β 0 of the QCD β-function is given in the Appendix. Note that in the pNRQCD power counting, α s and r in Eq. (3.9) and (3.11) are considered as order β. Worthy of particular attention is that in Refs. [53][54][55], the pNRQCD Lagrangian is derived in the rest frame of the quarkonium, where the heavy quark pair is recoiled by ultrasoft momenta. This is different from our case of tth production, where the top quark pair is recoiled by the Higgs boson. However, since the LO and NLO potentials in Eqs. (3.9)-(3.11) only involve the relative momentum between the heavy quark and antiquark, in our case the LP and NLP Lagrangians take the same form. It should be stressed that this fact is not expected to hold beyond NLP. For example, as shown in Appendix B, new structures depending on the recoil momentum appear in the NNLP Lagrangian.
To derive the factorization formula, we first match the QCD amplitudes onto an effective Hamiltonian constructed out of the SCET and pNRQCD fields. Generically, we write where I labels different color structures and m for Lorentz structures, O Im LP and O Im NLP are leading power and next-to-leading power effective operators describing the scattering process under consideration, while C Im LP and C Im NLP are their Wilson coefficients arising from the hard region contributions. Note that the NLP effective Hamiltonian H NLP actually does not contribute to the cross section at next-to-leading power. The reason is that such a contribution would be given by the interferences between O Im LP and O Im NLP , which vanish due to angular momentum conservation.
To obtain the Wilson coefficients C Im LP , we need to calculate on-shell scattering amplitudes in the limit √ŝ = 2m t + m h using both QCD and the effective Hamiltonian. In this limit, the loop integrals in the effective theories are scaleless and vanish in dimensional regularization. Therefore, we only need to calculate the QCD amplitudes up to NLO. For this calculation we employ the program packages FeynArts [32], FeynCalc [33,34] and FIRE5 [58] to generate the amplitudes and perform the reduction to master integrals. The resulting master integrals can be evaluated to analytic expressions using Package-X [59,60]. The interface connecting Package-X, FIRE5 and FeynCalc is provided by FeynHelpers [61]. We renormalize the top quark mass in the on-shell scheme, and the strong coupling α s in the MS scheme. The Wilson coefficients can then be extracted after renormalizing the effective operators (which is equivalent to subtracting the infrared poles from the QCD amplitudes). For the purpose of this paper, we don't need the explicit forms of individual Wilson coefficients and effective operators, but the combinations of them entering the cross section for tth production. These combinations will be given in the next subsection as the "hard functions".

Factorization at leading and next-to-leading power
The discussion in the last subsection tells us that for the cross section up to NLP, we only need to consider the amplitudes of H LP . At leading power, we use H LP together with the LP Lagrangians L 0 SCET and L 0 pNRQCD to calculate the cross section where ij = qq, gg labels the initial state, and dΦ f represents phase-space integration over the momentum p f of the final state particle f . In the LP Lagrangians (3.5) and (3.9), the interactions of the ultrasoft gluon with the collinear fields and heavy quark fields are encoded in the covariant derivative D us . Such interactions can be removed by the decoupling transformations [27,48] where S v (x) and S q n(n) are ultrasoft Wilson lines in the fundamental representation along the directions implied by the subscripts, while S g n(n) are ultrasoft Wilson lines in the adjoint representation.
After the decoupling, the factorization of the LP cross section follows along the same line of arguments as the tt production [27]. The only difference comes from the appearance of the Higgs momentum p h . The factorization formula therefore readŝ where The potential function is given by where {a} = {a 1 , a 2 , a 3 , a 4 } and P α {a} are the projectors to the singlet-octet color states of the tt system, which are given by Note that since the interactions in the LP Lagrangian L 0 pNRQCD are spin-independent, the two ψ fields in the definition of the potential function share the same polarization index s 1 , and similar for the two χ fields. It should be stressed that the potential function here is different from that in tt production, due to the presence of the recoil momentum p h .
The soft function in Eq. (3.15) is defined as where the color basis for the qq channel is given by 20) and that for the gg channel is gg,{a} = 1 2 This soft function is the same as that for tt production in [27], since it does not feel the presence of the recoil momentum. It is diagonal in the color basis we have chosen. In practice, it is more convenient to quote its Laplace transform where Up to the NLO, the result reads [27] s αII ij (L, µ) = 1 + where C i = C F for i = q,q, C i = C A for i = g, C α = 0 for α = (1) (singlet), and C α = C A for α = (8) (octet). Note that the soft function actually does not depend on the index I. Finally, the hard function H IJ ij is defined as the product of LP Wilson coefficients C Im * LP C Jm LP projected onto the S-wave spin structure determined by the potential function. In general, the hard function is not diagonal. However, since the soft function is diagonal, only the diagonal entries of the hard function contribute to the cross section. Furthermore, since the soft functions αII ij does not depend on the index I, we can project the diagonal entries of the hard function onto the singlet-octet basis as We have calculated these entries H α ij at NLO explicitly using the method described in the last subsection, and write the result as (3.27) The one-loop hard anomalous dimensions γ H,α ij,0 are given by To achieve NLL ′ accuracy for the resummation, we need to further consider nextto-leading power corrections to the cross section. These amount to contributions from the NLP interactions in L 1a,1b,1c SCET and L 1a,1b pNRQCD to the squared-amplitude in Eq. (3.13). In analogy to the arguments in [26,27] for tt production, it can be shown that single insertions of L 1a,1b SCET give vanishing results due to angular momentum conservation, while L 1c SCET does not contribute due to baryon number conservation. The terms in L 1b pNRQCD involve subleading potentials between the top and anti-top quarks. These contributions can be incorporated by upgrading the potential function J α (q) to the NLO, which we will discuss in the next subsection. Finally, we need to consider the corrections induced by L 1a pNRQCD . The terms in L 1a pNRQCD involve extra interactions between ultrasoft gluons and potential modes, which are not removed by the decoupling transform (3.14). In [27], it was proved that these interactions do not contribute to the tt cross section at NLP. However, the arguments there rely on the fact that the partonic center-of-mass frame and the tt rest frame are the same, and therefore the potential function does not depend on an external 3-momentum. However, for tth production, the recoil momentum p h from the extra Higgs boson spoils the proof, and we need to reinvestigate the contributions from L 1a pNRQCD here. We begin with an explicit diagram depicted in Fig. 6. Its contribution to the partonic cross section in the threshold limit can be written as (up to overall factors due to coupling constants, color factors, Wilson coefficients, etc.) where A( p h ) is given by where E J is defined in Eq. (3.16), and Note that we have suppressed the imaginary part +iε in the propagators. We now observe that the last factor in the above expression does not depend on k, p t , pt and p h , while the other factors do not depend on n. Together with the fact that v = 0, we can conclude that after integrating over k, p t , pt and p g , the function A( p h ) must be proportional to n · p h (multiplied by a function of | p h | 2 and other scalar quantities). As a result, after performing the integration over p h as in Eq. (3.30), the contribution of this diagram to the partonic total cross section must vanish. The argument above can be generalized to all contributions from a single insertion of L 1a pNRQCD in a more formal way. The cross section induced by L 1a pNRQCD can be written aŝ where T denotes time-ordered product. We can perform the usual decoupling transforms (3.14) to remove the leading power interaction between ultrasoft and potential modes. The remaining interaction is of the x · E us form from L 1a pNRQCD . As a result, we can write the cross section aŝ where we have suppressed all color indices for simplicity, while E J and p J are given in Eq. (3.16). The subleading potential function and soft function are defined as where again we have ignored all color structures which are not important for the arguments here. Note that j 1a (E q , q, k 0 ) must be proportional to q since this is the only 3-vector it can depend on. For tt production, there is no recoil momentum and p J = 0. Therefore in the integrand for the cross section one has j 1a (E, 0, k 0 ) = 0, and one can conclude that the contribution from L 1a pNRQCD to the cross section vanishes. This is essentially the argument in [27]. For tth production, due to the presence of a recoil momentum, j 1a can depend on p J = − p h , and is not zero in general. However, note that the whole integrand in Eq. (3.35) is an odd function of p h . Consequently, after integrating over the phase space dΦ h of the Higgs boson, the contribution still vanishes. Therefore, we arrive at the same conclusion as in the tt case that the only NLP contribution to the total cross section comes from L 1b pNRQCD . We emphasize that this fact only holds at the level of total cross section, and extra corrections may be present if one does not integrate over the momentum of the Higgs boson.
In summary, up to the next-to-leading power, the cross section can be factorized aŝ where the hard function H and soft function S only receive leading power contributions. The potential function J α (q) contains both LP and NLP contributions, which we present in the next subsection. Note that this simple form of the factorization formula is not expected to hold at higher powers in β, as we'll discuss in Appendix B.

The potential function with a recoil momentum
As introduced in the last subsection, a non-trivial difference between tth production and tt production is the dependence of the potential function J α (E q , q) on the recoil momentum p h of the extra Higgs boson. Up to the NLO, we write the potential function as where the LP term is defined in Eq. (3.17), and the NLP term is given by

(3.40)
In calculating the potential function, one needs to consider the interactions induced by the leading power Lagrangian L 0 pNRQCD to all orders in α s . In this way, one resums all terms of the form (α s /β) n and α s (α s /β) n . According to [27], the potential function can be related to the imaginary part of the pNRQCD Green function G α ( r 1 , r 2 ; E q , q) of the tt pair at origin: where the Green function is defined as where we have suppressed the color and Lorentz indices for simplicity.
We would now like to relate the potential function with a recoil momentum J α (E q , q) to the zero-recoil one J α (E, 0) given above. From the perturbative point-of-view, we can write the potential function as a sum of diagrams with arbitrary numbers of insertion of L 0 pNRQCD and up to one insertion of L 1b pNRQCD : where the ellipsis denotes more insertions of heavy quark propagators and LO potential terms. Note that the perturbative expansion of J α (E, 0) is the same as the above if we identify E = E q − | q| 2 /(4m t ). This identity can also be seen if we evaluate the potential function (3.17) and (3.40) in the tt rest frame. We can therefore deduce the relation from which we obtain the potential function we need.

Resummation
After deriving the factorization formula (3.38), we will now perform the resummation of both 1/β and ln β enhanced corrections at the NLL ′ accuracy. Schematically, the resummed result takes the form where the counting rule is in accordance with [41,45]. The 1/β terms in the above expression are contained in the potential function J α (E q , q) in the factorization formula. To resum the ln β terms, we need to evaluate the hard function H α ij (µ) and the soft function S α ij (ω, µ) separately at the hard scale µ h and the soft scale µ s , and then use their renormalization group equations (RGEs) to evolve them to the factorization scale µ f .
The Laplace-space soft functions α ij satisfies the RGE where the cusp anomalous dimension γ cusp (α s ) is given by [63] γ cusp (α s ) = α s π + α s 4π 50) and the soft anomalous dimension γ s,α (α s ) is [45] γ s,α (α s ) = − α s π C α + O(α 2 s ) . (3.51) Here we recall that C α = 0 for α = (1) (singlet) and C α = C A for α = (8) (octet). The RGE for the hard function is given by where the hard anomalous dimension can be expanded as with the one-loop coefficients given in Eq. (3.28). The method to solve the RGEs (3.49) and (3.52) is standard [64,65]. Plugging the result back to the factorization formula (3.38), we obtain the resummed cross section aŝ The explicit expressions for the QCD β-function and the evolution factor U α ij (µ f , µ h , µ s ) are collected in Appendix A.
An important validation of the resummed formula (3.54) is to compare its fixed-order expansion to explicit calculations. We define the coefficients of the expansion up to the NLO according toσ with L S = ln(β 2 (m h + 2m t )/µ). The above results contain the leading terms at LO and NLO in the threshold limit, and can be compared to explicit calculations performed in [16].
We have found complete agreement between the two results.

Matching to the NLO
In this section, we present some numeric results based on our resummation formula Eq. (3.54).
Since the resummation formula contains only up to NLP terms, it is necessary to include higher power contributions whenever possible. For this reason we need to match the resummation formula to fixed-order calculations in order to extend the validity of our predictions beyond the threshold limit. To take into account the higher power effects contained in the exact LO and NLO results, we use the following matching formulâ This will be the main quantity for which we present numeric results later.

The choice of the scales
The resummed formula (3.54) involves a set of unphysical scales, which are the hard scale µ h , the soft scale µ s and the factorization scale µ f . In additional, although the potential function J α (q) is formally independent of a renormalization scale, its finite order truncation has a residue scale-dependence. We denote this scale as µ J . These scales must be chosen appropriately to improve the convergence of the perturbation theory. This subsection is devoted to this issue. From the explicit form of the hard function, one observes the appearance of the logarithms ln n ((2m t + m h )/µ) at higher orders in perturbation theory. It is therefore natural to choose µ h ∼ 2m t + m h to make these logarithmic corrections under control. For the potential function, one observes logarithms of the form ln(4m t β 2 /µ 2 ). However, the natural choice µ J ∼ 2m t β requires an infrared cut-off, since the variable β is integrated over from 0. We therefore choose scale µ f should be chosen to be appropriate in fixed-order calculations, since it enters the matching formula (4.1). For that we take the conventional choice µ default f = m t + m h /2. The choice of the soft scale µ s is more subtle. In the literature there are two kinds of methods for that purpose. One is to choose µ s in the Laplace moment space, e.g., µ s ∼ (2m t + m h )/N where N is the moment variable entering the Laplace transform Eq. (3.22). Another is to choose µ s in the momentum space, where µ s ∼ (2m t + m h )β 2 . We will use the latter method, for which we need to impose an infrared cut-off in order to avoid the Landau pole when integrating over β. Namely we have The value of µ cut s should be much larger than Λ QCD , but should not be too large since that will reintroduce large logarithms into the soft function. To study the effect of varying µ cut s , we define the following quantity represents the contributions from the soft logarithms to the cross section in the threshold limit. We further define the ratio κ soft as where the denominator is simply the leading order part of the numerator. The above quantity represents the relative size of the soft corrections, and can be used to motivate an appropriate choice for µ cut s . In Fig. 7, we show the numeric values of κ soft as a function of µ cut s . It can be seen that the size of the soft corrections stays stable when µ cut s < 20 GeV, and increases dramatically when going beyond. Therefore, we choose the default value of µ cut s to be 20 GeV in our numeric study.

Results and discussions
In this subsection, we present the numeric results for the total cross section at 13 TeV and 14 TeV LHC. For readers' convenience, we list here again the parameters we use: m t = 173.5 GeV, m h = 125.09 GeV and v = 246.22 GeV. We have employed the MMHT2014 (N)LO PDFs [39] with the corresponding α s (m Z ). We begin with the scale dependence of the total cross section at 13 TeV LHC. The result at 14 TeV LHC is similar and we do not show it here. The LO and NLO cross sections depend on the factorization scale µ f , where the strong coupling α s and the PDFs are evaluated. The matched NLL ′ +NLO cross section depends in addition the hard scale µ h , the soft scale µ s and the potential scale µ J . In the left plot of Fig. 8, we show the dependence of the NLL ′ +NLO cross section on µ h , µ s and µ J . We observe that the dependence is rather mild. This can actually be expected since these scales only affects the region β < β cut , which does not make dominant contributions. In the right plot of Fig. 8, we show the dependence of the LO, NLO and NLL ′ +NLO cross sections on the factorization scale µ f . It can be seen that the µ f dependence is significantly reduced when going to higher orders in perturbation theory. At NLL ′ +NLO, the residue µ f dependence is merely about 2%. To estimate the theoretical uncertainties of the NLL ′ +NLO predictions, we vary the 4 scales up and down by a factor of 2, and add the resulting variations of the cross sections in quadrature.
The predictions for the total cross sections are summarized in Table 1. The K-factor is defined as the ratio of the NLL ′ +NLO cross section to the NLO one. Although we do not expect huge effects from the resummation due to the β 4 suppression of the threshold region, we still find that the higher order threshold corrections enhance the cross section by 6%. This should be taken into account for the high precision HL-LHC run. We also observe a significantly reduction of the scale dependence from NLO to NLL ′ +NLO. This should be contrasted to the result of [22], where the NLL resummation was performed for the ln β corrections, while the 1/β corrections were added at fixed-order (i.e., not resummed).
The difference between our result and the result in [22] resides in the fact that we have computed the NLO hard function exactly for each color channel, and we have resummed the 1/β terms to all orders into the potential function. The reduction of the scale dependence shows that these additional efforts are important phenomenologically.

Conclusion
In this work, we have generalized the resummation framework for tt production to investigate the associated production of a Higgs boson with a pair of top quarks. A major difference between the two processes is that for tth production, the tt pair is recoiled by the Higgs boson and has a residue momentum of the potential scaling; while for tt production, the residue momentum of the top quark pair is of the ultrasoft scaling. The presence of the recoil momentum leads to several complications in the derivation of the factorization formula. We have shown that the next-to-leading power interaction between the ultrasoft mode and the potential mode does not contribute to the total cross section when the momentum of the Higgs boson is integrated over. We have also argued that the contributions from the potential mode can be resummed into a potential function, which is related to that in tt production via a boost. The final outcome of these considerations is that the total cross section for tth production admits a similar factorization formula up to NLP as that for tt production. This similarity relies on subtle cancellations of the ultrasoft-potential interactions in the integrated cross section and on the same form of pNRQCD Lagrangians up to NLP, which may not hold at higher powers in β. An important ingredient entering the factorization formula is the hard function, which was not known in the literature beyond the LO. We have explicitly calculated the NLO corrections to the hard function, decomposed into singlet and octet color configurations. We have validated our factorization formula by expanding it to the NLO in α s and comparing with the explicit calculations in [16]. Based the factorization formula, we have derived a resummation formula at NLL ′ accuracy using RG equations. By matching to the NLO result, we are able to provide numeric predictions for the total cross sections at the NLL ′ +NLO accuracy. We find that the resummation effects enhance the cross sections at 13 TeV and 14 TeV LHC by about 6%, and reduce the scale dependence significantly.
We emphasize that the resummation framework in this paper can as well be applied to tth production at a future e + e − collider, where β is fixed by the collider energy instead of being integrated over. The impact of resummation is expected to be more important in that case if the collider energy is not too far beyond the production threshold.

A Ingredients in the renormalization group evolution
In this appendix, we list the ingredients entering the resummation formula (3.54). We begin with the QCD β-function, whose perturbative expansion is defined as where the one-loop and two-loop coefficients are given by [66] The evolution factor in the resummation formula is defined as B New structures at next-to-next-to-leading power In this work, we have only been concerned with up to next-to-leading power contributions. We have seen that at this order, the total cross sections for tth production and tt production have a similar form of factorization in the threshold limit. However, this similarity is in general not expected to hold at higher powers in β. In this appendix, we discuss some next-to-next-to-leading power (NNLP) contributions which may lead to new structures in the factorization formula. Firstly, the vanishing of the contribution from the x · E us term strongly relies on the fact that at O(β), the integrand for the cross section must be linear in p h . This will no longer be true at NNLP, where one can have corrections quadratic in p h , which do not vanish even after integrating over all phase space. This may lead to new terms in the factorization formula which are not present in tt production.
Secondly, at NNLP the third-order terms in the effective Lagrangian will come into play. For example, the NNLP pNRQCD Lagrangian contains where ψ ≡ ψ(x 0 , x + r) and χ ≡ χ(x 0 , x). Note that the term V nr will also appear in tt production at NNLP, but the term V rc is one power higher in tt production, since it involves the 3-momentum of the tt pair. The different counting of the V rc term may spoil the simple relation (3.47) at NNLP. Finally, note that the above discussions are just speculations, and we cannot draw a definite conclusion unless we perform a more thorough analysis, which is beyond the scope of the current work.