One-loop QCD contributions to differential cross-sections for Higgs production at N$^3$LO

We present one-loop contributions to the fully differential Higgs boson gluon-fusion cross-section for Higgs production via gluon fusion. Our results constitute a necessary ingredient of a complete N$^3$LO determination of the cross-section. We perform our computation using a subtraction method for the treatment of soft and collinear singularities. We identify the infrared divergent parts in terms of universal splitting and eikonal functions, and demonstrate how phase-space integrations yield poles (up to $1/\epsilon^6$ ) in the dimensional regulator $\epsilon=(4-d)/2$. We compute the coefficients of the $\epsilon$ expansion, including the finite part numerically. As a demonstration of our numerical implementation, we present the corrections at N$^3$LO due to one-loop amplitudes in the rapidity and transverse momentum of the Higgs boson.


Introduction
Particle physics has entered an era of precision phenomenology which, at its core, aims to probe the interactions of the Higgs boson and other known or undiscovered particles. With the second run of the LHC, the statistical accuracy of the experimental measurements will increase significantly, allowing a precise determination of a variety of differential cross-sections and kinematic distributions. Precise theoretical predictions for fully differential cross-sections are highly desired. Their comparison to the measurements will offer valuable tests of the Standard Model and will set constraints to physics beyond the Standard Model.
Recently, the inclusive Higgs boson cross-section was computed through next-tonext-to-next-to-leading order (N 3 LO) in the strong coupling perturbative expansion [1,2]. Important achievements have also been accomplished towards differential Higgs cross-section at N 3 LO. The N 3 LO gluon-fusion Higgs production cross-section with a jet-veto has been obtained [3] by combining the fully differential cross-section for pp → H + 1 jet at NNLO [4][5][6][7][8] with the N 3 LO inclusive cross-section [1,2]. Other differential Higgs cross-sections have been computed by means of a threshold expansion [9][10][11] and the q T subtraction formalism which is being extended to N 3 LO [12].
We envisage a calculation of the fully differential Higgs cross-section with a direct subtraction of infrared divergences from the phase-space integrations over the partonic radiation associated with the Higgs boson production. At N 3 LO, one must add the fully differential partonic cross-sections for radiative processes of the type: parton + parton → Higgs + n partons, n ≤ 3. (1.1) The complete set of phase-space integrations in the above has been achieved inclusively, as it was required for the determination of the total Higgs cross-section. For a fully differential Higgs cross-section, the integrals over the phase-space must be performed as functionals of a generic infrared-safe measurement function. This can be achieved numerically, after the subtraction and cancelation of soft/collinear divergences. Processes with one parton in the final state (n = 1), represent the simplest nontrivial case of phase-space integrations. At N 3 LO, they receive contributions from the "real-virtual-virtual (RVV)" interference of two-loop amplitudes and tree-level amplitudes as well as the square of one-loop amplitudes (real-virtual-squared (RV) 2 ). The RVV contributions to the fully differential Higgs boson have been studied in Ref. [13].
In this article, we consider the (RV) 2 contributions, making a modest step towards a complete determination of the fully differential cross-section at N 3 LO. In particular, we revisit the soft and collinear singular limits of the one-loop amplitudes in terms of universal splitting amplitudes and soft-currents at one-loop. Then we isolate the -2 -singular terms of the partonic cross-sections with the aid of appropriate counterterms. This allows us to compute the coefficients of the expansion in the dimensional regulator ε = 2 − d 2 numerically for arbitrary measurement functions.
As explicit examples, we present these coefficients differentially, in bins of the Higgs transverse-momentum or its rapidity.

Setup
We consider processes where i, j, k denote quark, antiquark or gluon partonic flavors, p 1,2 are the momenta of the incoming partons, p h is the momentum of the Higgs boson and p 3 is the momentum of the radiated parton. The Mandelstam variables are: with s 12 + s 23 + s 13 = m 2 h . We parameterise the final-state momenta in terms of dimensionless positive variables z, λ ≤ 1 as in: We evaluate perturbatevely the amplitudes for the processes (and the ones related to the above by crossing symmetry and/or charge conjugation) in the Standard Model and in the limit of a very heavy top-quark. The leading contribution to this limit in the strong-coupling sector is described by the Lagrangian density: where L QCD is the QCD Lagrangian density (with n f = 5 massless quark flavours and N c number of colours), h is the Higgs boson field, G µν the gluonic field-strength tensor -3 -and C 1 is the Wilson coefficient [14][15][16] that arises from matching the effective theory to the full Standard Model. We renormalise the bare strong coupling constant α s ≡ g 2 s 4π and the Wilson coefficient in the MS renormalisation scheme: where the multiplicative factors Z α (µ 2 ) and Z C (µ 2 ) are given by: The renormalised Wilson coefficient is given by: We compute the required one-loop amplitudes as well as their soft/collinear limits in conventional dimensional regularisation (CDR). The "form-factors" of the amplitudes as computed in CDR suffice to determine fully the amplitudes for the scattering of partons of definite helicity. The universal collinear and soft limits of helicity amplitudes are known in the literature [17][18][19][20][21] and we verify that our results agree with them.

Tree and one-loop amplitudes
In this section, we present the tree and one-loop amplitudes which are required for the gluon-fusion Higgs production cross-section at N 3 LO in perturbative QCD.

The gg → h amplitude
Let us start first with the gluon-gluon scattering process with p h = p 1 + p 2 . For physical external polarisations, ε(p) · p = 0, we can write the amplitude as The coefficient A h admits a perturbative expansion in the bare strong coupling constant α s , In the following, we will be concerned only with the first two terms in the expansion, which read

The gg → gh amplitude
We now the consider the process with p h = p 1 + p 2 + p 3 . If we choose polarisation vectors ε µ i ≡ ε µ (p i ) for the external gluons which satisfy we can cast the amplitude in the form: Squaring and summing over polarisations and colour, we find pols cols We can further relate the form factors A 1 , A 2a , A 2b , A 2c to helicity amplitudes [22]: (3.12) , (3.14) where the coefficients α i 's are related to the amplitude coefficients A i 's via Notice that all of the above relations are valid at any order in the perturbative expansion in the strong coupling constant. We now expand the form-factors perturbatively: The leading order amplitude coefficients A i are where the final-state momenta are given by Eq. (2.3). At one-loop, the amplitude coefficients A 2c are linear combinations of the bubble and box integrals, which are defined as In the above, (q 1 + q 2 ) 2 = s, (q 2 + q 3 ) 2 = t, (q 1 + q 3 ) 2 = u, q 2 1 = q 2 2 = q 2 3 = 0 and s + t + u = m 2 . The hypergeometric function can be expanded in ε in terms of polylogarithms, The arguments of the master integrals which appear in the amplitudes are (3.26) The expressions for the amplitudes A i written in terms of these Master Integrals are given in the appendix A. In order to simplify the notation, we set s 12 = 1 in these expressions and the rest of this work. We also denote for any integer a, bearing in mind that the quantity in the parenthesis has a small negative imaginary part. The mass dimensions of any quantity can be recovered easily with dimensional analysis.

The qg → qh and qq → gh amplitudes
The amplitudes for the q(p 1 )+g(p 2 ) → q(p 3 )+h(p h ) and the q(p 1 )+q(p 2 ) → g(p 3 )+h(p h ) processes are related by crossing symmetry. We will therefore present here only the amplitude for the former. For physical gluon polarisations, it takes the form Squaring and summing over spins, polarisations and colour, we find spin cols Each of the amplitude coefficients in Eq. (3.28) can be expanded as a power series in the strong coupling constant At leading order 3 amplitude coefficients are presented in the appendix A.

Infrared divergences of one-loop amplitudes
The one-loop amplitudes of the previous section are divergent in d = 4 dimensions, due to singularities when the momenta of two adjacent massless particles become collinear, or when a massless particle is soft. These divergences cancel in a complete hadronic cross-section computation. In the singular limits, the amplitudes exhibit universal factorisation properties. We will exploit them in order to isolate the divergent parts and to facilitate the integration of the one-loop contributions to the partonic crosssections, which we are computing here, into a future complete hadronic cross-section computation.
In the following subsections, we will recall the factorization of the amplitudes in the limits where two external partons become collinear or an external parton becomes soft.

Collinear limits
In the limit where two external-particles become colliner, colour ordered amplitudes factorise in a universal way [13,17,18,22,[24][25][26][27][28]. In this section we will compute explicitly the collinear limits of the gg → gh and gq → gh amplitudes and cast them in terms of universal functions related to the tree and one-loop splitting amplitudes.

Collinear limits for the gg → gh amplitude
Let us first consider the limit for p 3 ||p 2 becoming collinear (the other limit p 3 ||p 1 can be derived from this by symmetry At leading order, we find that At one-loop, we find the following collinear limits for the amplitude coefficients, where the universal functionsr 1 andr 2 are given bỹ and the functions f 1 and f 2 are defined as Our results are in agreement with Ref. [17,18].

Collinear limits for the qg → qh and qq → gh amplitudes
The collinear limit p 3 ||p 1 is common to both one-loop amplitude coefficients A 2 and A (1) 3 for the qg → qh amplitude. Specifically, we find where The above results are in agreement with Ref [17,18]. The other limit for p 2 ||p 3 (λ → 0) is not singular. Similarly, the one-loop amplitudes for the process qq → gh are not singular in the above collinear limits.

Soft limit
We now turn our attention to the factorisation of the amplitudes in their soft limits. At tree-level, the soft limit,z → 0, of the gg → gh amplitude is 19) and, at one-loop, The above results are in agreement with Ref. [18,21]. The one-loop amplitudes for both qg → qh and qq → gh processes are not singular in the soft limit.

Hadronic cross-section and subtraction of infrared divergences
We now consider the hadronic production of a Higgs boson in association with a parton i in the final state, The proton momenta P i in the hadronic centre of mass frame are given by with √ S being the collider centre of mass energy. The hadronic cross-section is given by the integral x 1 and x 2 are the Bjorken fractions and f i (x i , µ F ) are the renormalised parton distribution functions. The sum runs over all the initial-state parton pairs. For the computation of the fully differential cross-section at N 3 LO, we need, among other contributions, to integrate the matrix-elements of Eq. (3.10) and Eq. (3.29) over the phase-space of the final-state partons, weighted with an infrared-safe measurement function J O (p h ). The matrix-elements are independent of the azimuthal angle φ. The momentum of the Higgs boson can be conveniently parametrised in terms of its rapidity Y h and transverse momentum p T , defined in the plane transverse to the beam axis as In terms the partonic variables of Eq. (2.3), these variables can be rewritten as The partonic cross-sectionσ ij , differential in the variable λ, for the different initial state parton contributions is given by where M ij are the corresponding matrix elements, while the phase space for the production of the Higgs plus one final state parton is given by In order to be able to compute a finite differential cross-section, we need to subtract the collinear and soft singularities that arise in the corresponding limiting kinematics. More specifically, we find that all the partonic cross-sectionsσ ij have only single poles in the variables z and λ of the form where the coefficients a i are integer numbers. Therefore, we can regulate the divergencies by subtracting the corresponding limiting contribution at the integrand level and adding back the integrated counterterm, as shown schematically in the following example for a general singular variable x, For the gg initial states, the collinear limiting behaviour for the matrix elements squared can be cast in the following form in terms of the universal functions r 1 and r 2 , pols cols and similarly for pols cols |M (1) gg | 2 (z, λ → 1, m 2 h ), with the exchange λ ↔λ. Here N gg is the initial averaging factor over the spins and colours, and we have already normalised the coupling constant with the factor e εγ E (4π) ε . The soft subtraction term, instead, is given by pols cols For the qg initial states, the matrix elements squared in the collinear limit λ → 1 take the form pols cols with the initial averaging factor over the spins and colours, We recall that the partonic cross-section for the qq initial states does not present any collinear nor soft singularity, and hence no subtraction is needed for this contribution.
With the aid of Eq. (5.11) and the above infrared limits, we can expand the partonic cross-sections in the dimensional regulator ε. The resulting expressions are lengthy but straightforward to derive and we refrain from presenting them here. Our results are in agreement with an independent calculation performed for the purposes of Ref. [9].
-13 -We implement the results for the integrated counterterms and the finite part of the unrenormalized cross-section in a private numerical code and compute distributions for the Higgs rapidity and transverse momentum. Note that in the subtraction procedure we adopt to remove the soft and collinear divergences, we need to add back integrated counterterms with explicit poles up to ε −6 .
We consider a centre of mass energy of 13 TeV. We use the NNLO PDF4LHC15 set [29] for the parton distribution functions, as available from LHAPDF [30] and evolve the strong coupling constant α s to NNLO. In addition, we set the Higgs boson mass to the value of m h = 125 GeV and take the same value of µ = m h for the renormalisation and factorisation scales.
The distributions we obtain, albeit non-physical as they correspond only to a part of the complete N 3 LO calculation, serve as a validation of our results and of the subtraction method employed to remove the infrared divergences. In particular, this procedure provides a stable and reliable numerical implementation.
In figure (1) we show the distributions for the transverse momentum of the Higgs boson p T for all the six ε poles and the finite part of the cross-section. In the ones corresponding to the deepest poles ε −6 and ε −5 we recognise the expected behaviour, as these contributions are proportional to the born cross-section. For what concerns the other poles, as there is not a direct interpretation in terms of specific kinematic configurations, the behaviour is anyway unphysical, with rather big cancellations occurring between the first bins. Nevertheless, we have checked that the sum of all the bins returns the total cross-section for the corresponding pole in the dimensional regulator ε.
In figure (2) we show the rapidity distribution for all the six ε poles and the finite part of the cross-section for positive values of the rapidity Y h , since it is symmetric in We have presented the computation of the real-virtual squared contributions to the fully differential cross-section for Higgs production via gluon fusion. This is a part of the complete N 3 LO fully differential cross-section. One of the main results of our computation is the analytic expression for all the soft and collinear counterterms needed for the subtraction of the infrared divergencies. In particular, we have been able to express these counterterms in terms of universal functions related to the QCD soft currents and splitting amplitudes at one loop. Therefore, the results obtained here can be directly used in other processes containing a colourless final state particle and allow for an easier and general numerical implementation of differential cross-sections for such processes.
In order to validate our results and the subtraction procedure employed, we implemented the contributions corresponding to the six poles in the dimensional regulator ε of the integrated counterterms and the finite part of the cross-section in a numerical code and obtained distributions for the transverse momentum and rapidity of the Higgs boson. These results, although non-physical, reproduced the expected behaviour for the deepest poles ε −6 and ε −5 corresponding to born kinematic configurations, while providing a stable numerical evaluation for all the higher terms through the finite part of the ε expansion.