Two-loop splitting amplitudes and the single-real contribution to inclusive Higgs production at N3LO

The factorisation of QCD matrix elements in the limit of two external partons becoming collinear is described by process-independent splitting amplitudes, which can be expanded systematically in perturbation theory. Working in conventional dimensional regularisation, we compute the two-loop splitting amplitudes for all simple collinear splitting processes, including subleading terms in the regularisation parameter. Our results are then applied to derive an analytical expression for the two-loop single-real contribution to inclusive Higgs boson production in gluon fusion to fourth order (N3LO) in perturbative QCD.


Introduction
Scattering amplitudes in massless QCD diverge if one or more of the external momenta become soft or if two or more external momenta become collinear. In these infrared limits, one observes process-independent factorisation properties: the divergent behaviour of the amplitude is described by universal unresolved factors, multiplying a reduced amplitude with lower partonic multiplicity. These unresolved factors are expanded in a two-fold manner in perturbative QCD: in the number of unresolved partons and in the number of virtual loops correcting the factors at fixed multiplicity. The eikonal factors describing the simple soft behaviour and the collinear splitting amplitudes were established already long ago [1,2]. They are sufficient to understand the infrared singular structure of real radiation contributions at next-to-leading order (NLO) in QCD, and have been instrumental in devising systematic subtraction schemes for NLO calculations [3,4,5]. They also form the starting point for a systematic expansion around unresolved limits to all orders in perturbation theory [6,7].
At next-to-next-to-leading order (NNLO), tree-level contributions with up to two unresolved particles contribute, again displaying universal factorization properties in collinear and soft limits [8,9]. At the same order, one-loop corrections to simple collinear [27] and simple soft [10] emissions need to be accounted for. The behaviour in all these limits is understood in detail, and served as a construction principle for subtraction methods at NNLO [11,12,13]. In an actual calculation of higher order corrections to a collider observable, the universal soft and collinear factors need to be integrated over their appropriate phase spaces, including a regulator for the infrared divergences. Consequently, subleading terms in the regulator will yield finite contributions to the final result: simple soft configurations demand two subleading orders, and simple collinear configurations one subleading order.
With advances in multi-loop calculations, calculations of next-to-next-to-next-to-leading order (N 3 LO) corrections to benchmark observables are now becoming feasible. Multiple unresolved limits at zero and one loop order were understood already some time ago [14,15,16]. Two-loop corrections to simple soft and simple collinear limits were derived only to finite order in the regulator [17,18] by taking the appropriate limits of the two-loop three-parton decay matrix elements [19,20]. With the two-loop soft gluon current recently derived to all orders in the regulator [21,22], N 3 LO calculations at the two-loop single real level are now missing only the two-loop corrections to simple collinear limits, described by the two-loop splitting amplitudes. It is the aim of this paper to derive these, based on a new derivation of the two-loop matrix elements for three-parton decays, expanded to higher orders in the regulator.
The most important application of our newly derived results are the N 3 LO corrections to inclusive Higgs production in gluon fusion. Following the derivation of the three-loop virtual corrections [23], of the required renormalization and factorization counterterm contributions [24], and of expansions of single and multiple real radiation contributions around the soft limit [25,21,22], the threshold approximation of this coefficient function was computed earlier this year [26]. We use our results for the two-loop splitting amplitudes to -1 -derive a closed analytic expression for the two-loop single-real contribution to the inclusive N 3 LO coefficient function for Higgs production. Depending on the availability of results for the other channels (one-loop double-real and triple-real), our expression can be used to derive further terms in the threshold expansion, or to obtain an exact expression for the full coefficient function.
Our paper is structured as follows. In Section 2, we establish the notation and discuss the behaviour of multi-loop amplitudes in single collinear limits. Section 3 describes the extraction of the two-loop splitting amplitudes from the calculation of two-loop matrix elements in three-particle decay kinematics. The analytical expressions for the splitting amplitudes are quite lengthy (especially in the term subleading in the dimensional regulator), and are collected in Appendix A. These results are then applied in Section 4 to compute the two-loop single-real contribution to Higgs production at N 3 LO. The two-loop matrix elements relevant to this process were previously known only for external four-dimensional helicity, and are re-derived in conventional dimensional regularization in Appendix B, their contribution to the inclusive coefficient function is documented in Appendix C. We conclude with an outlook in Section 5.

Collinear limits of multi-loop QCD amplitudes
We consider an ℓ-loop QCD amplitude with a certain number of coloured partons in the final state. Our example will be the decay of a heavy colourless state X into n partons with momenta p i , i = 1, . . . , n, and all other cases can be obtained by crossing symmetry. In the limit where a pair of partons become collinear, say i and j, then the corresponding Mandelstam invariant vanishes, s ij = 2p i · p j → 0. The amplitude diverges in the limit, and the divergence is characterised by a simple pole at s ij = 0. The residue at the pole can be described by the factorisation formula where |M (ℓ) (. . . , p i , p j , . . .) is the ℓ-loop amplitude and α 0 the bare strong coupling constant and P = p i + p j the four-momentum of the parent parton. We work in dimensional regularisation in D = 4 − 2ǫ dimensions, and with γ E = −Γ ′ (1) the Euler-Mascheroni constant. The '≃' sign denotes 'equality up to terms that are power-suppressed in the limit'. The function P (ℓ) ij (z) is the ℓ-loop splitting amplitude, and depends on the colour and spin of the particles involved in the splitting, as well as on the momentum fraction z carried by particle i, p i = z P . Note that, in contrast to the well-known Altarelli-Parisi splitting functions, the splitting amplitudes considered -2 -in this paper only describe the collinear behaviour of a system of two real partons with momenta p i and p j , i.e., they do not include purely virtual 1PI corrections and multiple real radiation corrections. The (spin-averaged) tree-level splitting amplitudes are given by [1,9] and P In the remainder of this section we summarise the main properties of splitting amplitudes. It is obvious that Bose symmetry implies The limits where z → 0 or z → 1 correspond to the limits where one of the particles involved in the splitting becomes soft. The splitting amplitude for g → qq does obviously not have any soft singularity, and so the residue of P (ℓ) qq (z) at z = 0 and z = 1 vanishes. The splitting amplitude P (ℓ) gq (z) only develops a soft singularity if the gluon becomes soft, i.e., z → 0, and the residue at the pole is given by the QCD soft current, known up to two loops in perturbative QCD [2,10,18,21,22], S denotes the ℓ-loop QCD soft current (in the normalisation of ref. [21]). Similarly, the gluon splitting amplitude becomes singular if any of the two final-state gluons becomes soft, (2.6) The pole structure of P (ℓ) ij (z) in dimensional regularisation is completely fixed up to two loops. At one-loop order, we have with where β l is the l-loop QCD β function. At two-loop order we have [7,18] P (2) ij (z) = I (2) where the H and with N f the number of flavours. We recall that the above expressions refer to the splitting amplitudes describing the collinear behaviour of the unrenormalized one-loop and two-loop matrix elements, as in [18]. Formulae for the factorization of the pole structure can also be derived after renormalization [7]. Splitting amplitudes have been computed in the literature up to two-loops. The treelevel and one-loop splitting amplitudes are known to all orders in the dimensional regulator ǫ (in a variety of schemes) [1,27]. The two-loop splitting amplitudes have so far only been computed up to finite terms [18,17]. Our goal is to compute the two-loop splitting amplitudes to O(ǫ) in conventional dimensional regularisation (CDR). This result is needed in order to construct the collinear counterterms for the double-virtual-real contributions to the inclusive N 3 LO cross section for the production of a heavy colourless state.

Computation of the splitting amplitudes
In this section we describe the computation of the two-loop splitting amplitudes to O(ǫ) in CDR. The splitting amplitudes can be extracted most conveniently from the two-loop -4 -amplitude for a heavy colourless state decaying into three massless partons. In particular, we have H→gg , H→gg , with z = s 23 /Q 2 , where Q 2 is the virtuality off the initial state. The hard matrix elements are given by where λ 0 and e q denote the bare coupling of the Higgs boson to the gluons (in the large top mass limit) and the bare electromagnetic coupling of the quarks, and F (ℓ) g and F (ℓ) q denote the ℓ-loop gluon and quark form factors, known up to three loops in QCD [23]. As a consequence, if we know the the matrix elements in the left-hand side in CDR through O(ǫ) up to two loops, we can immediately extract the corresponding splitting amplitudes through O(ǫ). The two-loop matrix element for γ * → qq g in CDR was presented to all orders in ǫ in terms of the master integrals in ref. [19]. In ref. [20] the two-loop Ddimensional tensor coefficients as well as the two-loop helicity coefficients (in the 't Hooft Veltman) scheme for H to three partons were presented. In Appendix B we show how to construct the matrix elements in CDR from the tensor coefficients of ref. [20]. The two-loop matrix elements are given, to all orders in the dimensional regulator, in terms of the master integrals of ref. [28,29], where the complete set of master integrals was evaluated in terms of harmonic polylogarithms (HPLs) and two-dimensional harmonic polylogarithms (2dHPLs) [30,28] up to transcendental weight four. This is sufficient to compute the corresponding matrix elements up to finite terms (some of the master integrals are known to all orders in ǫ, see, e.g., ref. [31]). Since our goal is to compute the splitting amplitudes to O(ǫ), we need to compute all the masters integrals to one order higher in the ǫ expansion, i.e., up to terms of weight five. Moreover, it turns out to be convenient to have the master integrals in a form where all the divergent logarithms in the collinear limit are resummed in ǫ, because the ℓ-loop splitting amplitude can then immediately be read off from the coefficient of s −ℓǫ ij in the residue of the matrix element at s ij = 0 (and dividing out the hard matrix element, cf. eq. (3.1)). We have therefore explicitly computed all the master integrals of ref. [28,29] up to transcendental weight five using the method of differential equations [31,32]. We choose a canonical basis of master integrals where all the master integrals are of uniform transcendental weight [33]. In this basis, every master -5 -integral is a function of three variables x ij (after dividing out an overall scale), defined by and constrained by x 12 + x 13 + x 23 = 1. Every uniformly transcendental master integral then admits a representation of the form 1 with u i,m,n,k (x 13 , x 23 ) a linear combination of 2dHPLs with rational coefficients that are analytic in a neighbourhood of (x 13 , x 23 ) = (0, 0). We have explicitly computed the functions u i,m,n,k (x 13 , x 23 ) up to weight k = 5, which is sufficient to compute the two-loop matrix elements up to order ǫ. The initial conditions for the differential equations were obtained from the leading term in the expansion of the master integrals around the soft limit [21]. We checked that our results agree with the first few subleading terms in the soft expansion of the master integrals. Once the master integrals are known up to terms of transcendental weight five, we can immediately compute the residue at x 13 = 0 by expanding in this variable (while keeping the logarithmic divergencies resummed in the form x −nǫ 13 ). The variable x 13 only enters the matrix elements through rational functions, as well as through 2dHPLs of the form G( a; x 13 ), where a is independent of x 13 , and the last entry in a is non zero. 2dHPLs of this type can easily be expanded into a power series close to x 13 = 0 by using the series representation of multiple polylogarithms, with 0 m = (0, . . . , 0 m times ) and the Z m 2 ,...,m k denote Z sums [34]. The splitting amplitudes up to two-loops and up to O(ǫ) can then immediately be read off by comparing the expression to eq. (3.1). At one loop we find where C i+j is the Casimir of SU (N ) in the representation of the parent parton, and A (1) (z, ǫ) and B (1) (z, ǫ) are functions that are analytic close to z = 0. They are independent of the identities of the partons involved in the splitting and can be written as a combination of harmonic polylogarithms of uniform weight, where we used a shorthand for harmonic polylogarithms, H i,...,j ≡ H i,...,j (z). The function R (1) ij (z, ǫ) are rational functions of z order by order in the ǫ expansion and depend on the identity of the partons involved in the splitting, (3.8) We emphasise that the previous expressions are valid for a time-like splitting where the collinear pair is in the final state. The analytic continuation to a space-like splitting can easily be performed by the replacement z −kǫ → e −iπkǫ z −kǫ (3.9) in eq. (3.6). At two loops, we write where P (2,k) ij (z) are analytic in a neighbourhood of z = 0. Their analytical expressions are given in Appendix A.
We checked that the tree-level and one-loop are correctly reproduced from our two-loop computation. Moreover, we checked that our results have the correct symmetry properties and soft limits, eq. (2.4 -2.6), as well as the the correct pole structure in dimensional regularisation, eq. (2.10).

Two-loop single-real emission contributions to Higgs production at N 3 LO
As an application of the splitting amplitudes up to two loops computed in the previous section, we present in this section the integration over phase space of the matrix elements -7 -for the scattering processes These processes are important because they contribute to the inclusive Higgs production cross section at N 3 LO in perturbative QCD. Indeed, the final-state parton may become collinear to any of the two incoming partons, and the matrix elements develop singularities which are described by the splitting amplitudes we have just computed. The knowledge of the spitting amplitudes therefore enables us to write down explicit counterterms that we can subtract from the matrix element in order to render the phase-space integrations finite. Moreover, the counterterms will turn out to be trivial to integrate over the unresolved phase space, giving rise to a pole in dimensional regularisation. This pole, in turn, is at the origin of the appearance of the O(ǫ) part of the splitting amplitudes computed in the previous section in the finite terms which enter the cross section. We start by giving some generic considerations about the phase-space parametrisation we will use in the following to perform the integration, and we discuss the different channels separately in subsequent subsections.
We denote the momenta of the two incoming partons by q 1 and q 2 , and the momentum of the final-state parton is denoted by q 3 . Furthermore, we define Note that s, t, u > 0 and s = m 2 H + t + u. The phase space integrals we want to compute are N X 2s where X denotes the channel under consideration, and N X denotes the averaging factor over the initial-state colours and spins, , The D-dimensional phase-space measure is given by We parametrise the phase space by In this parametrization the phase-space measure becomes where Θ denotes the step function. In the following we discuss the integration over phase space of the three different processes in eq. (4.1), and we separately discuss the processes with different initial states.

The qq initial state
The two-loop matrix element for qq → H g can be written as, with s = µ 2 = 1 and M (2) qqg is a Laurent series in ǫ whose coefficients are rational functions of the x ij multiplied by 2dHPLs. The expression of M (2) qqg in terms of the tensor coefficients of ref. [20] can be found in Appendix B. Note that M (2) qqg develops an imaginary part, because the master integrals have branch points at x t = 0 and x u = 0, The functions u i,m,n,k (x t , x u ) are real, because they are analytic in a neighbourhood of (x t , x u ) = (0, 0). In the previous section, we argued that the amplitude is only singular if the quark pair becomes collinear. This implies that M qq→Hg is finite everywhere inside the phase space, and so we can simply expand in ǫ and perform the integral over λ order by order in ǫ. The integrand contains 2dHPLs of the form G( a(x u ); −x t ), where x t and x u are rational functions of λ andz (cf. eq. (4.6)). Interpreting the 2dHPLs as special instances of multiple polylogarithms [35], and using the techniques of ref. [36], we can easily convert all the 2dHPLs of the form G( a(x u ); −x t ) into multiple polylogarithms of the form G( b(z); λ) and multiple polylogarithms inz. As an example, the following relation holds (4.11) Similar identities can easily be derived in all other cases, and the phase-space integration over λ can easily be performed using the the definition of the multiple polylogarithms, After a final massaging, we find that all the phase space integrals can be expressed in terms of HPLs in 1 −z, in agreement with the expectations.

The qg initial state
The two-loop matrix element for q g → H q can be written as, with s = µ 2 = 1, (4.13) is identical to the corresponding quantity in the qq → H g case, up to analytic continuation. Note that there is an additional minus sign, coming from the crossing of a single fermion to the initial state. The analytic continuation of the function M can be performed in a similar way as in the qq → H g case.
The integrand becomes singular when the initial and final-state quarks are collinear, and the singularity is controlled by the splitting amplitudes computed in Section 3. More precisely, after inserting the parametrisation (4.2), the integrand has a simple pole as λ → 0, and from eq. (3.1), we have Note that the splitting amplitude is evaluated at λ = 0, so that we only subtract the residue of the matrix element at λ = 0. The phase-space integral then becomes (4.16) The integral in the second line is convergent, and so we can expand in ǫ under the integration sign and integrate order by order in ǫ. At this point, however, we need to face a technical difficulty: performing the integrations using eq. (4.12) requires the integrations to be performed term by term. In general, these integrals will be divergent, because only the sum in the second line of eq. (4.16) is finite. We therefore introduce a cut-off δ for the lower integration limit 2 and perform all the integrations term by term. As a result, all the integrals are convergent, but they explicitly depend on the cut-off. We observe that all divergences cancel when we expand the result around δ = 0, in agreement with the expectation that the integral in the second line of eq. (4.16) is convergent.

The gg initial state
The two-loop matrix element for g g → H g can be written as The integrand does not develop poles for λ → 1, and so no cut-off is required for the upper integration limit.
-10 - The expression of M (2) ggg in terms of the tensor coefficients of ref. [20] is given in Appendix B. Its analytic continuation can be performed in the same way as in the previous cases. The phase-space integral diverges in the collinear limits where either t or u vanish. Moreover, the gluon amplitude also develops a singularity when the final state gluon becomes soft, corresponding to a simultaneous vanishing of t and u. In terms of our phase space parametrisation (4.6) this corresponds toz → 0, and in this limit the matrix element behaves like We define the soft-regularized amplitude, defined as the amplitude with the soft singularity atz = 0 subtracted, .
(4.20) The soft-regularised amplitude is free of soft singularities, but still contains non-overlapping collinear divergences. The residues at the collinear poles at λ = 0 and λ = 1 can be obtained in a way similar to the q g → H q case considered in the previous section. By Bose symmetry, the residues at the two poles are identical. In particular, the residue at λ = 0 is given by and The phase-space integral can now be written as The remaining integral is convergent and can be performed order by order in ǫ using eq. (4.12). Individual terms in the integrand can develop poles as λ approaches either 0 or 1. We therefore introduce a cut-off for each integration limit, and we observe that all the singularities cancel in the sum over all terms when the cut-offs are removed, confirming the expectation that the remaining integral in eq. (4.24) is convergent.

Conclusion
The single collinear factorization of scattering amplitudes in massless QCD is determined by so-called splitting amplitudes, which can be expanded perturbatively in the number of virtual loops [6]. We derived the two-loop corrections to all splitting amplitudes, finding agreement with previously known results [17,18] up to finite terms in the dimensional regulator, and deriving the subleading terms for the first time. When integrated over the unresolved single collinear phase space in actual N 3 LO calculations, these subleading terms produce finite contributions. We used our result to analytically compute the twoloop single-real contribution to the coefficient function of inclusive Higgs boson production in gluon fusion, improving upon the previously available soft approximation [21,22] of this function.
When combined with either expanded or full results for the other channels (one-loop double-real and triple-real), our expression will allow to derive further terms in the threshold expansion [26] or to complete an exact calculation of the full N 3 LO coefficient function for Higgs production in gluon fusion.

A. Splitting amplitudes at two loops
The two-loop splitting amplitudes can be expanded according to the powers z −kǫ , with coefficients that are analytic around z = 0, as defined in (3.10). We write these coefficients as: Note that these expressions are again valid for a time-like splitting, and the space-like case can again be obtained by the replacement (3.9). The functions U (2,k) ij (z) are of uniform weight, and are given by The functions V (2,k) ij (z) are not of maximal weight, and are given by      In this section we show how to obtain the two-loop matrix element in CDR from the Ddimensional tensor coefficients of ref. [20]. The amplitude for H → ggg can be written as where the gluon tensor is given by, where λ 0 and α 0 denote the bare coupling constants, and where γ E = −Γ ′ (1) denotes the Euler-Mascheroni constant. After tensor decomposition, we get where T ijk are the tensors given in eq. (3.7) of ref. [20], and the A The gauge dependent terms are proportional to the gluon momentum. The tensors T ijk ≡ T µνρ ijk ǫ µ 1 ǫ ν 2 ǫ ρ 3 of ref. [20] are transverse, i.e., they vanish whenever they are contracted with an external gluon momentum, As a consequence, the gauge dependent terms will always drop out, and we can use the 'naive' polarisation sum pol.
We get (B.11) 3 We slightly changed the normalisation of the A (ℓ) ijk compared to ref. [20], in order to have the colour tensor and the overall coupling constant explicit.

B.2 The matrix element for H → qqg
The amplitude for H → qqg can be written as where as in the gluon case we have factored out explicitly the overall coupling and colour structure. We have The colour sum is trivial and the quake spin sum give 14) The tensors T i ≡ T iρ ǫ ρ 3 are transverse, T iρ p ρ 3 = 0, and so we can again use the 'naive' polarisation sum for the gluons. We get, qqg (x 12 , x 13 , x 23 ; ǫ) , where we used the fact that A where v ≃ 246 GeV is the vacuum expectation value of the Higgs field.

C.1 The gg initial state
The contribution of the two-loop amplitude for gg → Hg to the inclusive Higgs cross section at N 3 LO can be written aŝ The first term represents the soft-virtual term, which is entirely determined by the QCD soft current,σ is defined in eq. (4.19). The second term is regular as z → 1, and can be written asσ -37 - -39 -

C.2 The gq initial state
The contribution of the two-loop amplitude for gq → Hq to the inclusive Higgs cross section at N 3 LO can be written asσ Note that, unlike for the pure-gluon initial state, the result is regular in the limit z → 1, and so we do not need to separate off the contribution from the soft limit. The coefficient σ (3) gq→Hq can be written aŝ

C.3 The qq initial state
The contribution of the two-loop amplitude for qq → Hg to the inclusive Higgs cross section at N 3 LO can be written asσ qq→Hg . (C.85) The result is regular in the limit z → 1, and so we do not need to separate off the contribution from the soft limit. The coefficientσ (3) qq→Hg can be written aŝ