Two-loop infrared singularities in the production of a Higgs boson associated with a top-quark pair

The associated production of a Higgs boson and a top-quark pair is important for probing the Yukawa coupling of the top quark, and calls for better theoretical modeling. In this paper, we calculate the two-loop infrared divergences in $t\bar{t}H$ production at hadron colliders. To do that we compute the one-loop amplitudes to higher orders in the dimensional regulator $\epsilon$. Numeric results for the infrared poles are given as a reference at several representative phase-space points. The result in this work servers as a part of the ongoing efforts towards the $t\bar{t}H$ cross sections at the next-to-next-to-leading order.


Introduction
The associated production of a top quark pair and a Higgs boson is one of the most important processes to study the Yukawa coupling of the top quark at the Large Hadron Collider (LHC) and the next-generation experimental facilities. The Yukawa coupling is crucial to understanding the origin of the large mass of the top quark. It can also probe the possible violation of the CP symmetry in the top quark sector [1,2]. Such a violation is required to generate the matter-anti-matter asymmetry in our observable universe. With an integrated luminosity up to 139 fb −1 , the LHC Run 2 has measured the cross section for this process to a relative accuracy of about 20% [1][2][3][4]. With the accumulation of much more data in the near future, the experimental precision is expected to be significantly improved. In order to extract the top Yukawa coupling from the high precision cross section measurements, it is necessary to have high precision theoretical predictions for the relevant observables.
In quantum chromodynamics (QCD), the total and differential cross sections for ttH production at the next-to-leading order (NLO) have been known since almost twenty years ago [5][6][7][8][9][10]. Approximate next-to-next-to-leading order (NNLO) as well as soft gluon resummed results are also calculated in [11][12][13][14][15][16][17][18]. These approximations are only valid in certain kinematic limits, and it is highly desirable to have a complete NNLO calculation for this process, in order to control the theoretical uncertainties. However, such a calculation is still out-of-reach due to the obvious obstacles from the complicated two-loop amplitudes. A part of the NNLO contributions that do not involve two-loop amplitudes is recently available [19].
A prominent property of gauge theory amplitudes is the existence of infrared (IR) singularities. Understanding their structure is crucial for designing IR safe quantities that can be compared to experimental measurements. Even in IR safe observables, in certain kinematic regions, there can be large logarithms which need to be resummed to all orders in perturbation theory. The structure of these logarithms is, again, governed by the IR behaviors of scattering amplitudes. In the past couple of years, significant progress has been achieved in the understanding of the IR singularities in non-abelian gauge theories, both in massless  and massive [44][45][46][47][48][49][50][51][52][53][54][55][56][57][58] cases. Nevertheless, applying these universal behaviors to a given scattering process still requires a considerable amount of work.
In this paper, we apply the method in [52,53] to calculate the IR poles in the two-loop amplitudes for the ttH production process. The biggest challenge in this calculation is that we need to compute the one-loop amplitude to higher orders in the dimensional regulator = (4 − d)/2, where d is the dimension of spacetime. Unlike the tt case, the one-loop ttH amplitudes involve many 4-point and 5-point integrals whose higher-order coefficients in are not known in the literature. Hence, a major part of this paper is devoted to the systematic calculation of these integrals. These integrals at higher orders in are also required for the finite part of the NNLO cross sections. This paper is organized as follows. In Section 2 we introduce our notations and review the generic structure of IR singularities of two-loop scattering amplitudes in non-abelian gauge theories. In Section 3 we show the details of the calculation of the one-loop am-plitudes. In particular, we demonstrate how to construct canonical differential equations for the master integrals. In Section 4, we give our results for the IR poles at several representative phase-space points, and briefly summarize our work. We leave some lengthy expressions in the Appendix.

Notations and structure of IR singularities
For the production of a Higgs boson associated with a top-quark pair, we consider the partonic processes where α, β, k, l, a, b are color indices. We define the following kinematic variables: where σ i = +1 if p i is incoming, and σ i = −1 if p i is outgoing. We use the color space formalism [59,60] where the amplitudes are vectors |M q,g . The subscript q or g specifies the quark-antiquark annihilation channel or the gluon fusion channel, respectively. For the qq channel, we choose the independent color structures as For the gg channel, we use the color basis The UV divergences in the amplitudes are renormalized according to where Z g , Z q and Z Q are the on-shell wave-function renormalization constants for gluons, light-and heavy-quarks, respectively. We have suppressed the dependence of the amplitudes on other kinematic variables. The Yukawa coupling g Y is defined as .
We renormalize the top quark mass in the on-shell scheme: m 0 t = Z m m t , and the Yukawa coupling is renormalized accordingly. The strong coupling constant α s is renormalized in the MS scheme with n f = n l + n h active flavors. The relations between the bare couplings and the renormalized ones are given by The renormalization constants are After UV renormalization, the remaining IR divergences can be subtracted by a multiplicative factor Z −1 ( , m t , µ), where the bold symbol denotes an operator in color space. More precisely, we have The Z factor satisfies a renormalization group equation (RGE) of the form where Γ is a universal anomalous-dimension operator, which has been calculated up to order α 2 s in [51][52][53]. In the ttH production process, the anomalous-dimension matrices for the qq and gg amplitudes are rather similar as those in the tt production process considered in [53,61]. They are given by [12] Γ qq = C F γ cusp (α s ) log −s 12 µ 2 + C F γ cusp (β 34 , α s ) + 2γ q (α s ) + 2γ Q (α s ) 1 (2.12) The cusp angle β 34 is defined by (2.13) The perturbative expansions of γ cusp , γ q , γ g and γ Q can be found, for instance, in the Appendix of [61]. The only difference of the anomalous-dimension matrices here (with respect to those in tt production) is thats 13 =s 24 ands 14 =s 23 due to the 2 → 3 kinematics.
Both the UV-renormalized amplitudes and the IR subtraction factors can be expanded in powers of α s : (2.14) We may then extract the IR singularities of the amplitudes order-by-order in α s : Note that to predict the IR poles at the two-loop order, one must calculate the UVrenormalized one-loop amplitudes to O( 1 ). They multiply the 1/ 2 terms in Z (1) q,g and give rise to 1/ divergences. The next section is devoted to this non-trivial task.
3 Calculation of the one-loop amplitudes to higher orders in

Setup
As is clear from the last section, in order to predict the IR poles at two loops, we need the one-loop amplitudes up to order 1 . We generate the amplitudes using FeynArts [62], and manipulate the expressions with FeynCalc [63]. We then need to express the amplitudes in terms of scalar integrals, there are two ways to achieve this. Since we are interested in the IR poles in the interference of M (2) with M (0) , we may readily multiply the oneloop amplitudes by the tree-level ones. The Lorentz contractions and Dirac traces can now be easily performed while keeping the color information. Alternatively, we may also apply a complete set of projectors to the one-loop amplitudes, and extract the coefficients as a linear combination of scalar integrals. 1 The second method is more complicated for our purpose, but the results can be useful if one wants to obtain the one-loop amplitude squared. We have performed the calculation in both ways, and the results agree. The one-loop scalar integrals can be categorized into 12 families (topologies), four of which are independent (and the others can be obtained with exchanges of external momenta). Each family is defined by 5 propagator denominators denoted as D i (i = 1, . . . , 5): where a ≡ a 1 + a 2 + a 3 + a 4 + a 5 . The prefactor is introduced such that the integrals are dimensionless and do not contain γ E in their series expansions in . The propagator denominators for the four independent families are chosen as The corresponding diagrams are depicted in Figure 1. The remaining 8 topologies can be obtained by the exchanges p 1 ↔ p 2 and/or p 3 ↔ p 4 . There are 18, 20, 22 and 22 master integrals in the topologies A, B, C and D, respectively. To calculate the master integrals, we adopt the method of canonical differential equations [64]. Namely, we construct linear combinations of the master integrals which satisfy a set of differential equations of the -form. We denote such a "canonical basis" as f . Integrals in such a basis have the property of uniform transcendentality (UT), and hence are also dubbed "UT integrals". We introduce the dimensionless kinematic variables The differential equations can then be written as where x denotes a set of independent kinematic variables chosen from x ij and x h (note that the choices of independent variables are different for each topology). The matrix dA takes the d log-form: where C i are matrices consisting of rational numbers, and W i ( x) are algebraic functions of the kinematic variables. The functions W i are called the "letters" for this topology, and the set of all independent letters is called the "alphabet". The canonical differential equations (3.4) can be solved order-by-order in . To this end, we expand the (suitably normalized) integrals as Taylor series where the nth order coefficient function can be written as a Chen iterated integral [65] In certain cases, these iterated integrals can be solved analytically (either by direct integration or by bootstrapping). The results can often be written in terms of generalized polylogarithms (GPLs) which allow efficient numeric evaluation [66][67][68]. When an analytic solution is not available, it is straightforward to evaluate them numerically either by numerical integration or by a series expansion [69,70]. In the rest of this Section, we discuss the construction of the canonical basis f and the matrix dA( x).

The canonical master integrals
We use the method of [71] to construct the canonical bases using the Baikov representation [72]. We present the results with generic external momenta and internal masses. The results for the ttH process can be obtained by inserting the momenta and masses associated with the propagator denominators for each topology. Consider a generic one-loop integral family with N = E + 1 external legs, where E is the number of independent external momenta. Integrals in this family can be written as where z i are the propagator denominators given by The corresponding Baikov representation is given by where z = {z 1 , . . . , z N } is the collection of the Baikov variables (i.e., propagator denominators). The function G N (z) is a polynomial of the N variables, while K N is independent of z. They are given by where the Gram determinant is defined as The UT integrals g N for any N is obtained in [71]. For the purpose of this work, we need them up to N = 5. They are given by Note that g 3 and g 4 can be straightforwardly identified as Feynman integrals in 4 − 2 dimensions. On the other hand, g 1 and g 2 can be naturally regarded as creatures in 2 − 2 dimensions, while g 5 lives in 6 − 2 dimensions. Namely we have: (3.14) The 2 − 2 and 6 − 2 dimensional integrals can be expressed as Feynman integrals in 4 − 2 dimensions using the dimensional recurrence relations [73,74]. Applying the above to all sectors of a family, we build a complete canonical basis satisfying -form differential equations. As a final remark, we note that there is a freedom in multiplying a UT integral by some complex number, and the result remains UT. We will use this freedom when writing down the canonical basis for ttH production listed in the Appendix.

Bootstrapping the coefficient matrices in the differential equations
Given the UT integrals in (3.14), it is straightforward to calculate their derivatives with respect to some kinematic variable x i : where the elements in the matrix A i ( x) have the property that they only contain simple poles. We would like to combine these derivatives into a total derivative as in Eq. (3.4). We will achieve this by bootstrapping. According to Eq. (3.5), we can write the total derivative as Since A i ( x) are known, it is easy to extract the coefficient matrices C j once we know the letters W j ( x). We obtain the full alphabet using the method described in [75][76][77][78]. We write the differential equation satisfied by an N -point UT integral g N (see Eqs. (3.13) and (3.14)) as 18) where g N ( x, ) and g In the following, we present the generic form of the letters W N and W (i) N,m obtained in [78]. For each m we take g m to be the UT integral with the denominators z 1 , . . . , z m , and show the corresponding W N,m . The other ones can be obtained by rearranging the order of denominators.
The self-dependent letter W N is given by 20) where G N ≡ G N (0). As an example, we consider a 5-point integral in the topology A of ttH production. The Gram determinants entering the letter W 5 are given by: Here we have omitted the dimension of Gram determinants since they only appear in dimensionless letters. The letter W N,N −1 with odd N is given by where the extended Gram determinant is defined as G(q 1 , . . . , q n ; k 1 , . . . , k n ) = det (3.24) -10 -We again show an example in topology A for ttH production. Consider the contribution from F 1,1,1,1,0 to the differential equation of F 1,1,1,1,1 , we need the Gram determinants and K 5 has been shown previously. Plugging these into Eq. (3.22), we readily obtain the letter W 5,4 . The letter W N,N −1 with even N is given by (3.26) We give an example in topology D. Consider the contribution from the 3-point integral F 1,0,1,1,0 to the derivative of the 4-point integral F 1,1,1,1,0 . The Gram determinants we need are: We can then easily obtain the letter W 4,3 by plugging the above into Eq. (3.26). The letter W N,N −2 with odd N is given by where We give an example in topology A for ttH production. Consider F 1,1,0,1,0 's contribution to F 1,1,1,1,1 , we need

30)
and K 5 is already given before. For the letter W N,N −2 with even N , there are two possible cases. In the first case, both the two (N − 1)-point integrals between g N and g N −2 are masters, and the letter is given by where D N = D N (0) and As an example, we show the letter in the contribution from F 1,0,0,0,1 to F 1,0,1,1,1 in topology D. The Gram determinants are The second possibility is that one (or both) of the two (N − 1)-point integrals is not a master and can be reduced to lower-point integrals. Here we only show the letter for the case where the integral g (1) N −1 with the denominators z 1 , . . . , z N −1 is reducible, while the other (N − 1)-point integral g (2) N −1 is a master. In this case, the letter is given by where we use the superscript (i) to label the Gram determinants associated to the integral g (i) N −1 . As an example, we consider the contribution from F 0,1,0,1,0 to F 0,1,1,1,1 in topology D. The Gram determinants are given by It happens that G 2 = K (1) 3 in this case, and the letter is simply given by W 4,2 = K 4 / G (2) 3 . We finally consider the letter W 5,2 . Such dependence can only be present if at least one of the 3-point integrals g (i) 3 is reducible. In this case W 5,2 is the same as (or a combination of) W 5,3 in Eq. (3.28).
The full alphabets for all four topologies in ttH production are collected in electronic files attached with this paper. With the alphabets at hand, we reconstruct the dA matrices which are also attached. These completely fix the differential equations for the master integrals.

Boundary conditions and solution to the differential equations
Given the canonical bases and their differential equations, we still need to know the value of the integrals at some boundary point x 0 . The boundary points should be chosen for each family such that the integrals become simpler. Our strategy is to utilize the spurious singularities in the differential equations. At these points, many UT master integrals vanish; or equivalently, there are further relations among the master Feynman integrals. It turns out that the most difficult boundary conditions are the 5-point integral in topology A and B, and we discuss their determination in the following.
For topology A we choose the boundary point x 0 to be Note that all 4-point integrals can be reduced at the boundary. The lower-point integrals can be easily evaluated as a function of x 12 using Feynman parameters. The 5-point integral is more difficult, and we obtain its boundary value using a small trick. First of all, the 5-point integral only appears in one UT master integral (f 18 in Eq. (A.1)). This corresponds to g 5 in Eq. (3.14), which is UV/IR finite. Hence we have To determine the coefficient c 3 , we exploit the differential equation with respect to x 24 , while keeping the other variables (collectively denoted as x ) at the boundary: where the list A(x 24 ) is given by The differential equation is singular at x 24 = 1/2. On the other hand, this is not a genuine singularity of the integrals. Hence, we know that the expansion coefficients f It is then required that lim which allows us to determine f  18 ( x 0 , 0). The results are expressed as GPLs of the argument 1/2 (with indices 0, 1/2 and ±1), which can be converted to polylogarithms according to [79]. The result for c 3 is given by We now turn to topology B. There are 20 master integrals in this family, and only f 20 involves the 5-point integral F 1,1,1,1,1 . We choose the boundary point x 0 to be The lower-point integrals at the boundary can be easily calculated. We again need to determine the order 3 coefficient of b 20 ( ) ≡ f 20 ( x 0 , ). We now exploit the differential equation with respect to x h (while keeping the other variables at the boundary), and use the fact that F 1,1,1,1,1 is reducible at the point x h = −1. The boundary condition b 20 ( ) can then be obtained by evolving from x h = −1 to x h = 0 using the differential equation. However, the resulting expression contains GPLs with complicated indices and arguments such as On the other hand, the only letter in the alphabet that survives at the boundary is d log(2), and it is clear that the boundary conditions can only contain powers of log(2) in additional to transcendental constants such as zeta values. At weight 3, the only possibilities are π 2 log(2), ζ 3 , log 3 (2) (there can be no imaginary part at the boundary we choose). We use GiNaC to compute the GPLs at the boundary to very high precision, and then use the PSLQ algorithm [80] implemented in PolyLogTools [81] to fit the rational coefficients. Finally we find in topology B that With the boundary conditions at hand, we are now ready to solve the differential equations to get the values of the master integrals at any phase-space point. This step can be done analytically for not-so-complicated integrals. In general, the analytic form is not always easy to obtain, especially for integrals involving many square roots in their alphabets. In these cases, we employ the program package DiffExp [70], which can take a set of -form differential equations and compute the solutions numerically using series expansions along a path in the phase space. For a balance between computation speed and accuracy, we aim at results with a relative precision of 10 −10 . We cross-check the numeric results for the master integrals with analytic expressions whenever the latter are available. In phase-space regions where the sector decomposition method can get reasonably accurate results, we also cross-check our results to those of FIESTA [82] and pySecDec [83,84], and find agreement within the accuracy. Assembling the integrals into the one-loop amplitudes, we obtain the amplitudes up to O( 1 ). These serve as an important building block in the two-loop IR divergences we are going to show in the next Section. We have compared the amplitudes up to O( 0 ) against the results of the program packages GoSam [85] and OpenLoops [86][87][88][89], and find complete agreements. 2

Numeric results and summary
We now come to the main results of this paper, namely the predictions for the IR poles in the two-loop amplitudes for ttH production. In practice, it is more convenient to show the interference between the two-loop amplitudes with the tree-level ones, which is of phenomenological interest. We decompose the color-and spin-summed interference terms into several color coefficients according to [90]: In Tables 1, 2, 3 and 4 we list the numeric values of the IR poles as color coefficients at four representative phase-space points. The first point corresponds to the bulk region M ttH ∼ 550 GeV where the differential cross sections are large. The second and third points are in the high energy region where M ttH ∼ 2 TeV. At the second point, the Higgs boson and the top/anti-top quarks are all moderated boosted. At the third point, the top/anti-top quarks are highly boosted while the Higgs boson is produced at relatively low transverse momentum. Finally, the fourth point is near the production threshold where all final state particles have small energies. These results provide a strong check for a future calculation of the two-loop amplitudes. Even before a full two-loop calculation is available, it is possible to study the amplitudes in various kinematic limits, such as the boosted limit or the threshold limit. Our results in these kinematic regions are therefore useful to validate such calculations.
Comparing to the case of tt production, we observe that there are slight differences in ttH production due to the 2 → 3 kinematics. In particular, the 1/ 2 coefficient of F g h and the 1/ coefficients of G g h and I g lh vanish for µ = m t in tt production. However, they are all non-zero in the case of ttH production.
In summary, we have calculated the two-loop infrared divergences in ttH production at hadron colliders. To do that we have employed the universal anomalous dimensions obtained in [52,53]. We compute the one-loop amplitudes in dimensional regularization up to order 1 , which are important building blocks in the two-loop IR structure. We show the numeric results for the two-loop IR poles at several representative phase-space points. These serve as references for future calculations at this order.
The result in this work is an important part of the ongoing efforts towards the ttH cross sections at NNLO. The one-loop amplitudes calculated in this work can be easily extended to order 2 , which are essential ingredients in the NNLO cross sections. It is interesting to study in more detail the behavior of the IR divergences in the high-energy boosted limit and the low-energy threshold limit, where the amplitudes admit further factorization properties. We leave these to future investigations.    Table 4. IR poles decomposed as color coefficients for the phase-space point x 12 = 3249/400,