Differential Higgs production at N3LO beyond threshold

We present several key steps towards the computation of differential Higgs boson cross sections at N$^3$LO in perturbative QCD. Specifically, we work in the framework of Higgs-differential cross sections that allows to compute precise predictions for realistic LHC observables. We demonstrate how to perform an expansion of the analytic N$^3$LO coefficient functions around the production threshold of the Higgs boson. Our framework allows us to compute to arbitrarily high order in the threshold expansion and we explicitly obtain the first two expansion coefficients in analytic form. Furthermore, we assess the phenomenological viability of threshold expansions for differential distributions. In addition, we report on an interesting obstacle for the computation of N$^3$LO corrections with LHAPDF parton distribution functions and our solution. We provide files containing the analytic expressions for the partonic cross sections together with the arXiv submission.


Introduction
The discovery of the Higgs boson in 2012 by the ATLAS [1] and CMS [2] experiments at the Large Hadron Collider (LHC) was a major scientific breakthrough. Confirming for the first time conclusively the presence of the Higgs field and explaining the origin of the masses of elementary particles, it renders the standard model mathematically self-consistent, allowing it to be used to formulate credible predictions at high energies. At the same time it is clear that the standard model needs to be extended to explain the puzzle of neutrino masses and cosmological observables as well as the nature of dark matter.
If these open questions are related to the origin of masses it is likely that the precise nature of the Higgs boson will differ quantitatively from what is expected in the standard model. Lacking direct observations of particles beyond the standard model, the Higgs boson is the most promising probe into possible ultra-violet completions of the standard model. As such, a precise study of the properties of the Higgs boson, such as mass, spin/parity, branching ratios and production rates, is a central part of the present and future precision Higgs physics program at the LHC that has been initiated since the discovery. The remarkable performance of the LHC, delivering unprecedented luminosity, as well as the unrelenting efforts of the experiments have led to vast improvements of Higgs analyses resulting in more and more precise extractions of the properties of the Higgs boson.
In order to truly exploit the potential of these excellent experimental results, it is imperative to confront them with equally precise theoretical predictions. This demand for precise predictions from theory has lead to a flurry of calculations in recent years at next-to-leading (NLO) and next-to-nextto-leading (NNLO) order in perturbative QCD. Higher order predictions are particularly important for Higgs phenomenology, in part due to the slow convergence of the perturbative expansion in the strong coupling constant for the dominant mode of Higgs hadroproduction via gluon fusion. The phenomenological importance of these corrections can be seen from the large size of the NLO corrections [3,4] which almost double the leading-order prediction [5]. The magnitude of these perturbative corrections indicate potentially significant contributions from even higher orders in the perturbative series, leading to a substantial uncertainty on the gluon-fusion cross section. Even with the inclusion of the NNLO [6][7][8] corrections to the gluon-fusion cross section the perturbative series seems to converge slowly keeping the perturbative uncertainty at a significant level.
Recently, the gluon-fusion cross section has been computed through next-to-next-to-next-toleading order (N 3 LO) in perturbative QCD [9][10][11] in the limit of an infinite top mass. This calculation has lead to a reduction of the perturbative uncertainty, making the theoretical predictions competitive with current experimental analyses. At this level of precision, effects that go beyond the leading approximation of an infinite top mass or a treatment in pure QCD, neglecting effects from quark masses or electroweak loops, become important. The state of the art predictions for Higgs hadroproduction have been combined in a consistent way in [12] and are also compiled in [13,14].
The Higgs cross section is measured in dedicated regions of phase space as required by the detector geometry and optimized by carefully designed experimental event selections. Within these acceptance regions, the experiments have excellent capabilities to measure a plethora of kinematic distributions for the Higgs boson and its decay products that can be used to characterize the properties of the Higgs boson. As such it is imperative for theoretical predictions to not only inclusively describe the total cross section but also provide high precision theoretical predictions for differential cross sections. Recently, the pp → H + 1 jet fully differential cross section has been computed at NNLO [15][16][17]. In combination with the inclusive N 3 LO cross section this enables the computation of the jet-vetoed Higgs cross section at N 3 LO [18]. Fully differential parton-level Monte-Carlo simulations at N 3 LO will enable the study of efficiencies of many other event selection criteria at the same accuracy in perturbation theory as the jet-veto efficiency.
One way to achieve this goal would be to attempt to generalize any of the methods available at NNLO (sector-decomposition [19][20][21][22], slicing [23][24][25][26], subtraction [27][28][29][30], reweighting [31] and other methods [32]) to N 3 LO. However, the adaptation of any of these methods to N 3 LO in full generality is a formidable task. Another avenue is to focus on some specific differential observables at N 3 LO. To this end we introduced the Higgs-differential method [33] and studied its feasibility in obtaining differential distributions without a sophisticated NNLO subtraction method. The method is fully-differential in all components of the Higgs momentum and its decay products, extending ideas first used to obtain rapidity distributions at NNLO [34,35]. At the same time, it treats the additional QCD radiation inclusively (i.e. it integrates over the unrestricted phase space of all final state partons). This enables us to compute differential Higgs boson observables, i.e. the Higgs boson rapidity and transverse momentum distributions, as well as any distributions of decay products. At a later stage the results obtained with this method may serve as a key ingredient to a fully differential computation, for example in combination with q T -subtraction [25].
A complete analytical computation of N 3 LO Higgs differential cross sections presents a formidable challenge. Key ingredients for such a computation are analytic expressions for Feynman integrals in kinematic limits that can serve as boundary data as well as counter terms that render the Higgs differential cross section finite in all its physical limits.
In this article we achieve a significant step towards this computation by performing an expansion of the complete partonic cross sections around the production threshold of the Higgs boson. The analytic information we obtain is comprised of the first two terms of in the threshold expansion and represents the foundation of a complete future computation. Furthermore, we explore the phenomenological potential of threshold expansions to approximate differential Higgs boson observables. Computations at very high orders stress the reliability of conventional tools for higher order computations to their limits. We identify a pitfall with the standard treatment of parton distribution functions within the framework of LHAPDF and present our solution. Furthermore, we report on the complete computation of contributions to the Higgs differential N 3 LO corrections involving explicit logarithms of the perturbative scale.
The threshold prediction for the Higgs boson rapidity distribution at N 3 LO was also obtained in [36][37][38] in the double threshold expansion. The logarithmic contributions to the tranverse momentum spectrum for small transverse momenta were also obtained in [39,40] within the framework of SCET. We would like to point out that our method does not rely on the simplification of the definition of the physical observables in any kinematic limit. This allows us at least in principle to compute the Higgs differential N 3 LO cross sections to any order in the threshold expansion. The results that we obtain in this note are primarily a proof of principle, demonstrating that it is possible, to extend our Higgs differential method to N 3 LO. Consequently, we do not claim any phenomenological significance of the distributions presented in the following. This is applies in particular to the transverse momentum distributions that we show in section 4. It is of course a well known fact that transverse momentum distributions for Higgs production in gluon fusion in the effective field theory have a fairly limited range of applicability. In the low end of the transverse momentum spectrum, the cross section is dominated by large logarithms of the transverse momentum that need to be resummed to all orders in perturbation theory to obtain stable predictions. At the high tail of the transverse momentum spectrum the cross section becomes sensitive to finite quark mass effects and the accuracy of the infinity top mass approximation decreases rapidly. This can be countered by systematically adding corrections to the infinite top-mass limit to the calculations. A detailed discussion of both of these issues is premature at this point, as our first goal is to obtain a fixed order prediction for the Higgs differential distributions in the infinite top-mass limit through N 3 LO in QCD. Once this is achieved, it will be a separate issue to systematically improve upon this result and combine it with known important effects as the ones mentioned above.
This article is organized as follows. In section 2, we review our method of Higgs-differential calculations. In section 3 we describe how threshold expansions can be performed for partonic Higgs-differential cross sections and obtain analytic results for the N 3 LO coefficient function. In section 4 we critically asses the quality of results obtained with a threshold expansion for differential distributions at NNLO. On the way to obtain numerical results for distributions at N 3 LO we encounter in section 5.1 an issue that arises when interpolating parton density functions. Naive usage of LHAPDF leads to a drop in accuracy when computing N 3 LO results due to a lack of smoothness in the default interpolation routines. We illustrate our findings and discuss our way of avoiding this issue. Yet another crucial ingredient for N 3 LO Higgs-differential phenomenology are the complete contributions due to explicit dependence on the perturbative scale which we obtain in section 5.2. Then, in section 5.3 we demonstrate the impact of the newly obtained approximations to the N 3 LO coefficient functions on the rapidity distribution of the Higgs boson and discuss their validity. Finally, we conclude in section 6.

Higgs Differential Cross Sections
In this section we briefly review the definition of Higgs differential cross sections introduced in ref. [33]. Within this framework we consider scattering process of two protons that produce at least a Higgs boson. Proton(P 1 ) + Proton(P 2 ) → H(p h ) + X, (2.1) P 1 and P 2 are the momenta of the colliding protons and p h the momentum of the Higgs boson. The framework of Higgs differential cross sections allows to compute the scattering probability for any observable that is solely dependent on the four momentum of the Higgs boson. Typically such observables are related to the rapidity Y , transverse momentum p T and mass m h of the Higgs boson.
E, p z and p T are the energy of the Higgs boson, its momentum along and transverse to the beam axis in the laboratory rest frame, respectively. The master formula for a Higgs differential cross section for an observable O is then given by Here, we employed the parton model and factorization of long and short range interactions into parton distribution functions f i (x) and partonic differential cross sections. The momenta of the colliding partons are related to the proton momenta by p 1 = x 1 P 1 and p 2 = x 2 P 2 = τ x1z P 2 . The variable τ is given by The sum over i and j ranges over all contributing partons. The variables x, λ and φ parametrize the momentum of the Higgs boson. φ is the azimuthal angle of the Higgs boson with respect to the collision axis. x and λ are related to the more familiar Higgs boson rapidity and transverse momentum by Here,x = 1 − x,λ = 1 − λ andz = 1 − z. The partonic Higgs differential cross section is given by The observable we are interested in is specified by the measurement function J O (x 1 , z, x, λ, φ, m 2 h ), that filters the regions we are interested in. Assume we are interested in computing the probability for a Higgs boson to be produced in the rapidity interval Y ∈ [1, 2] the measurement function would take the form In ref. [33] the partonic Higgs differential cross sections were computed in heavy quark effective theory for the gluon fusion production mode to NNLO in QCD perturbation theory in terms of analytic functions of the variables z, x and λ. Higgs differential cross sections can be easily combined with subsequent decays of the Higgs boson in order to allow for the prediction of fiducial cross sections for Higgs boson decay products as demonstrated in ref. [33].

Threshold Expansion for Higgs-differential N LO
In this section we present the analytic computation for the first and second term of the threshold expansion of the partonic N 3 LO coefficient functions. In order to derive our results we strongly rely on techniques recently developed in refs. [11,41,42,42,43]. We begin by clarifying the ingredients for our partonic coefficient functions. Next, we set-up the notation for the partonic phase space integrals we perform to obtain Higgs differential partonic cross sections. Finally, we explain how threshold expansion for Higgs differential cross sections can be performed at the integrand level.

Setup of the Calculation
The gluon fusion cross section in the standard model is mainly mediated by a top quark loop, coupling the gluons to the Higgs boson. For the production of a Higgs boson with the observed mass of m H = 125GeV with a transverse momentum below the top-pair threshold, the process can be safely described in the limit where the top quark is infinitely heavy. In this limit, the top quark can be integrated out, which induces an effective theory that directly couples the gluons to the Higgs boson via an effective dimension five operator. The Lagrangian for this effective theory is given by, where H is the Higgs field, G µν a is the gluon field strength tensor and L SM,5 denotes the standard model Lagrangian with N f = 5 light flavors. The bare Wilson coefficient C 0 is obtained by matching the effective theory to the standard model in the limit of an infinitely heavy top quark [44][45][46][47]. The inclusive cross section σ P P →H+X has been computed at NLO [3,4,48] as well as at NNLO [6][7][8].
Within the effective theory, we can write the Higgs-differential partonic cross section as, Dividing out the Born cross section,σ we can write the partonic coefficient functions as, The initial state dependent prefactors N ij are given by Here, g, q andq indicate that the initial state parton is a gluon, quark or anti-quark respectively. dΦ H+m is the phase space measure for the production of a Higgs boson and m partons and is explained in more detail below. M h s tends to zero as we approach the production threshold and the partonic center of mass energy becomes equal to the Higgs boson mass. In this work we perform a systematic expansion of the partonic coefficient functions around the production threshold.
We separated the leading term in the expansion that is indicated by the superscript (SV ) and that is commonly referred to as the soft-virtual contribution. This particular term is singular asz → 0 and acts as a distribution on the measurement function and the parton distribution functions as we explain below. All higher power terms depend on z in the form of polynomials of logarithms log(1 − z). The individual terms η (n,i) ij (z, x, λ) depend on the threshold variablez in the form of polynomials of logarithms of the form log(z). In this article we obtain the first and second term in the threshold expansion for all required partonic coefficient functions.
The purely virtual matrix elements are independent of the expansion parameter and were computed in refs. [49]. Matrix elements with two loop and one emission were computed in refs. [42,50,51] and recomputed for the purpose of this article and are known to all orders inz.
In order to obtain the required matrix elements with two or three additional partons in the final state we followed the techniques developed in refs. [11,42,42,43]. We generate Feynman diagrams with QGRAF [52] and perform spinor and color algebra with a private c++ code based on GiNaC [53]. Next, we perform a threshold expansion on the integrand level and subsequently integrate out the momenta of all radiation produced in addition to the Higgs boson. We discuss this step in greater detail below.
Our coefficient functions contain single poles in the variablesz, x and λ. These poles correspond to kinematic singularities of the Higgs boson cross section where the kinematic degrees of freedom degenerate. Explicitly, they correspond to vanishing transverse momentum of the Higgs boson or vanishing virtuality of the system of all radiation produced in association with the Higgs boson. Specifically, before expansion in the dimensional regulator these singularities are of the form where the coefficients a i are small integer numbers. When we compute Higgs differential cross sections as in eq. (2.3) we integrate in the variables z, λ and x and the singularities may lie within our integration range, depending on the observable under consideration. For example, to compute the inclusive cross section we integrate the variables x and λ within the interval [0, 1]. Consequently, regularization of these divergences is required and we proceed as outlined in ref. [33] (and repeated in appendix A). The analytic computation of the first two terms in the threshold expansion of the N 3 LO coefficient functions represents one of the main results of this article. We present this result in Mathematica readable form as an ancillary file submitted together with the arXiv version of this article.

Higgs Differential Phase Space
The integration measure for the production phase space of a Higgs boson and m additional partons is given by We want to develop a parametrization of the above phase space that allows us to compute observables that are differential in the Higgs boson four momentum. Consequently, it seems natural to separate the integration over the momenta of the Higgs and the final state parton momenta. We can achieve this by inserting a unity This identity allows us to write the H plus m parton phase space measure as an integral over a phase space measure for two massive particles and m massless partons.
We choose the rest frame of the initial state partons We introduce the following definition of Lorentz invariant scalar products.
With this we have s = s 12 . We refer to the last component of our d-dimensional vectors as z components and all other space like components as transverse components. The z-axis is parallel to the collision axis of the incoming partons. In this frame we can now express the energy-and z-component of any vector p i in terms of scalar products of this vector and the two incoming momenta.
We start by parametrizing the phase space for two massive particles.
In general we can write Integrating out p h and using the above parametrization we find We exploit the on-shell condition of k to perform the k ⊥ integration to find . (3.18) and find that Consequently, In order to make our live a little easier we are going to exploit the fact that we will later on integrate over µ 2 and re-parametrize µ 2 = sz 2 λλx: Solving the θ constraints and integrating outλ we find We now can combine the previous result and rewrite eq. (3.11) as The remaining massless parton measure is given by Analytic partonic coefficient functions for Higgs differential cross sections are now obtained by performing the integration over m-parton final state squared matrix elements using the above measure.

Threshold Expansions for Higgs Differential Cross Sections
In this section we describe our method to perform a threshold expansion at the integrand level for the required matrix elements for the N 3 LO coefficient function that involve two or more partons. We start by regarding matrix elements that correspond to purely real radiation diagrams and contain no closed loops and then we address the case of virtual radiations.
Consider as an example the following scalar phase space integral.
In the above picture solid lines correspond to scalar propagators, doubled line to massive external legs and lines crossed by the dashed line represent the on-shell constraints of the phase space integration measure.
As described above, the threshold limit corresponds to the kinematic configuration where all radiation produced in association with the Higgs boson is uniformly soft. It is thus natural to perform the variable transformation p f →zp f [41]. Here, p f indicates the momentum of any final state parton, which can be identified in (3.25) by the momenta crossing the dashed line, namely p 3 , p 4 and p 5 . This rescaling induces a transformation on the phase space measure dΦ 0−3 → z 2d−6 dΦ 0−3 contained in (3.25) given that k →zk.
Performing a series expansion of the integrand yields Every term in the soft expansion of a Feynman integral can be written in terms of a linear combination of soft integrals as already observed in ref. [43,54]. Soft integrals are Feynman integrals that are independently homogeneous under rescaling of the initial momenta p 1 , p 2 or all momenta in the integral simultaneously.
Consider for example what happens to our example integral I (0) as we rescale p 1 → λ 1 p 1 .
The associated rescaling dimension can be easily read off the integral and in the specific case of our example we find γ 1 = −2.
In general, The last line in the above equation indicates a simultaneous rescaling of all momenta in the curly bracket. Note that the respective scaling dimensions γ i depend on the specific integral in question but for simplicity we write them without any argument. We realize that the integrated soft integrals are functions of four Lorentz invariant scalar products s, k 2 , 2p 1 k and 2p 2 k. We can use the scaling behavior of our soft integrals to determine its functional dependence on three of the four scalar products. Consequently, the soft integrals depend on one variable that is invariant under any of the three rescaling symmetries. This invariant cross ratio is given by the dimensionless variable x = k 2 s 2kp12kp2 , introduced in eq. (2.3). Combining these properties we are able to write: Once our integrand is expressed in terms of integrals that only depend on the cross ratio x we can use standard phase space integral techniques to compute these functions. In particular, we employ the framework of reverse unitarity [7,34,35,55,56] to express our differential partonic cross sections in terms of a few soft master integrals. Subsequently, we make use of the method of differential equations [57][58][59] to compute our soft master integrals. The soft expansion greatly simplifies these steps as in the expression above we only need to maintain functional dependence on x. The resulting functions are harmonic polylogarithms [60] in the x variable. Note, that so-far we did not apply the relation among the invariants that arises due to the on-shell constraint for the Higgs boson, p 2 h − m 2 h = s + 2kp 2 + 2kp 3 + k 2 − sz = 0. This relation is inhomogeneous under rescaling the final state momenta k and will introduce sub-leading terms in thez expansion. In ref. [42] this issue was solved by performing a systematic expansion of the on-shell constraint by including it in the reverse unitarity method. Here, we choose to impose the aforementioned on-shell constraint only after reduction to master integrals.
An additional complication arises when loop integrals are part of the phase space integrand. The loop momentum can take arbitrary values that are parametrically smaller or larger compared to the parameterz we want to expand our cross section in. The obstacle is easily illustrated by regarding a single propagator containing a loop momentum l and a rescaled final state momentum is simply not convergent if for example l is uniformly smaller thanz.
To be able to perform a systematic threshold expansion at the integrand level for the partonic cross section we rely on the strategy of regions [61] that allows to consistently treat the problem. In particular, we want to expand contributions with one loop and two additional partons in the final state around the threshold limit. Exactly this issue was addressed in much detail in ref. [43] and we refrain from repeating the procedure here. Once the loop and phase space integrand is expanded we perform a reduction of the loop integrals to loop master integrals. The initial step of performing the reduction of loop integrals allows us to determine the rescaling behavior of the loop integrals under the scaling transformations introduced in eq. (3.28). Next, we embed the loop master integrals again in terms of their Feynman propagator representation into the phase space integration and we subsequently perform the combined loop and phase space integral in the same fashion as the pure phase space integrals discussed above (i.e. via reverse unitarity and differential equations). Again, we benefit from having only to maintain functional dependence on the cross ration x by inferring the dependence of our integral on the other variables from its behavior under scaling transformations.
With the techniques summarized in this section we can perform a threshold expansion of the partonic Higgs differential cross sections at arbitrary order in the strong coupling constant and to arbitrary power inz. Specifically, we perform the computation of tree level partonic cross sections with three partons in the final state and partonic cross section with one loop and two partons in the final state to first and second order in the threshold expansion. Extending the threshold expansion to higher powers is technically challenging and is left for future work. As a result we obtain all ingredients to compute the first and second term in the threshold expansion of the N 3 LO Higgs differential cross section.

Validating the Threshold Expansion for differential Observables at NNLO
In proton collisions with a fixed center of mass energy the probability for two constituent partons to collide can be understood as a function of center of mass energy of the colliding partons. This probability is falling as the center of mass energy of the colliding partons rises. In particular, if the two partons under consideration are gluons this probability is falling faster as a function of their energy than the probability for two quarks or one quark and one gluon. As the main source for the production of a Higgs boson at the LHC are two colliding gluons this implies that there is a kinematic enhancement for the Higgs boson to be produced from a system of gluons with as little energy as possible. The lowest possible energy to produce a Higgs boson is referred to as its production threshold and corresponds to the Higgs boson mass. In general, it can be expected that the bulk of Higgs bosons produced in proton collisions are produced at their production threshold and the cross section to find more energy in the produced system is falling with energy.
In the past this simple kinematic consideration was exploited to derive simplifications for the prediction of Higgs boson cross sections. Perturbative NNLO [6] and N 3 LO [11] corrections were approximated in the form of an expansion around the threshold. Factorization properties of the leading term in the threshold expansion are commonly exploited to perform all order resummation of threshold enhanced terms. In this section we will analyze the performance of an expansion of partonic Higgs differential cross sections around the production threshold.
As the partonic differential cross sections were computed analytically as a function ofz in ref. [33] we can easily perform a threshold expansion a posteriori. In the following we want to study the quality of this approximation as higher and higher terms in the expansion are included for differential observables. In order to do so we compute the rapidity and transverse momentum distribution of the Higgs boson. We perform a threshold expansion for all matrix elements occurring at NNLO and keep lower order matrix elements exact. We then truncate the expansion at different orders and compare with the exact results.
To derive numerical values we numerically perform the remaining integrals in eq. (2.3) over the partonic cross sections in conjunction with MMHT2014 parton distribution functions [62] in a private C++ implementation. We perform renormalization and convolutions with the mass factorization counter terms numerically. We only expand the partonic NNLO matrix elements around the production threshold. This leads to a mismatch in the cancellation of infrared and ultra violet divergences which are treated in the framework of dimensional regularization. Specifically the cancellation of poles in the dimensional regulator is only given up to the respective order in the threshold expansion at which we truncate. Throughout this section we choose the perturbative scale to be µ = m h .
Let us first consider the inclusive cross section produced at a perturbative scale µ = m h . In refs. [63] similar studies were performed for the inclusive cross section at NNLO and our findings agree. We show the inclusive cross section through NNLO in fig. 1 for different truncation orders. The first few terms in the threshold expansion (in red) significantly deviate from the exact result (in blue). After the first five terms the expansion stabilizes and subsequent terms gradually improve the result. The agreement after including five terms in the expansion is fairly good. Further improvement is achieved at a comparably slow rate. This slow convergence of the remaining difference to the exact result can be attributed to explicit divergences of the partonic coefficient functions at the high energy limit z = 0. A similar behavior was observed for the expansion of the N 3 LO corrections to the inclusive cross section in refs. [11,12,63]. In fig. 2 we show the rapidity and transverse momentum distribution. Increasing truncation order of the threshold expansion is indicated by increasingly dark shades of red and the exact result in blue. As for the inclusive cross section we observe that the first few terms display large variations from the full result. After about five terms the expansion stabilizes and adding higher terms shows gradual improvements on the approximation. Note that in order to derive a physical prediction for the first couple of bins in the transverse momentum distribution logarithms of the transverse momentum should be resumed to all orders in perturbation theory. No such procedure was applied here as this is beyond the scope of this article.
The actual quality of the expansion can be studied in more detail by analyzing the relative deviations of the expanded distributions from the full result. In fig. 3 we show the rapidity and transverse momentum distribution normalized to the unexpanded respective distributions. We note that by including only the third term in the expansion the rapidity distribution at NNLO is approximated to a level better than five percent. The transverse momentum distribution is improved as higher terms in the expansion are included. However, even with ten terms in the expansion the overall agreement between the exact result and the expansion is merely at the level of ten percent. At large rapidity the quality of the threshold expansion deteriorates as the Higgs boson is produced with a larger boost along the beam axis and thus on average more energy is in the final state system. The stark difference between the behavior of the rapidity and of the transverse momentum distribution can be understood by considering the structure of the partonic coefficient functions. The transverse momentum of the Higgs boson is identically zero at leading order as there is now parton produced for the Higgs boson to recoil against. At the kinematic threshold all radiation produced in association with the Higgs boson is soft and does not provide any recoil either. Adding terms in the threshold expansion only gradually builds up the functional dependence of the matrix elements on the transverse momentum. At the same time, the partonic matrix elements are singular as the transverse momentum of the Higgs boson vanishes. The partonic transverse momentum distribution contains kinematic singularities at finite values ofz that are expanded by the threshold expansion. This leads to a slower convergence compared to the rapidity distribution.
We can conclude that the quality of a threshold expansion is subject to the particularities of individual observables. If such an expansion is to be used to approximate cross section predictions a dedicated analysis of the quality of the approximation has to be performed specific to every observable. Complementing threshold expansions with expansions in the rapidity of the Higgs boson would be an interesting way forward and such studies are left for future work.
It is commonly the case that predictions for inclusive cross sections become available prior to predictions for exclusive observables. Suppose this was the case for NNLO Higgs differential cross sections. In this case we could approximate differential cross sections at NNLO by performing a threshold expansion. To further improve the expanded results we could ensure that the inclusive cross section is reproduced by the differential approximation and only shapes are computed by the approximated result. We define the improved Higgs differential cross section to be Here, σ Full and σ Expanded are the inclusive cross section without and with expanding around the threshold limit respectively. σ Expanded P P →H+X [O] is the Higgs differential cross section based on partonic coefficient functions approximated by a threshold expansion.
We show predictions for the rapidity and transverse momentum distribution based on the improved approximation in fig. 4 normalized to their exact counter parts. The effect of the rescaling  is largest for low truncation order where the difference between the inclusive cross section based on the expansion and the full result is largest. We observe that the rapidity distribution is approximated at a level significantly better than two percent including just three terms in the threshold expansion throughout the rapidity interval [0, 3]. The significant deviations of the shape of the transverse momentum distribution based on the threshold expansion are only mildly impacted by rescaling the distribution to the correct inclusive cross section.

Numerical Results for approximate N 3 LO Cross Sections
In the following we discuss the remaining ingredients required to promote our analytical results for approximate N 3 LO differential partonic cross sections to distributions. As in any calculation of hadronic observables we require parton distributions functions as an external input. Special care has to be taken not to introduce artifacts due to the interpolation of the input parton grids. Next we discuss scale dependence of the cross section that we extract to all orders in the threshold expansion. Combining all ingredients obtained, we then show differential distributions.

Curious Encounters with Parton Distribution Functions
We compute predictions for the total cross sections and distributions by performing Monte-Carlo integrals over the remaining variables of the Higgs phase space given in eq. (2.3). This requires numerical values for the parton distributions which are obtained by accessing the grids of standard PDFs [62,[64][65][66][67] through LHAPDF [68]. The parton distributions functions f (x, Q 2 ) are functions of the partonic momentum fraction x and the scale Q 2 . Internally, the parton distributions are stored as finite grids in (x, Q 2 ) space. These grids are interpolated on-the-fly by LHAPDF to provide the value of the parton distribution at the requested point in (x, Q 2 ) space to the user. By default, this interpolation is performed using a log-cubic spline.
In the following, we want to study the numerical evaluation of the soft-virtual term of the rapidity distribution at N 3 LO. In the strict soft limit, the rapidity distribution takes a particular simple form, since the partonic rapidity becomes unity. Thus the hadron level rapidity Y is just a function of the parton momentum fractions, which can be seen from taking the limitz → 0 in eq. (5.1). The differential partonic cross section in that limit is therefore simply the soft-virtual inclusive cross section [9] and we can write the hadronic cross section as, Evaluating the soft-virtual term of the rapidity distribution at N 3 LO using the NNLO sets from NNPDF 3.0 at α s (m z ) = 0.118 obtained using the default LHAPDF setup, we observe an unexpected loss of accuracy. As can be seen in fig. 5, the rapidity distribution displays strong oscillations. Clearly, the soft-virtual limit of the rapidity distribution in eq. (5.1) has no structures that would warrant this oscillatory behaviour. Comparison with lower orders in perturbation theory suggest strongly that these features are numerical artifacts which should be suppressible by more careful numerics. The origin of these numerical artifacts can be understood when analyzing the influence of the interpolator used in LHAPDF to obtain continuous values from the discrete (x, Q 2 ) grids. Figure 5: Soft-virtual term of the N3LO correction to the absolute rapidity distribution evaluated with various ways of interpolating the PDF grid. The blue line is obtained using a linear interpolator, the yellow line was obtained using a log-linear interpolator, the green line was obtained using a cubic interpolator. The solid orange line was obtained using the default LHAPDF setting, which uses a log-cubic interpolator. The solid purple line was obtained using our custom fit to the LHAPDF grid, which is described in the text. Note that the excursions of the blue and yellow line from the central value obtained through the fit were divided by a factor of five to be able to visualize the oscillations. Fig. 5 shows how the oscillations of the rapidity distribution change with different interpolation orders of the parton distributions. Clearly, these oscillations are artifacts of the interpolators lacking smooth higher order derivatives as can be seen from the fact that the magnitude of the oscillations increases when using lower rank splines for interpolation. These artifacts seem to be caused by the appearance of high powers of logarithms of z in the partonic cross section. As can be seen in eq. (5.1), the soft-virtual N 3 LO contribution to the rapidity distribution contains terms of form for n ∈ [0, 5]. In comparison, the soft-virtual NNLO contribution to the rapidity distribution contains only terms for n ∈ [0, 3]. It seems therefore that the n = 4 and n = 5 terms pick up contributions from discontinuous higher moments of the interpolator. This hypothesis can be tested at NNLO. With the default log-cubic interpolator the soft-virtual contributions to the NNLO rapidity distributions are smooth, as expected. Switching however to a log-linear interpolator, we see the same kind of oscillatory behavior appearing in the NNLO rapidity distribution. It is clear that we need a smoother way to interpolate the LHAPDF grids in order to obtain useful predictions at N 3 LO. One way to obtain smooth values from the parton distributions, is to evolve the parton distributions to a fixed value of Q 2 using LHAPDF and fit the resulting grid in x space with an analytic function that is sufficiently smooth. We make an empirical ansatz for the function as, We fit this ansatz to points obtained from evolving the gluon NNLO NNPDF 3.0 grid to a scale Q 2 = (125GeV) 2 , finding a χ 2 /ndof of 1.9 × 10 −7 with the parameters: Evaluating the rapidity distribution with the fitted x-dependence for the PDFs, we obtain the smooth line in fig. 5.
Performing a fit for parton distribution functions in terms of a smooth functions entails the disadvantage that this procedure has to be repeated for every required PDF set and ensuring a high fit quality sufficient for arbitrary observables is a non trivial task. Furthermore, assessing the goodness of such a fit should be subsequently incorporated in the analysis of uncertainties of cross section predictions. As such it seems advantageous to use instead a higher order interpolator for LHAPDF grids. We therefore implement a custom interpolator that is able to interpolate with splines of varying polynomial degree. We test the interpolation with degree six and degree twelve Legendre polynomials. The results can be seen in fig. 6. As one can see in fig. 6, the interpolation with log-polynomial splines at order 6 smoothes the oscillations in comparison to the default LHAPDF setup, however some artifacts still remain. Using log-polynomial splines at order 12 then leads to acceptably smooth results.
The numerical results at N 3 LO in the remainder of the paper are obtained using the order 12 interpolator.
It should of course be noted that although the loss of accuracy cannot be directly observed in inclusive calculations at N 3 LO, we can still expect an effect to appear from the integral over the oscillations. To test this, we integrate the threshold expansion through 37 terms of the inclusive N 3 LO cross section, as obtained in [11], with two different ways of obtaining numeric values for the parton distributions. We integrate the cross section with values for NNPDF 3.0 directly taken from LHAPDF and compare with integrating the cross section using the fit obtained in eq. (5.3). We observe a deviation of about 1.3% of the full N 3 LO coefficient. We can therefore conclude that the effect on the total cross section through N 3 LO is negligible.
We want to stress here the generality of the appearance of this loss of accuracy. Even though we show here numbers obtained using a particular PDF provider, we have investigated all common PDF sets from the big collaborations and find that these artifacts appear for any PDF set. We should also point out that the appearance of these oscillations is not specific to our calculation. These features had not been observed before, as the interpolators provided by default are smooth enough for NNLO calculations, however any calculation at N 3 LO that relies on LHAPDF will be susceptible to this effect. Clearly, it is desirable to study the effect of interpolator choices in more detail and to provide a flexible means of interpolating LHAPDF grids at sufficiently high orders to obtain smooth predictions for future N 3 LO phenomenology.

Exact Scale Variation at N 3 LO
After ultraviolet (UV) renormalization of the coupling constant and the Wilson coefficient, and after suitable redefinition of the parton distribution function, the Higgs differential cross section takes its final and finite form. The coefficient of α 3 S can be written as correspond to the Laurent series coefficients of the sum of UV renormalization counter term and mass factorization counter term. They are constructed in the usual way in terms of lower order cross sections and universal anomalous dimensions and splitting functions [33]. The renormalized coefficient function is finite as the residual poles of the partonic coefficient function and the counter terms cancel. Consequently we find It is thus easy for us to construct these coefficients explicitly. Utilizing the above identity we may write the finite term of the N 3 LO coefficient function as With this all contributions explicitly depending on the perturbative scale µ of the N 3 LO coefficient functions are known. Additional dependence on the perturbative scale µ arises due to the multiplication of the partonic coefficient functions with the Wilson coefficient and due to the dependence of the strong coupling constant, the Wilson coefficient and the parton distribution functions on the perturbative scale. We now present the impact of the all contributions of the α 3 S coefficient on the Higgs differential cross sections. Specifically, we compute contributions to the rapidity distribution for the Higgs boson given by all ingredients that explicitly contain a logarithm L µ . We include the partonic coefficient function at N 3 LO asη (3,0) ij, RGE (z, x, λ, L µ ) =η as well as all contributions to the α 3 S coefficient of the cross section containing renormalization group logarithms that arise due to the multiplication of the Wilson coefficient with lower order coefficient functions. We thus obtain all contributions to the N 3 LO correction to the Higgs differential cross section involving explicit RGE logarithms. We show numerical impact of the contributions involving RGE logarithms on the N 3 LO corrections to the rapidity distribution of the Higgs boson in fig. 7. We choose µ = m h 2 as a central scale which results in the prediction given by the solid line. The red bands correspond to a variation of the perturbative scale in the interval µ ∈ m h 4 , m h . The contribution is monotonously rising as we increase the perturbative scale. At µ = m h the argument of the RGE logarithm is one and the contribution considered here is identically zero. The corresponding inclusive cross section agrees for all scales with the results of ref. [11]. The contribution presented here on its own does not allow for an improved prediction at N 3 LO since the finite coefficient functions without any RGE logarithms are still missing. However, it represents one more essential stepping stone towards Higgs differential cross sections at N 3 LO.

Numerical Results for approximate differential Distributions at N 3 LO
In section 4 we discussed the phenomenological implications of performing a threshold expansion for Higgs differential cross sections at NNLO. The findings clearly indicate that several terms in the expansion are required. Particularly, predictions made by performing the threshold expansion to only the first or second order displayed sizable deviations from the true result. Nevertheless, in section 3 we went on to compute the first and second term in the threshold expansion of the N 3 LO coefficient function. The main motivation is that this result provides key ingredients for the full analytic computation of the N 3 LO coefficient functions. Furthermore it represents the complete soft counter term for the Higgs differential cross section at N 3 LO. In this section we will demonstrate that the same pattern as observed for the first two terms in the expansion at NNLO proliferates at N 3 LO. We implemented the analytical results we obtained for the partonic coefficient functions in terms of the first and second term in the threshold expansion into a private c++ code. Furthermore, we combine our new results with our computation of the exact scale variation contributions obtained in section 5.2. We show our results for the N 3 LO corrections to the rapidity distribution of the Higgs boson in figure 8a. Results including only the first term in the threshold expansion are shown in blue and including also the second term in red. The bands correspond to variations of the perturbative scale around the central value µ = m h 2 by a factor of two. It is evident that the two predictions based on the first and second order expansion wildly differ which confirms our expectation from the analysis at NNLO. While the scale variation of the correction to the rapidity distribution based on the leading term in the threshold expansion is monotonously increasing with the scale the second order approximation is not. The negative contributions arising from the explicit RGE logarithms are, depending on the exact value of the scale, compensated by the positive contributions arising from the N 3 LO coefficient function and the Wilson coefficient. Clearly, the scale variation can in no way describe the uncertainty due to truncation of the threshold expansion after a finite number of terms. If we were to derive phenomenological predictions from the threshold expansion at N 3 LO, especially with only so few terms, we would have to carefull study the progression of the threshold expansion and derive a measure of uncertainty from e.g. analyzing the threshold expansion at an order where the full result is known, similar to what was done in [12] for the inclusive cross section.
In figure 8b we combine the predictions for the corrections to the rapidity distribution at N 3 LO based on the first and second term in the threshold expansion with lower order results (in red). Exact lower order results are shown for LO, NLO and NNLO in green, yellow and blue respectively. We observe a fairly large impact of the approximate N 3 LO corrections on the rapidity distributions.
The inclusive cross section obtained with our current next-to-soft coefficient functions differs significantly from the inclusive cross section obtained in ref. [11]. Our approximate results show large differences between the two newly computed terms. This confirms our NNLO analysis that demonstrates that an approximation based only on the first and second term in the threshold expansion is insufficient in order to improve on currently existing phenomenological predictions. Further subleading terms in the threshold expansion or a complete computation are required. Improvements towards this goal will be part of future work.

Conclusions
In this article we achieve several key steps towards predicting differential observables to N 3 LO in QCD perturbation theory. We illustrate how a systematic expansion around the production threshold of Higgs-differential cross sections can be performed to arbitrary order. We apply this method to obtain the first and second term in the threshold expansion of the N 3 LO coefficient functions in analytical form. This analytic data represents a corner stone of a complete N 3 LO calculation as it constitutes the complete soft limit of the cross section and contains vital boundary information for the computation of master integrals via differential equations. Furthermore, the obtained information may in future work serve as data to extract anomalous dimension for the resummation of logarithms in the transverse momentum of the Higgs boson.
Furthermore, we analyze the performance of a threshold expansion for the Higgs-differential cross section. We start by considering corrections at NNLO, for which also the exact result is known. Analyzing the analytic structure of the NNLO coefficient functions shows that their threshold expansion is convergent within the unit interval of the threshold parameterz. When studying the inclusive cross section in the threshold expansion, we observe that a oscillatory behaviour of the series when only including the first three powers inz. The series then stabilizes within inclusion of the first five coefficients, resulting in a difference of 3% compared to the full NNLO result. Further improvement due to including even higher order terms is comparably slow.
Next, we analyze the quality of differential predictions obtained with threshold expansions. The rapidity distribution of the Higgs boson at NNLO displays similar behaviour as the inclusive cross section. Including about five terms stabilizes initial oscillatory pattern and leads to good approximations of the full result. The distribution starts to deviate from the full result at high rapidities as more and more energy is required in the final state. The second observable we analyze is the transverse momentum distribution of the Higgs boson. The quality of the approximation obtained including the same amount of terms in the threshold expansion as for the rapidity distribution is greatly reduced. While including higher and higher terms in the expansion is improving the approximation the convergence is so slow that even with ten terms in the expansion the deviations from the exact result are the level of ten percent. In general, the approximations based on threshold expansions can be improved by normalizing differential cross sections such that their cumulant reproduce the exact inclusive result. Our analysis at NNLO shows that the threshold expansion for Higgs-differential cross sections can be a powerful tool. The quality of the approximation has to be carefully assessed for every observable under consideration. Even for comparatively inclusive observables as the total cross section or the rapidity distribution several terms in the threshold expansion are required to obtain a reliable approximation.
We study the numerical impact of the newly obtained terms in the threshold expansion of the N 3 LO coefficient function. The resulting rapidity distribution displays a similar pattern as we observed for the corresponding NNLO coefficient function. For improved phenomenological predictions more terms in the threshold expansion are required. Already now we obtain the full corrections at N 3 LO due to terms with explicit dependence on the perturbative scale.
When computing corrections at N 3 LO to the rapidity distribution of the Higgs boson we observe that the widely used framework for parton distribution functions LHAPDF needs to be modified. The routines used by the tool to interpolate underlying grids for the parton distributions are insufficient to produce smooth distributions at N 3 LO. As a consequence we observe an oscillatory pattern that is modulating the N 3 LO correction to the rapidity distribution obtained with the soft-virtual approximation. We advocate to implement a log-polynomial interpolator of order twelve or a smooth fitting procedure. DE-AC02-76SF00515. BM is supported by the European Commission through the ERC grant HICCUP.

A Regularization of Coefficient Functions
Consider a function f (x) = x −1+aǫ f h (x), for some integer a and with f h (x) holomorphic around x = 0. We are interested in integrating the function over a test function φ(x) on the range [0, 1]. In the case of our Higgs-differential cross section, the test function φ(x) corresponds to the product of the parton luminosity and the measurement function. We can explicitly subtract the divergence at x = 0 and integrate by parts to obtain We now want to give an expression for the partonic cross section that is finite even if all inclusive integrations are performed. To this end we define in a slight abuse of notation, Here the δ distribution is to be understood as acting only on the test function and not on its coefficient in the square bracket. It is easy to see that f s (0) integrates to zero. We can therefore regulate the integrand f (x) by subtracting f s (0), so that every term of its ǫ expansion can be integrated numerically.