{'fraction_non_alphanumeric': 0.039814946619217084, 'fraction_numerical': 0.013551601423487545, 'mean_word_length': 4.316809203057595, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
NLO QCD corrections to t anti-t H production in hadron collisions

The Higgs boson H of the Standard Model can be searched for in the channels p pbar / pp ->t tbar H + X at the Tevatron and the LHC. The cross sections for these processes and the final-state distributions of the Higgs boson and top quarks are presented at next-to-leading order QCD. To calculate these QCD corrections, a special calculational technique for pentagon diagrams has been developed and the dipole subtraction formalism has been adopted for massive particles. The impact of the corrections on the total cross sections is characterized by K factors, the ratios of the cross sections in next-to-leading order over leading order QCD. At the central scale mu_0 = (2 m_t + M_H)/2 the K factors are found to be slightly below unity for the Tevatron (K \sim 0.8) and slightly above unity for the LHC (K \sim 1.2). Including the corrections significantly stabilizes the theoretical predictions for total cross sections and for the distributions in rapidity and transverse momentum of the Higgs boson and top quarks.


Introduction
Even though the Standard Model (SM) has been tremendously successful in describing matter and forces in particle physics, the third component of the model, the Higgs mechanism [ 1], has not been established yet experimentally. This component, however, is crucial for the closure of the SM in a mathematically consistent form, thereby allowing for systematic perturbative expansions and precise predictions in the electroweak sector [ 2]. The search for Higgs bosons [ 3,4] is therefore one of the most important experimental programs in present-day high-energy physics [ 5].
For a light Higgs boson the electroweak interactions can be extended perturbatively up to the grand unification scale. The validity of this extrapolation is backed by the fact that it provides a qualitatively correct estimate of the electroweak mixing angle. Demanding, therefore, the perturbative extrapolation to hold up to the grand unification scale, the mass of the Higgs boson is bounded to be below 180 GeV in the SM [ 6]. An upper limit in the same range is also obtained experimentally by evaluating radiative corrections to electroweak precision observables, i.e. M H ≤ 196 GeV at the 95% confidence level [ 7]. A lower experimental mass limit of 114.4 GeV has been set by the LEP experiments [ 8].
In the near future, the search for Higgs bosons will continue at hadron colliders, the proton-antiproton collider Tevatron [ 3] with a centre-of-mass energy of 2 TeV, followed by the proton-proton collider LHC [ 4] with a centre-of-mass energy of 14 TeV. Various channels can be exploited at hadron colliders to search for Higgs bosons in the intermediate mass range. Among these channels, Higgs-boson radiation off heavy top quarks [ 9,10,11] plays an important rôle: pp/pp → ttH + X. (1.1) Although the expected rate is low at the Tevatron, a sample of a few but very clean events could be observed for Higgs masses below 140 GeV [ 12]. At the LHC, Higgsboson radiation off the top quarks is an important search channel for Higgs masses below ∼ 125 GeV [ 13,14]. Moreover, analyzing the ttH production rate at the LHC can provide valuable information on the top-Higgs Yukawa coupling in relation to other Yukawa couplings [ 15,16]. Precision measurements of the absolute value of the ttH coupling can be completed later at e + e − colliders [ 17,18,19,20,21,22]. Theoretical predictions for cross sections that are based merely on leading-order (LO) QCD are notoriously imprecise. They are plagued by considerable uncertainties owing to the strong dependence on the renormalization and factorization scales, introduced by the QCD coupling and the parton densities. Therefore, higher-order QCD corrections are needed for satisfactory theoretical predictions.
The cross sections for the processes (1.1) are found in next-to-leading order (NLO) QCD by convoluting the cross sections dσ ab of the parton subprocesses ab → ttH(+k) with the parton distributions f a and f b of the initial-state hadrons, both evaluated in next-to-leading order: dσ(h 1 h 2 → ttH + X) = a,b=q,q,g dξ a f h 1 a (ξ a ; µ 2 F ) dξ b f h 2 b (ξ b ; µ 2 F ) dσ ab (ŝ = ξ a ξ b s; µ 2 F , µ 2 R ).
(1.2) In this convolution the partons a and b carry fractions ξ a and ξ b of the original momenta of the hadrons h 1 (= p) and h 2 (=p/p), respectively. The total hadronic centre-of-mass 1 energy is denoted by √ s, whereas the partonic subprocess energy is given by √ŝ . The factorization and renormalization scales are denoted by µ F and µ R , respectively. In general these two scales will be identified.
A brief report on the radiative QCD corrections for the processes (1.1) in NLO QCD has been presented for both the Tevatron and the LHC in Ref. [ 10]. The results of that study were found to be in agreement with a parallel investigation for the Tevatron in Ref. [ 11], which was based exclusively on the dominant partonic process qq → ttH. As expected, in complete calculations of the NLO QCD corrections the scale dependence is reduced significantly, and stable theoretical predictions can be derived for the cross sections.
At the technical level there are two main obstacles in calculating the NLO QCD corrections. First, the one-loop pentagon (5-point) diagrams, involving both soft/collinear singularities and massive particles, have to be calculated in dimensional regularization. To this end, we have developed a calculational technique, exploiting the fact that the singularity structure of the pentagon diagrams can be derived in a universal way. This enables us to switch to regularization schemes that are more suitable for analyzing pentagon diagrams, without losing regularization-scheme-dependent constant terms. Second, the extraction of the singularities in the real part of the NLO QCD corrections has to be performed in a numerically stable way. For this purpose we have adopted the dipole subtraction formalism [ 23] for massive particles [ 24]. [This is the first application of this particular method to a complex NLO QCD calculation with massive quarks.] As an additional cross-check, the partial results for qq → ttH have been recalculated by means of the phase-space slicing method.
Besides the total cross sections, final-state distributions of the Higgs particle and the top quarks in rapidity y and transverse momentum p T will be presented. As in the case of the total cross sections, including the NLO corrections significantly reduces the scale dependence of these distributions. It should be noted that based on this complete calculation, the NLO QCD corrections to arbitrary distributions can be investigated.
The report is divided into four sections. In the next section we describe the calculation of the NLO QCD corrections at the partonic level. In the third section the hadronic results for the SM reactions pp/pp → ttH are presented, comprising the total hadronic cross sections as well as the distributions in rapidity and transverse momentum. The analytical and numerical results of our study are summarized in the last section. Finally, the appendices provide some useful formulas, such as scalar one-loop integrals, results on colour algebra, and splitting functions.
To introduce a compact notation for the SME, the tensors Γ qq {1,µ,µν,...} =vq(p 2 ) {1, γ µ , γ µ γ ν , . . .} u q (p 1 ), Γ tt {1,µ,µν,...} =ū t (p 3 ) {1, γ µ , γ µ γ ν , . . .} vt(p 4 ), (2.4) are defined with obvious notations for the Dirac spinors u q (p 1 ), etc. Furthermore, symbols like Γ p are used as shorthand for the contraction Γ µ p µ . For the qq channel we introduce the following set of SME, The LO amplitude M qq LO of the qq channel is proportional to the colour factor where the first Gell-Mann matrix λ c belongs to the qq state, the second matrix to the tt state. The one-loop amplitudes M qq 1−loop also involve other colour operators, such as 1 ⊗ 1 and (λ a λ b ) ⊗ (λ a λ b ), but they can be reduced to multiples of C qq and tensors that are orthogonal to C qq and therefore do not contribute in the interference colours (M qq LO ) * M qq 1−loop . Thus, C qq is the only relevant colour structure at one loop, and the produced tt pair is in a pure colour-octet state in the qq channel. Both the LO and the (relevant part of the) one-loop amplitude can be written as for the independent crossed SMEM gg i with i > 8. For the gg channel it is convenient to introduce the three colour operators where c n is the colour index of gluon n, and the matrices 1 and λ c act on the tt state. The totally antisymmetric and symmetric SU(3) structure constants f abc and d abc are defined in the usual form. Parts in amplitudes that are proportional to different C gg j do not interfere with each other owing to the orthogonality relations: Obviously C gg 1 corresponds to a colour-singlet state of the tt system, while C gg 2 and C gg 3 describe the two different octet states. Both the LO and the one-loop amplitude can be written as The terms in Eq. (2.18) proportional to M direct correspond to the three "direct graphs", where two parallel gluons couple to the top-quark line (see 3rd graph in Fig. 1). The terms proportional to M crossed result from the direct graphs by crossing, i.e. by interchanging the two gluons (see 4th graph in Fig. 1). The terms proportional to M fusion correspond to the gluon-fusion graphs (see 2nd graph in Fig. 1) which do not receive crossed counterparts. Some comments on the SME and their evaluation ought to be added. The M qq,gg i defined above lead to a unique representation of all LO and one-loop diagrams, as indicated in Eqs. (2.7) and (2.16), if only the Dirac equation, the transversality of the polarization vectors, and the gauge conditions (2.10) are used. The number of SME could actually be reduced further by exploiting discrete symmetries and by using the four-dimensionality of space-time, inducing relations among the SME given above. However, eliminating some of the SME does not simplify the numerical evaluation significantly. The explicit calculation of the SME was carried out in two different ways. First all interference terms spins (M qq,gg i ) * M qq,gg k were calculated by making use of the general polarization sums for Dirac spinors and polarization vectors. In a second independent approach we directly expressed the M qq,gg i in terms of Weyl-van der Waerden spinor products [ 25] for each helicity configuration.
The final step in the calculation of the partonic cross sections and distributions involves the 4-dimensional integration over the final-state three-particle phase space. By introducing an intermediate "tt-state" with virtuality s 34 , one obtains with The labels ab indicate the specific initial state of the partonic process. The solid angle Ω * t of the top quark is defined in the rest frame of the "tt-state", using the direction of flight of the "tt-state" in the initial-state centre-of-mass frame as z-axis. The angle θ CM H denotes the polar angle of the Higgs boson in the initial-state centre-of-mass frame, using the direction of the initial-state beams as z-axis. The corresponding azimuthal angle is trivially integrated to 2π in view of rotational invariance about the initial-state beams. The factor k ab av results from the average over the initial-state spins and colours: where k qg av and k gq av are given for future use. The differential cross section including virtual NLO corrections can be obtained in a straightforward way from Eq. (2.20) by substituting the correct matrix element.

One-loop diagrams and calculational framework
Among the one-loop QCD diagrams, three different gauge-invariant subsets can be distinguished, for both the qq and gg channels. Representative diagrams of these subsets are shown in Fig. 2. Subset (i) comprises all closed-quark-loop graphs where the Higgs boson couples to the external top-quark line, (ii) includes all graphs where the Higgs boson couples to a closed quark loop, and (iii) is formed by all other diagrams, which include at least one gluon in the loop. Counting closed quark loops only once per quark flavour, we get 2, 2, and 25 graphs of the respective subsets in the qq channel; and 6, 20, and 108 graphs in the gg channel. The counting is performed in the Feynman gauge, omitting tadpole-like graphs and external self-energy corrections consistently. Of course, the subsets (i) and (iii) receive counterterm contributions from renormalization, which are considered in Section 2.2.2 in more detail; the second subset is not renormalized in NLO, since Hgg and Hggg couplings do not exist in LO.
In NLO the one-loop amplitudes M qq,gg 1−loop contribute via interference with the corresponding LO amplitudes M qq,gg LO . Using the colour decomposition rules of Section 2.1, the square of the sum of the LO amplitude and its virtual NLO correction reduces to spins, colours The evaluation of the virtual NLO correction was performed in two completely independent calculations, leading, in particular, to two independent computer codes. In the first evaluation, the Feynman graphs have been generated by using FeynArts [ 26]. These amplitudes have subsequently been reduced in terms of SME with the help of Mathematica. Following standard techniques of one-loop calculations [ 27,28,29,30,31], the tensor integrals are algebraically reduced to scalar integrals. The treatment of the so-called "pentagon" diagrams, involving five propagators in the loop, requires particular care and extensions of the standard methods given in Refs. [ 27,28,29,30,31]. The treatment of the pentagons is described in Section 2.2.3 in detail, where the results on the scalar 5-point functions are listed explicitly. The results on the divergent scalar 3-and 4-point functions are collected in Appendix A. Both ultraviolet (UV) and soft/collinear infrared (IR) divergences are treated in dimensional regularization. The renormalization, which removes the UV divergences, is discussed in Section 2.2.2. The explicit structure of the IR singularities occurring in the virtual corrections is given in Section 2.2.5, and the cancellation against their counterparts in the real corrections is described in Section 2.3. The algebraic output of the first calculation has been implemented in Fortran for numerical evaluation. The second calculation basically follows the same strategy, with a conceptual difference that is described in Section 2.2.3 (v). The amplitudes have been generated by hand and reduced with Form [ 32], and the numerical evaluation is carried out in C++.
For later use, we introduce the following definitions for the scalar and tensor integrals in analogy to Refs. [ 27,28,29,30,31], Here +i0 is the infinitesimal imaginary part as required by causality, D = 4 − 2ǫ is the space-time dimension, and µ is the arbitrary reference scale of dimensional regularization, which we may identify with the renormalization scale. The UV and IR divergences appear as poles in ǫ; it is convenient to use the abbreviations for these divergences. [The top-quark mass has been used to define a characteristic energy scale.] If we distinguish between UV and IR divergences, we explicitly write ∆ UV,IR n (µ).

Renormalization
The renormalization is performed in the MS scheme with the top-quark mass defined on shell. The top quark is decoupled from the running of the strong coupling α s (µ). Denoting the bare top-quark mass and the bare strong coupling as m 0 t and α 0 s , respectively, these conditions fix the renormalization parameters in the transformations m 0 t = m t + δm t and α 0 s = α s + δα s as follows, 26) where N f = 5 is the number of light flavours. The divergence ∆ UV 1 (m t ) = Γ(1 + ǫ)(4π) ǫ /ǫ is independent of µ and originates from light-quark and gluon loops. The term proportional to ∆ UV 1 (µ), on the other hand, originates from the top-quark loop in the gluon self-energy that is subtracted at zero-momentum transfer. In this scheme the running of the coupling α s (µ) is generated solely by the finite contributions of the light-quark and gluon loops, while the top-quark contribution is absorbed completely in the renormalization condition and thus decouples effectively. The transformation for m 0 t also fixes the renormalization of the Yukawa coupling, g 0 We renormalize the fields of the gluons, G a,µ , of the light quarks, ψ q , and of the top quark, ψ t , all in the on-shell scheme, i.e. the wave-function renormalization constants δZ G,q,t , defined by the transformations are adjusted to cancel the external self-energy corrections exactly. Distinguishing between divergences of UV and IR origin, these constants can be written as The contributions of the counterterms and the external self-energies to the one-loop matrix elements M qq,gg 1−loop can be derived easily. For the qq channel they explicitly read (2.29) The first term accounts for the external self-energy corrections as well as the counterterms of the strong and Yukawa coupling renormalization, while the remaining terms proportional to δm t come from the renormalization of the top-quark mass. For the gg channel the former corrections have the same structure (proportional to the LO amplitude), but the mass renormalization involves cumbersome expressions.

Pentagon diagrams (i) General strategy
Four different types of pentagon diagrams occur in the calculation; a representative of each type is shown in the last row of Fig. 2. Since the incoming quark q is massless, three different types of 5-point functions are generated which involve one, two, or three massless propagators in the loop. Each of the three types is IR divergent, so that we need explicit results for scalar and tensor 5-point functions in D = 4 dimensions. As known for a long time [ 33], in D = 4 dimensions all 5-point functions can be expressed in terms of 4-point functions, simplifying the calculation considerably. In the following we describe an elegant way how this reduction can be generalized for singular 5-point functions to D = 4 up to terms of order Ø(D − 4) in dimensional regularization.
In order to make use of the reduction of 5-point functions to related 4-point functions, which works in four space-time dimensions, we first translate the dimensionally regularized integral E (D) into another regularization scheme that is defined in D = 4 dimensions. To this end, we first endow all massless propagators in the loop with an infinitesimal mass λ and call the new integral E (mass,D) , which is identical to E (D) if λ = 0. Next we determine a related well-defined integral, denoted E , has a uniquelydetermined finite integrand and is therefore regularization-scheme independent, i.e. the limits D → 4 and λ → 0 commute in this quantity. In total, we have generated in this way the relation these divergences can be controlled by distinguishing the space-time dimensions D IR and D UV for the IR and UV domains explicitly. As shown below, the singular parts are given in terms of 3-point functions, and E (mass,D=4) can be expressed in terms of 4-point functions. Consequently, solving Eq. (2.30) for E (D) , this integral is expressed (up to irrelevant terms indicated by the dots) in terms of 3-point and 4-point functions.
(ii) Determination of singular subintegrals in 5-point functions Having described the general strategy, we now determine the singular terms in the relevant 5-point functions. The following derivation is valid for the integrals E (mass,D) , which are defined for arbitrary D in the vicinity of D = 4 and which include infinitesimal mass regulators λ for internal massless particles. The final formulae given for E sing below are thus valid for both E . We do not write the regulator λ explicitly, since the relations for E sing also hold for other choices of mass regularizations. 1 We first consider the case of one vanishing internal mass, as shown in the second diagram of the last row in Fig. 2, for example. There is only a purely soft (logarithmic) divergence which is located at q → 0 in the integrand if the vanishing mass is m 0 = 0. In this soft limit the numerator of the integral and the denominators d 2 and d 3 are regular. As a result, the singularity structure of the integral remains unchanged if we set q = 0 in the numerator and the denominators d 2 and d 3 , since the divergence is logarithmic. This, in particular, shows that all such tensor integrals E µ... are finite, and only the scalar integral E 0 develops a soft singularity, which is already contained in a scalar 3-point function. There are two different types of 5-point integrals with one internal massless line, corresponding to the second diagram of the last row in Fig. 2 and the additional diagram generated by shifting the Higgs line along the top-quark line. Explicitly the singular parts of the corresponding integrals are given by Next we consider the case with two internal massless lines, as generated, for example, in the third diagram of the last row in Fig. 2. In this case, the denominators d i are given by where the imaginary i0 terms are suppressed in the notation. The integral becomes singular in the collinear limit q → xp 1 with arbitrary x, and in the soft limits q → 0 and q → −p 1 , which are included in the collinear case for the limits x → 0 and x → −1, respectively. The non-singular denominators d 2,3,4 behave in the collinear limit as This leads to the following decomposition of the singular part of the 5-point integral, This relation is valid for all tensors; in contrast to the case with only one vanishing mass, in this case also tensor integrals become singular. Note that E sing does not only include the singularities but also regular terms in a well-defined way, which is important when switching from one regularization scheme to the other.
Finally, we consider the case with three internal masses zero, as encountered, e.g., in the first and last diagram of the last row in Fig. 2. Here three soft and two collinear limits give rise to divergences. Identifying m 0 = m 1 = m 4 = 0, the collinear limits are realized by setting q → xp 1 or q → yp 2 , and the three soft limits are reached for particular values of x and/or y (i.e. x = y = 0, x = −1, or y = 1). The decomposition for the two collinear limits can be worked out in analogy to the previous case. When combining the two limits, however, double-counting of the soft singularity at x = y = 0 has to be avoided. The decomposition of the singular structure of the 5-point integral finally reads which is again valid for all tensors.
(iii) The scalar 5-point function In D = 4 space-time dimensions the scalar 5-point function E 0 can be expressed [ 33] in terms of the five related 4-point functions D 0 (i), obtained from E 0 by omitting the denominator d i , as defined in Eq. (2.24). Using the shorthand notations where the scalar integrals are defined for D = 4. Solving for E (mass,D=4) 0 we get where Y = (Y ij ), and Y i is obtained from Y by replacing all entries in the ith column by 1. Although Eq. (2.39) was derived for D = 4, it is valid (in all cases considered here) even for arbitrary D up to terms of order Ø(D − 4): This generalized relation simply follows from the result derived in (ii) that the same linear relation between E sing and 3-point subintegrals holds in all mass/dimensional regularization schemes.

(iv) Tensor 5-point functions
We evaluate all tensor integrals by calculating the coefficients to the Lorentz covariants that span the tensor. For the pentagons, tensor integrals E µν... appear up to rank 4. In practice, however, these integrals are needed only up to rank 3, since one momentum factor q in the numerator of a top-quark propagator can always be cancelled against a propagator denominator. The covariant decomposition of the relevant 5-point tensor functions reads Since the tensors E µν... are symmetric, the tensor coefficients E ij and E ijk are symmetric in the indices. The coefficients E i , E ij , etc., can be recursively reduced to the scalar integrals E 0 , D 0 , etc., by an algorithm known as Passarino-Veltman reduction, as described in Refs. [ 28,30], or by special techniques [ 34,35] developed for higher n-point tensor integrals. The resulting expressions are cumbersome and do not reveal any new insight; therefore we do not list the expressions explicitly. Note that the metric tensor appearing in Eq. (2.42) is D-dimensional. In contrast to the D-dimensional metric tensor g αβ , its four-dimensional part, g (D=4) αβ , is spanned by the linearly independent base vectors q i : where Z is the Gram matrix of the momenta q i . The difference ∆g αβ = g αβ − g (D=4) αβ is, thus, of Ø(D − 4). It is interesting to notice that the terms proportional to ∆g αβ do not contribute in our case. Since ∆g αβ is an object of Ø(D − 4) it could only contribute if an accompanying tensor coefficient, such as E 00 , is divergent. The explicit tensor reduction to scalar integrals, however, reveals that all such coefficients are finite. This fact can also be derived from the following argument. In tensor integrals the only covariants that come with divergent coefficients are built up by the singular regions in momentum space. Therefore, a collinear divergence occurring for q → xp can only show up in tensors built up by the momentum p alone. For purely soft (logarithmic) divergences occurring at q → 0 power-counting shows that only the scalar integral is divergent; shifting the integration momentum by p, i.e. q → q + p, only leads to divergent coefficients related to tensors of p alone. UV divergences could occur in covariants containing ∆g αβ ; however, all relevant 5-point integrals are UV-finite, which is easily proven by power-counting.

(v) Alternatives in practice
The results on 5-point functions described above admit two conceptually different ways to evaluate the pentagon diagrams. We followed both ways, as described below, and we found perfect numerical agreement.
In the first approach the scalar integral E 0 is evaluated in D dimensions using Eq. (2.41) and also the reduction of the tensor coefficients defined in Eq. (2.42) is performed consistently in D = 4. In this approach all divergent scalar 3-and 4-point integrals C 0 and D 0 are needed in dimensional regularization. We have listed them in Appendix A.
Alternatively one can make consistent use of Eq. (2.31) and its analogous version for 4-point functions, in order to replace all D-dimensional 4-and 5-point integrals by fourdimensional integrals and D-dimensional 3-point functions that contain the IR (soft and collinear) singularities. The reduction of the 4-and 5-point tensor coefficients then works in four IR dimensions 2 . An advantage of this procedure is the easy analytical control of all IR singularities which are contained in 3-point functions only. In this approach all divergent D 0 functions are needed in the mass-regularization scheme, and all divergent C 0 functions in both the dimensional and the mass-regularization schemes. The massregularized functions can be derived from the results of Appendix A. Another benefit of evaluating pentagon diagrams in four dimensions is the direct application of the tensor decomposition of Ref. [ 35]. This procedure consistently avoids leading Gram determinants (i.e. those of four momenta) that are potential sources of numerical instabilities, as described in the next section.

Numerical stabilization
As explained above, all tensor integrals occurring in the calculation can be algebraically reduced to scalar integrals following the well-known Passarino-Veltman algorithm [ 28]. This procedure is based on a decomposition of the tensor integral into a complete set of covariants that are spanned by the metric tensor and all independent momenta involved in the integral. The coefficients of the covariants are obtained by inverting a set of linear equations, which is non-singular as long as the momenta q i spanning the covariants are linearly independent, i.e. as long as their Gram determinant G = det(q i q j ) is non-zero. Each tensor rank adds a factor G −1 to the expression of the tensor integral. Near the boundary of phase space, some determinants G become arbitrarily small, since some momenta become linearly dependent at the boundary. The huge factors G −n are balanced by cancellations between different tensor coefficients that are calculated in terms of scalar integrals; the whole virtual correction remains well behaved at the phase-space boundary. These cancellations lead to severe numerical instabilities in the evaluation of the correction. The instabilities are most pronounced in the pentagon diagrams, which do not only possess the highest tensor ranks, but also involve additional sources of cancellations. For instance, the decomposition (2.43) of the four-dimensional metric tensor also leads to a Gram determinant in the denominator, and the reduction (2.39) of the scalar five-point function E 0 in terms of five four-point functions involves Caley determinants det(Y ... ) that can become very small.
If the virtual correction is not evaluated with particular care at the phase-space boundary, a phase-space point will soon be selected in the Monte Carlo integration where the result for the correction is totally wrong. The wrong result will typically be a large number owing to incomplete cancellations in the evaluation. Such unphysically large contributions occur more and more frequently in the adaptive integration procedure and they thus destroy the result of the phase-space integration completely .
When using the Passarino-Veltman procedure, we avoid the numerical instabilities by extrapolating the virtual correction from the numerically safe inner region of phase space to the numerically delicate boundary. To this end, we first divide the full phase space into two regions, one where the correction is evaluated directly and another where the correction is obtained from an extrapolation out of the first region.
The extrapolation procedure is applied to the full one-loop correction, since any subset of diagrams (such as the pentagons alone) in general shows stronger variations over the phase space than the full correction, owing to cancellations between different diagram types. Therefore, we have to inspect the whole boundary of phase space, which is characterized by the zero of the Gram determinant G of the four momenta p 1 , . . . , p 4 , where κ is defined in Eq. (2.21) and p 5 is fixed by momentum conservation, p 1 + p 2 = p 3 + p 4 + p 5 . The three-momenta p i and the angles θ 1 , χ are defined in the partonic centre-of-mass frame. The angle θ 1 is the angle between p 3 and the beam axis, and χ denotes the (azimuthal) angle between the p 4 -p 5 plane and the plane spanned by p 3 and the beam axis. We identify the "safe" region by demanding G/ŝ 4 as well as sin χ to be larger than certain numerical cut values. The extrapolation into the "problematic" region is done as follows: 1. Regarding the integrand as a function of the variable χ, thereby keeping p 3 and |p 4 | fixed, we expand it in terms of a Fourier series. As a set of basis functions we took trigonometric functions or Legendre polynomials. The first few coefficients of the series are determined by numerical integration, for which we took Simpson or Gaussian quadrature. If the extrapolation result is close to the directly calculated result for the correction, the direct result is taken over.
2. If the agreement between direct and extrapolated result is too bad, the first step is repeated again with twice as many coefficients in the Fourier series and twice as many points in the related integration for the coefficients. The first step is repeated at most twice. If the directly calculated result is still not in good agreement with the extrapolation, but instead the various extrapolations look convergent, the most precise extrapolated result is taken over.
3. If the procedure still has not stopped yet, step one is modified by changing the extrapolation direction by varying also p 3 and |p 4 | during the change in χ.
Step two is carried out as above.
4. The modified extrapolation of step three is repeated if still no success has been reached. Finally, if none of the extrapolations looks convincing, the integration point is discarded, but this does not happen very often.
The extrapolation procedure solves the problem of numerical instabilities at the expense of computing time. Of course, the described procedure has to be optimized by selecting appropriate cut values on G/ŝ 4 and sin χ and by choosing appropriate techniques for the calculation of the Fourier coefficients. Nevertheless the final phase-space integration has to be stable against any change in the extrapolation procedure. We have checked this stability by changing the cut parameters and by choosing different systems of orthogonal functions as well as different integration techniques, as mentioned above. In addition to adopting the usual Passarino-Veltman reduction of 5-point tensor coefficients, we have applied the method of Ref. [ 35] where the strategy [ 33] for scalar 5-point functions is generalized and 5-point tensor coefficients are reduced to 4-point integrals directly. In this approach inverse Gram determinants of four momenta are avoided completely. Applying this alternative renders the virtual correction near the phase-space boundary much more stable than in the usual Passarino-Veltman approach. The results obtained by the two methods mutually agree with each other.

One-loop soft and collinear divergences
There are two sources of IR (soft and collinear) singularities in the one-loop diagrams: the self-energy corrections of the external fields and various kinds of 3-point subintegrals whenever a (nearly on-shell) gluon is exchanged between external lines. The former singularities can be read off from the results of Section 2.2.2, and the strategy to determine the latter singularities has already been explained in Section 2.2.3 (v). Here we list the explicit results: The functions M ... are given in Eq. (2.19), and the divergent factors f qq ij and f gg ij read These results are in agreement with Ref. [ 36], where the general singularity structure of QCD amplitudes has been presented. The agreement becomes apparent after rewriting Eqs. (2.45) and (2.46) in the following generic form which is valid for both the qq and gg channels, Here T i are the colour operators defined in Refs. [ 23,24,36], and the symbol ⊗ denotes the colour correlations. Since we will make repeated use of this notation in the remainder of our report, a short description can be found in Appendix B.

Real corrections
The real corrections at relative order O(α s ) include processes with real-gluon radiation as well as ttH production in gluon-(anti)quark scattering. After a brief description of the parton reactions and their evaluation, we will discuss in some detail the isolation of IR divergences for soft and collinear parton configurations as well as the treatment of collinear initial-state singularities by means of mass factorization.

Parton processes and their evaluation
The evaluation of the O(α 3 s g 2 ttH ) cross section requires the calculation of the gluon bremsstrahlung processes q(p 1 ) +q(p 2 ) → t(p 3 ) +t(p 4 ) + H(p 5 ) + g(p 6 ), g(p 1 ) + g(p 2 ) → t(p 3 ) +t(p 4 ) + H(p 5 ) + g(p 6 ), (2.49) and the gluon-(anti)quark scattering reactions A representative set of Feynman diagrams, contributing to the real corrections, is given in Fig. 3. Helicity amplitudes for the processes (2.49, 2.50) have been generated and evaluated using the program packages MadGraph [ 37] and HELAS [ 38]. The result has been verified by an independent calculation based on standard trace techniques.

Soft and collinear divergences
The 2 → 4 cross sections contain singularities in the limits where the splitting g → gg, g → qq, q → qg, orq →qg becomes soft or collinear. We have used the dipole subtraction method [ 23], as formulated in Ref. [ 24] for massive QCD partons, to extract the singular part of the real cross section in D dimensions. The dipole subtraction method is based on the fact that soft and collinear limits of a real-emission matrix element can be expressed by convoluting process-independent splitting functions with the LO matrix element. By subtracting these splitting-function expressions as counterterms from the squared real matrix element, the IR (soft and collinear) singularities cancel and the phasespace integration can be performed numerically in four dimensions. The counterterms are constructed such that they can be integrated analytically in D dimensions over the phase space of the extra emitted parton, leading to poles in ǫ = (4 − D)/2. When the integrated subtraction counterterms and the virtual cross sections are combined, the soft singularities cancel. The remaining initial-state collinear divergences have to be absorbed into a redefinition of the parton distribution functions at NLO. Below we will describe the subtraction formalism in more detail and we will give explicit results for the soft and collinear contributions.
The NLO contribution of the cross section can symbolically be written as [ 23,24] σ ab This symbolic expression is valid under the assumption that no massless particles are identified in the final state, i.e. the additional gluon or light (anti)quark is integrated out. The indices a and b denote the flavours of the incoming partons, i.e. ab = qq, gg, qg, gq in the processes (2.49,2.50), and µ F is the factorization scale at which the redefinition of the parton distribution functions is performed at NLO. The notation for the integrals indicates that the real-emission contribution (labelled by R) involves four final-state particles [one gluon or light (anti)quark more than in LO], while the virtual contribution (labelled by V) has LO three-particle kinematics. The dipole subtraction counterterm is given in terms of universal splitting functions, dV dipole (see Ref. [ 24] for explicit expressions), and an appropriate colour, spin, and flavour projection of the LO cross section, dσ LO . The symbol ⊗ denotes such possible colour, spin, and flavour correlations. In particular, in the sum over "dipoles" the flavour correlations can change the incoming partons a or b in dσ ab LO . The final state of dσ LO is fixed to ttH, since no light (anti)quark can be radiated from the top quarks. As far as the four-particle phase space is concerned, each individual dipole counterterm is evaluated for a specific set of auxiliary momentap 1 , . . . ,p 5 . This involves combining the emitting particle (with momentum p j ) and the emitted light particle (with momentum p 6 ) into a single parent particle called "emitter" (with momentump j ). In order to balance the total momentum, one of the other external particles (with momentum p k ) is replaced by a "spectator" particle (with the same quantum numbers, but modified momentump k ). Apart from the situation where both emitter and spectator are incoming partons, the momenta of the remaining external particles are not altered 3 . This phase-space mapping is constructed in such a way that -the soft and (quasi-)collinear momentum configurations are reproduced point-wise; -both the emitter and spectator are on their respective mass shells; -momentum conservation is preserved; -the phase space can be factorized into a hard-scattering phase space and a singleparticle phase-space factor corresponding to p 6 .
Many of these phase-space mappings either have an initial-state emitter or an initialstate spectator. In those cases the corresponding initial-state momentum is rescaled by a longitudinal momentum fraction x, resulting in a boosted reference frame for the finalstate momenta and a phase-space factorization that involves a convolution over the boost parameter x. For explicit phase-space mappings we refer to the literature [ 23,24]. The readded subtraction counterterm, integrated over the phase space of the emitted extra parton, can be split into two parts. The first part, involving the universal insertion operator I ab , has LO kinematics and contains all the poles in ǫ that are necessary to cancel the soft singularities in dσ ab V . These so-called endpoint contributions are obtained by integrating the dipoles over the complete phase space of the emitted extra parton and by using the (. . .) + prescription to extract the singular terms related to the soft endpoint (x → 1) of the x-integration of the dipoles. The second part is a finite remainder that is left after factorization of initial-state collinear singularities into the parton distribution functions at NLO. It involves LO cross sections with boosted three-particle kinematics, x-dependent structure functions Γ(x, µ 2 F ), and an additional one-dimensional integration with respect to the longitudinal momentum fraction x of an incoming parton. The structure functions Γ(x, µ 2 F ) depend on the factorization scheme; only the terms proportional to ln(µ 2 F ) are factorization-scheme independent. The summations over the flavours a ′ and b ′ represent the afore-mentioned flavour correlations. Non-zero contributions are obtained if both dσ a ′ b LO and Γ a,a ′ ttH+b are non-zero, or if both dσ ab ′ LO and Γ b,b ′ ttH+a are non-zero. In the remainder of this subsection we give the explicit expressions for the readded subtraction counterterms, which are based on the results of Ref. [ 24]. The explicit construction of the dipole terms, which is not described here, can also be found there. We start with the insertion operator I ab , which depends on the colour charges, momenta, and masses of the final-state particles in the LO matrix element: for i = 1, 2 and j = 3, 4. (2.53) The variable β s 34 is defined in Eq. (2.47) and ρ is given by ρ = (1 − β s 34 )/(1 + β s 34 ). The dilogarithm Li 2 (x) is defined in the usual way as (2.54) The Casimir operators T 2 a and anomalous dimensions γ a read explicitly Inserting all this in Eq. (2.53), we observe indeed that the singular ∆ IR 2 (µ) and ∆ IR 1 (µ) terms in the functions J ab ij cancel against the singular terms in the coefficients f ab ij given in Eq. (2.47).
Finally, we give a general expression for the structure function Γ a,a ′ ttH+b (x, µ 2 F ) in the MS factorization scheme: . (2.56) The various Altarelli-Parisi splitting functions P aa ′ (x), P aa ′ reg (x), andP ′ aa ′ (x) can be found in Appendix C, and the kinematical variable µ 2 j1 is given by wherep j is the auxiliary momentum of particle j in the boosted reference frame, while p 1 is the original initial-state momentum of the incoming parton a. In our convention the (. . .) + prescription is defined by where the label (x) of µ (x) j1 indicates thatp j belongs to the phase space Φ(x) of the xboosted reference frame and therefore µ (1) j1 = m t / √ 2p 1 p j . In fact, this is the same as saying that µ j1 should be kept fixed during the x-integration, e.g. by using it as a phase-space integration variable in dΦ(x). The expression for the structure function Γ b,b ′ ttH+a (x, µ 2 F ) can be obtained from Eqs. (2.56) and (2.57) by interchanging (a, a ′ , µ j1 , p 1 ) and (b, b ′ , µ j2 , p 2 ).
The formulas presented above are spelled out for total cross sections, but they can also be used to evaluate any kind of IR-safe differential distributions by filling histograms during the Monte Carlo integration. In this context, the only subtlety in the dipole subtraction approach is due to the fact that contributions of different dipole terms in general contribute to different histogram bins. The IR safety of an observable guarantees that singular terms of the original integrand and their corresponding counterterms contribute to the same bin in the soft and collinear limits. 4

Slicing method for the qq channel
In the phase-space slicing approach (see Ref. [ 39] for a review) the soft and collinear regions are excluded from phase space by appropriate phase-space cuts, generically denoted by ∆ cut . The singular integration over these regions is then carried out analytically using factorization properties. The full real correction is obtained by the sum of these two contributions in the limit ∆ cut → 0. In practice, a plateau is required for small but finite values of ∆ cut .
In the following paragraphs we describe the so-called two-cutoff method for the qq channel, where the soft region for gluon bremsstrahlung is defined by a cut ∆E on the gluon energy and the collinear region for gluon radiation off the initial state by a cut ∆θ on the emission angle.

(i) Soft-gluon region
Denoting the momentum of the radiated gluon by k, the soft region is defined by where k 0 is the gluon energy in the partonic centre-of-mass frame. In the soft region, gluon radiation is described by an eikonal current (see e.g. Ref. [ 23]), and the contribution to the bremsstrahlung cross section is given by where g ij denote the soft integrals The integrals are obtained as follows: g 12 is calculated easily; g 34 can be taken over from Ref. [ 31], where this integral is given for an infinitesimal gluon mass m g that translates into 1/ǫ via the substitution ln(m 2 g ) → (4πµ 2 ) ǫ Γ(1 + ǫ)/ǫ; the remaining g ij can be derived using the auxiliary integrals (C.20) and (C.25) of Ref. [ 40]. The results are . (2.64) The 1/ǫ 2 poles, i.e. the ∆ IR 2 (µ) terms, of the functions g ij are the same as in the coefficients f qq ij given in Eq. (2.47). This implies that these terms cancel in the sum of virtual and soft corrections, as it should be.

(ii) Collinear gluon emission from the initial state
The region of collinear gluon emission from the initial state is defined by in the partonic centre-of-mass frame. According to the QCD factorization formula (see e.g. Ref. [ 23]), the contribution of collinear gluon radiation off parton i with i = 1, 2 reads where the Casimir operator C F = 4/3 has been defined in Eq. (2.55) and the Altarelli-Parisi splitting function P qq (x) can be found in Appendix C. These results can, for instance, be easily derived by evaluating the x and v integrals in Section 5.5 of Ref. [ 23] with the ∆E and ∆θ cuts. From the results given above it is easy to verify that all divergences ∆ IR 1 (µ) cancel in the sum of the virtual corrections (see Section 2.2.5), the soft corrections dσ qq soft , and the part of the collinear corrections dσ qq coll with the LO kinematics, i.e. the part proportional to dσ qq LO (p i ). The remaining P qq (x) ∆ IR 1 (µ) term is absorbed into the redefined parton densities.

(iii) Independent check within the mass-regularization scheme
The sum of dσ qq soft + dσ qq coll has been checked by deriving all formulae for infinitesimal quark masses. Apart from g 34 this calculation is completely independent to the procedure in paragraphs (i) and (ii). These are the steps to be followed: 1. Derive dσ qq soft and dσ qq coll for infinitesimal quark masses for q andq; see Ref. [ 31] for dσ qq soft , and Ref. [ 41], in particular Eq. (26), for dσ qq coll .
2. The translation of the continuum part of dσ qq coll [the part with P qq (x) and the (. . .) + prescription] to the dimensional regularization scheme can be performed by comparing the renormalization of the parton distribution functions within the two regularization schemes. Specifically, Eq. (6.6) of Ref. [ 23] has to be compared with Eq. (23) of Ref. [ 42] (see also Ref. [ 43]).
3. The translation of the endpoint part of dσ qq soft + dσ qq coll [this is the part proportional to dσ qq LO (p i )] is described in Eqs. (23) and (41) of Ref. [ 36].

Leading threshold corrections
In view of the complexity of the NLO QCD corrections, as detailed in the previous subsections, any additional cross-check is welcome. The partonic energy regime that is most suited for performing such cross-checks is the threshold region, Apart from providing explicit cross-checks, the leading threshold corrections can also serve as a simple, intuitive tool for understanding the differences between the hadronic NLO corrections at the Tevatron and the LHC. This should not be overrated, though. In ttH production, with its massive three-particle phase space, the influence of the threshold region is substantially less than for processes with massive two-particle final states, like tt production (see Section 3.1).

The Sommerfeld rescattering correction
For energies near the threshold the virtual QCD corrections are enhanced by Coulombic gluon exchange between the top and antitop quark in the final state. The corresponding generic diagram is shown in Fig. 4. When the gluon momentum tends to zero a so-called Coulomb singularity is generated. In this limit the contribution of the singularity to the matrix element, denoted by M Coul , is given by and it increases inversely proportional to the t/t velocity β s 34 in the tt centre-of-mass frame, as defined in Eq. (2.47). This is the well-known Sommerfeld rescattering correction [ 44]. Upon averaging 1/β s 34 over the phase space near threshold, we obtain the following singular contribution to the partonic cross section The coefficient β is proportional to the maximal relative velocity of the tt system. The rescattering coefficient C Coul depends on the production mechanism of the tt pair, In e + e − annihilation the tt pair is generated in a colour-singlet state in which the quark and antiquark attract each other, and with C e + e − Coul = + 4/3 the correction is positive. This leads to a strong enhancement of the e + e − → ttH annihilation cross section near threshold 26 [ 18]. By contrast, in qq annihilation the tt pair is generated in a colour-octet state in which the force is repulsive, and with C qq Coul = − 1/6 the correction is negative and relatively small. This configuration is realized in the dominant channel at the Tevatron, qq → ttH. Finally, in gg fusion the various colour channels contribute differently, and the relative weighting depends on the ratio M 2 H /4m 2 t . For not too large M H the rescattering coefficient C gg Coul is moderately positive, and for M 2 H /4m 2 t ≪ 1 the rescattering coefficient approaches C gg Coul (M H = 0) = 11/42 for gg → tt (see Refs. [ 40,45]). In the close vicinity of threshold, our numerical results for the NLO partonic cross sections indeed reproduce these Coulomb singularities (see Section 3.1).

Threshold logarithms from soft-gluon radiation
Owing to the lack of available phase space near threshold, enhanced real-radiation effects from the initial state (ISR) and the final state (FSR) are only noticeable in the vicinity of soft IR poles. In that case the phase-space suppression is compensated by the inverse dependence on the energy of the emitted massless particle, resulting in an integrated correction that is logarithmic in the phase-space volume -in contrast to the Rutherford pole in the Sommerfeld correction. The calculation of these so-called threshold logarithms simply amounts to a convolution of the LO cross section with the soft-pole parts [∝ 1/(1 − x)] of the dipole functions for real-gluon emission.
The LO matrix elements of the reactions qq, gg → ttH are basically constant near threshold. This leads to the following threshold behaviour of the LO cross section [using Eqs. (2.20) and (2.68)]: where σ ab LO represents the part of the LO cross section that remains constant near threshold. The strong suppression of the LO cross section near threshold ∝ β 4 is a consequence of the massive three-particle phase space.
The threshold logarithms at NLO are given by (2.72) The function Γ log (x) comprises all IR-sensitive parts of the structure functions Γ a,a ′ ttH+b and Γ b,b ′ ttH+a in the limit β → 0: (2.73) where µ F denotes the factorization scale. Exploiting charge conservation and keeping only the terms involving ln β, one obtains finally and (2.76) The contributions from final-state gluon radiation exhibit a behaviour in colour space that is parallel to the Coulomb singularities. In qq annihilation the tt pair is generated in a colour-octet state, resulting in an overall colour charge C qq FSR = 3. In the case of gg fusion both colour-octet and colour-singlet states are superimposed. However, soft gluons cannot be emitted from colour-singlet tt final states near threshold owing to the overall vanishing colour charge and the long wavelength of the emitted IR gluons -as evident from e + e − annihilation [ 46]. Therefore only the colour-octet state receives a logarithmic correction, causing the overall coefficient C gg FSR to depend on the ratio M 2 H /4m 2 t . In total, the threshold logarithms give rise to a positive NLO correction near threshold, as long as the factorization scale is not chosen too small. Furthermore we observe that the corrections in the gg channel are larger than the corrections in the qq channel by roughly a factor of two. A threshold analysis of our numerical results for the partonic NLO cross sections indeed confirms the above-given threshold logarithms (see Section 3.1).

Numerical Results
In this section we present the NLO QCD results for the total partonic cross sections and for the final hadronic cross sections at the Tevatron and at the LHC, including the differential distributions in transverse momentum and rapidity of the Higgs boson and of the top quarks. As discussed previously, renormalization and factorization have been performed in the MS scheme with the top-quark mass defined on-shell. The top quark is decoupled from the running of the strong coupling α s (µ) and from the evolution of the parton densities. For the evaluation of the hadronic cross sections we have adopted the MRST [ 47] parton densities at LO and NLO, corresponding to the QCD parameters Λ LO 5 = 167 MeV and Λ MS 5 = 239 MeV at the one-and two-loop level of α s (µ), respectively. The strength of the SM Yukawa coupling is fixed by g ttH = m t /v, where v = ( √ 2G F ) −1/2 is the vacuum-expectation value of the Higgs field and G F = 1.16639×10 −5 GeV −2 . The top-quark mass is set to m t = 174 GeV.

Partonic cross sections
To explore the NLO QCD results at the parton level for the production of ttH in quark and gluon collisions, we introduce the scaling functions f ab (η), ab (η) +f (1) ab (η) ln The centre-of-mass energy of the partonic reaction, √ŝ , is expressed in terms of the quantity η =ŝ/4µ 2 0 − 1, which is better suited for analyzing the scaling functions in the various regions of interest. The reference scale µ 0 = m t + M H /2 has been set to half the partonic threshold energy. The index pair a, b indicates the partonic q,q, g initial state of the reaction ab → ttH + X. We have identified the renormalization scale with the factorization scale: µ R = µ F = µ. The scaling functions are divided into the LO term f (0) , the scale-independent NLO term f (1) , and the scale-dependent NLO contribution f (1) ln(µ 2 /µ 2 0 ). The scaling functions are displayed in Fig. 5 for the quark-antiquark, gluon-gluon, and (anti)quark-gluon channels. The Higgs-boson mass is set to M H = 120 GeV. As described in Section 2.4, the parton cross section near the production threshold (η ≈ β 2 ≪ 1) is dominated by soft-gluon effects which can be calculated analytically. The inserts in Fig. 5 show the ratio of the threshold approximation, σ thres = σ Coul +σ log , Eqs. (2.68) and (2.74), and the numerically integrated, complete NLO scaling functions. It is evident that the complete NLO calculation approaches the threshold approximation as η → 0. However, because of the β 4 suppression of the massive three-particle phase space, see Eq. (2.71), the impact of the Coulomb singularity and the threshold logarithms from soft-gluon radiation is suppressed. This behaviour is significantly different from processes with massive twoparticle final states, where the phase space near threshold scales ∝ β and where soft-gluon radiation provides a source of large NLO corrections.
At high energies the NLO partonic cross sections in the gluon-gluon and (anti)quarkgluon channels approach non-zero limits (plateaus) asymptotically, resulting in a ratio f (1) gg /f (1) gq → 2C A /C F = 9/2. This has to be contrasted with the effective scaling behaviour ∝ 1/ŝ of the LO cross section in the invariant energy range of the partons far above threshold, √ŝ ≫ µ 0 (but at the same time ln √ŝ / ≫ ln m t , so that t-channel top-exchange diagrams are not yet logarithmically enhanced). The enhancement of the NLO diagrams is caused by the exchange of soft, nearly on-shell, space-like gluons, associated with the radiation of a hard parton from the initial state (see Fig. 6). This gives rise to a radiativereturn phenomenon, causing the gg-fusion subprocess induced by the virtual gluon to effectively take place at much lower energies, i.e. √ŝ eff = O(µ 0 ). In analogy with the Drell-Yan process, the cross section for the qq channel vanishes asymptotically in leading and next-to-leading order. Therefore we expect the hadronic cross sections at the Tevatron to be affected substantially less by these high-energy plateaus than the corresponding cross sections at the LHC.
The average CM energy √ŝ for the parton processes at the Tevatron and LHC colliders is approximately 650 GeV and 850 GeV, respectively, corresponding to η ≈ 1 − 2. From Fig. 5 we can conclude that in this range of η the NLO corrections at the central scale µ 0 are moderately negative for the quark-antiquark channel, which dominates the cross section at the Tevatron, and moderately positive for the gluon-gluon channel, which is the prominent production mechanism at the LHC. The (anti)quark-gluon channel is strongly suppressed at both the Tevatron and the LHC.

Hadronic cross sections
The numerical analyses of the hadronic cross sections have been performed for the Fermilab Tevatron pp collider with a centre-of-mass energy of √ s = 2 TeV, and for the

Total cross section and scale dependence
The total cross section at the Tevatron is displayed in Figs. 7 and 8. Representative results are listed in Table 1. For a Higgs mass between 100 and 150 GeV, the cross section varies between about 10 fb and 1 fb, if the central value µ → µ 0 = (2m t + M H )/2 is chosen for the renormalization and factorization scales. In NLO the theoretical prediction is remarkably stable with very little variation for µ between ∼ µ 0 /3 and ∼ 3µ 0 , in contrast with the Born approximation, for which the production cross section changes by more than a factor of 2 within the same interval. It has been checked that no accidental compensation of scale dependences is introduced by identifying renormalization and factorization scales. Although the cross section at the Tevatron is strongly dominated by the qq annihilation channel for scales µ ∼ µ 0 , the proper study of the scale dependence requires the inclusion of the gg, qg, and gq channels. If µ is chosen too low, large logarithmic corrections spoil the convergence of perturbation theory, and the NLO cross section would even turn negative for µ < ∼ µ 0 /5. As apparent from Fig. 8, the K factor, K = σ NLO /σ LO with all quantities calculated consistently in lowest and next-to-leading order, varies from ∼ 0.8 at the central scale µ = µ 0 to ∼ 1.0 at the threshold scale µ = 2µ 0 . The K factor is nearly independent of M H in the relevant Higgs-mass range (see Table 1). As argued in Ref. [ 10], the K factor below unity is in qualitative agreement with results obtained in the fragmentation picture proposed in Ref. [ 48].
The cross section at the LHC is dominated by the gluon-gluon process, which receives positive NLO corrections in the vicinity of the central renormalization/factorization scale µ 0 . At µ 0 we obtain K ∼ 1.2, increasing to ∼ 1.4 at the threshold value µ = 2µ 0 , see Figs. 9 and 10. Again these values are nearly independent of M H in the relevant Higgs-mass range and compatible at the qualitative level with estimates obtained in the fragmentation picture [ 48]. The NLO corrections significantly reduce the renormalization and factorization scale dependence and stabilize the theoretical predictions for the cross section at the LHC.
As argued in Section 3.1, the impact of the Coulomb singularity and the threshold logarithms from soft-gluon radiation is suppressed because of the β 4 behaviour of the massive three-particle phase space near threshold. Nevertheless the threshold region gives us insight into the different colour factors that determine the NLO corrections at the   Tevatron and the LHC. Just like the partonic high-energy region, where the fragmentation picture applies and the high-energy plateaus dominate the NLO corrections, also the lowenergy threshold region hints at K factors at the LHC that are substantially larger than the ones at the Tevatron (see Section 2.4).
We have also studied the uncertainty in the cross-section prediction due to the error in the parametrization of the parton densities. To this end we have compared the NLO cross section evaluated using the default MRST [ 47] parametrization with the cross section evaluated using the CTEQ6 [ 49] parametrization. The results are collected in the rightmost column of Table 1. Both MRST and CTEQ include a systematic treatment of partondistribution-function uncertainties and provide the means for a quantitative estimate of the corresponding uncertainties in the cross sections. Using the parton-distribution-error packages [ 47,49], we find that the uncertainty due to the parametrization of the parton densities in NLO is less than ≈ 5% at the Tevatron, where the quark-antiquark initial states dominate, and less than ≈ 10% at the LHC, where gluon-gluon initial states are more prominent.

Differential distributions
In Figs. 11-15 the normalized transverse-momentum and rapidity distributions of the Higgs boson and the top quarks [t andt not discriminated] are shown for pp → ttH at the Tevatron in LO and NLO approximation. For the transverse-momentum distributions we choose the transverse mass as the default renormalization and factorization scale, i.e. µ 2 = p 2 T,H + M 2 H and µ 2 = p 2 T,t/t + m 2 t for Higgs-boson and top/antitop distributions, respectively, which provides a more natural choice for large transverse momenta. 5 The scale dependence of the distributions in rapidity and transverse momentum of the Higgs boson and top quarks is greatly reduced at NLO. The improvement is exemplified in Fig. 11 which compares the transverse-momentum distribution of the Higgs boson, once calculated with the renormalization and factorization scales set to the transverse mass, and once with both scales fixed to µ 0 . We observe differences of up to 30% at LO, while the shape of the p T,H distribution at NLO is practically identical for the two choices of scales.
The ratio of the (normalized) NLO and LO transverse-momentum distributions of the Higgs boson, displayed in the insert in Fig. 12, reveals that the simple rescaling of the LO p T,H distribution by a constant K factor reproduces the NLO distribution only within ±10% in the relevant transverse-momentum range, if the transverse mass is chosen as the scale. The increase of the ratio dσ NLO /dσ LO with increasing p T,H can be understood in terms of the scale dependences of dσ NLO and dσ LO : with µ 2 = p 2 T,H + M 2 H rising with increasing p T,H , dσ LO decreases with increasing µ faster than dσ NLO (see Fig. 8). We note that choosing a fixed scale for the LO p T,H distribution provides a better description of the NLO shape.
For the transverse-momentum distribution of the top quarks, the increase of the ratio dσ NLO /dσ LO with increasing p T,H due to the scale dependence is balanced by gluon radiation off top quarks at NLO, which reduces the energies and momenta of the top particles and thereby depletes the region of large transverse momenta. As a net result, the shape of the transverse-momentum distribution of the top quarks is barely affected by the NLO corrections in the relevant range p T,t/t < ∼ 300 GeV, if the transverse top-quark mass is adopted for the renormalization and factorization scales (see Fig. 14).
The reduction in energies and momenta of the massive final-state particles due to gluon radiation at NLO is also expected to enhance the central rapidity region. However, as evident from Figs. 13 and 15, the effect on the (normalized) rapidity distributions of the massive final-state particles is small and becomes noticeable at the edges of phase space only.
Cum grano salis, the transverse-momentum and rapidity distributions of the Higgs boson and the top quarks for pp → ttH at the LHC, Figs. 16-19, show the same qualitative behaviour as observed for the Tevatron. Adopting the Higgs transverse mass as the scale, the difference between the NLO and LO Higgs-boson transverse-momentum distributions at the LHC is naturally more pronounced, with deviations of up to 30% at large p T,H . As for the Tevatron, choosing a fixed scale for the LO p T,H distribution provides a better description of the NLO shape.

Summary
In this report we have presented the theoretical analysis of the Standard-Model processes pp/pp → ttH + X in next-to-leading order QCD at the Tevatron and the LHC. We have focussed on two aspects.
(i) The technical basis for the prediction of the cross sections has been elaborated in detail. Particular emphasis has been put on evaluating pentagon diagrams, which involves the isolation of the singularities, the reduction to lower n-point functions, and the stable integration over the phase space. As a second important technical element we have applied the dipole subtraction formalism to treat massive final-state particles. In fact, the processes reported on here, represent the first complex examples that have been treated in this formalism in the case of massive quarks.
(ii) In expanding on results presented in an earlier letter, Ref. [ 10], we have analyzed not only the total cross sections for Higgs bremsstrahlung off top quarks at the Tevatron and the LHC, but also final-state rapidity and transverse-momentum distributions for the Higgs boson and the top quarks. Moreover, we have studied the uncertainty in the nextto-leading order cross section prediction due to different parametrizations of the parton densities.
When including higher-order QCD corrections, the K factor at the Tevatron has a value slightly below unity. It varies between ∼ 0.8 and ∼ 1.0 for renormalization and factorization scales between µ = µ 0 and 2µ 0 , almost independent of the Higgs-boson mass in the intermediate range. Here 2µ 0 = 2m t + M H denotes the threshold CM energy of the partonic subprocesses. Similarly, the K factor varies between ∼ 1.2 and ∼ 1.4 for the same scales at the LHC. Most importantly, the NLO predictions for the total cross sections and for the distributions in rapidity and transverse momentum of the Higgs boson and top quarks are stable when the renormalization and factorization scales are varied, in contrast to the Born approximation. The improved cross sections can therefore serve as a solid base for experimental analyses at the Tevatron and the LHC.       Figure 14: Normalized transverse-momentum distribution [ GeV −1 ] of the top and antitop quarks, (dσ/dp T,t + dσ/dp T,t )/2, for pp → ttH + X at the Tevatron in LO and NLO approximation, with the renormalization and factorization scales set to µ 2 = p 2 T,t/t + m 2 t .     Figure 18: Normalized transverse-momentum distribution [ GeV −1 ] of the top and antitop quarks, (dσ/dp T,t + dσ/dp T,t )/2, for pp → ttH + X at the LHC in LO and NLO approximation, with the renormalization and factorization scales set to µ 2 = p 2 T,t/t + m 2 t .   which agrees with the result in the literature [ 29]. So, by subtracting and adding welldefined singular integrals we are able to switch between arbitrary regularization schemes. For instance, to switch to an off-shell regularization scheme, with some of the external particles put slightly off-shell, it is sufficient to know the relevant 3-point integrals in that scheme.

B Colour operators and correlations
In order to make the results in this report as compact as possible, colour indices have been suppressed in the expressions by making use of colour operators. The definitions of these colour operators are given in this appendix (see also Refs. [ 23,24,36]).
Consider the matrix element for a process involving n coloured external particles, M a 1 ... an , where a 1 , . . . , a n are the colour indices of the coloured particles. These colour indices range from 1 to 8 for gluons [adjoint representation] and from 1 to 3 for (anti)quarks [fundamental representation]. Emission of a gluon with colour index c from the external particle j is represented by a specific colour operator (charge) T j acting on the matrix element: (M ⊗ T j ) c; a 1 ... an = b 1 , ... , bn δ a 1 b 1 · · · (T c j ) a j b j · · · δ anbn M b 1 ... bn = b j (T c j ) a j b j M a 1 ... b j ... an , (B.1) with (T c j ) ab = − if cab for incoming/outgoing gluons, (T c j ) ab = λ c ab /2 for outgoing quarks or incoming antiquarks, and (T c j ) ab = −λ c ba /2 for outgoing antiquarks or incoming quarks.

C The Altarelli-Parisi splitting functions at NLO
The so-called Altarelli-Parisi splitting functions P ab (x) [ 50] are a measure of the probability of finding the parton b with fraction x of longitudinal momentum inside the parent parton a. 6 At NLO these splitting functions can be written in terms of a regularized splitting function and an IR-sensitive part according to P ab (x) = P ab reg (x) + δ ab 2 T 2 The anomalous dimensions γ a are given by using the customary normalization T R = 1/2. The '+' distribution (. . .) + is defined in the usual way: If the Altarelli-Parisi splitting functions were calculated in D = 4 − 2ǫ dimensions by means of dimensional regularization, a term − ǫP ′ ab (x) would have to be added to Eq. (C.1). Such a term enters the final NLO expressions through the interference with the 1/ǫ collinear poles. The associated splitting functionsP ′ ab (x) amount tô Among other things, these splitting functions take into account the difference between the number of spin degrees of freedom of the quarks (2) and gluons (D − 2) in D dimensions.