Higgs boson production at the LHC using the $q_T$ subtraction formalism at N$^3$LO QCD

We consider higher-order QCD corrections to Higgs boson production through gluon-gluon fusion in the large top quark mass limit in hadron collisions. We extend the transverse-momentum ($q_T$) subtraction method to next-to-next-to-next-to-leading order (N$^3$LO) and combine it with the NNLO Higgs-plus-jet calculation to numerically compute differential infrared-safe observables at N$^3$LO for Higgs boson production in gluon fusion. To cancel the infrared divergences, we exploit the universal behaviour of the associated $q_T$ distributions in the small-$q_T$ region. We document all the necessary ingredients of the transverse-momentum subtraction method up to N$^3$LO. The missing third order collinear functions, which contribute only at $q_T$ =0, are approximated using a prescription which uses the known result for the total Higgs boson cross section at this order. As a first application of the third-order $q_T$ subtraction method, we present the N$^3$LO rapidity distribution of the Higgs boson at the LHC.


Introduction
The most straightforward and successful (as well as systematically improvable) approach to calculations for processes at high-momentum scales M in QCD is a perturbative expansion in the strong coupling α s (M 2 ). Cross sections are written as a series expansion in the parameter α s and an improvement in accuracy is obtained by calculating an increasing number of coefficients in the series. Until a few years ago, the standard for such calculations was next-to-leading order (NLO) accuracy. Recent years have seen a number of next-to-next-to-leading order (NNLO) results for many important processes of interest, such that the emerging standard for precision calculations relevant for LHC phenomenology is the second non-trivial order in the strong coupling α s .
Reducing the theoretical uncertainties remains one of the main motivations for the extension from NLO to NNLO accuracy. This is particularly relevant in two distinct situations. Firstly, NNLO corrections are mandatory for those processes where NLO corrections are comparable in size to the leading order (LO) contribution, both to establish the convergence of the perturbative expansion and to obtain reliable predictions. Secondly, many benchmark processes demand theoretical predictions with the highest possible precision to be able to fully exploit the extraordinary experimental precision that is achievable for this class of processes. Such "standard candles" are not only indispensable tools in detector calibration but also allow for a precise extraction of Standard Model (SM) parameters and parton distribution functions (PDFs).
Extending the perturbative accuracy of QCD calculations to one order higher implies developing new methods and techniques to achieve the cancellation of infrared (IR) divergences that appear at intermediate steps of the calculations. The past few years have witnessed a great development in NNLO subtraction prescriptions. The transverse momentum (q T ) subtraction method [1,2,3], the N -jettiness subtraction [4,5], projection-to-Born [8], residue subtraction [6,7], and the antenna subtraction method [9] have all been successfully applied for LHC phenomenology.
However, in view of the impressive and continuously improving quality of the measurements performed at the LHC, even NNLO accuracy is in some cases not sufficient to match the demands of the LHC data. Typically, these are processes in which the size of the NLO corrections are comparable with the LO, and where the NNLO corrections still exhibit large effects such that the size of the theoretical uncertainties remains larger than the experimental uncertainties.
This motivated a new theoretical effort to go beyond NNLO to include the next perturbative order: the next-to-next-to-next-to-leading order (N 3 LO). Sum rules, branching fractions [10] and deep inelastic structure functions [11] have been known to this order for quite some time. At present, the only hadron collider observables for which N 3 LO QCD corrections have been calculated are the total cross section for Higgs boson production in gluon fusion [12,13] and in vector boson fusion [14]. First steps have been taken towards more differential observables by computing the first N 3 LO threshold expansion terms to the Higgs boson rapidity distribution in gluon fusion [15]. Moreover, the projection-to-Born method has been most recently extended to compute fully differential distributions to N 3 LO, with a proof-of-principle calculation [16] of jet production in deep inelastic scattering.
In this paper we extend the q T subtraction method at N 3 LO to compute Higgs boson production differentially in the Higgs boson rapidity at N 3 LO accuracy. The paper is organized as follows: in Sec. 2 we recall briefly the main ideas of the q T subtraction formalism and we present the necessary ingredients up to N 3 LO, specifying which elements are known analytically and identifying the missing coefficients at N 3 LO. In Sec. 3 we present a prescription for approximating the missing collinear functions at N 3 LO based on the unitarity property of the integral of the transverse momentum distribution. In Sec. 4, we apply the q T subtraction formalism at N 3 LO to produce differential distributions in the rapidity of the Higgs boson. To validate our approach, Sec. 4.1 quantifies the quality of the approximations by repeating them at NNLO, where all of the ingredients to q T subtraction are known. We assess the magnitude of different sources of systematic uncertainties at N 3 LO in Sec. 4.2, yielding final results for the N 3 LO Higgs boson rapidity distribution and the associated theoretical uncertainty in Sec. 4.3. Finally, in Sec. 5 we summarize our results.

The q T subtraction formalism at N 3 LO
This section is devoted to present briefly the transverse-momentum subtraction formalism to N 3 LO in perturbative QCD. The method is illustrated in its general form and special attention is paid to the case of Higgs boson production through gluon-gluon fusion. The q T subtraction formalism presented in this section is the third-order extension of the subtraction method originally proposed in Refs. [1,2,3].
We consider the inclusive hard scattering reaction where h 1 and h 2 denote the two hadrons which collide with momenta p 1 and p 2 producing the identified colourless final-state system F , accompanied by an arbitrary and undetected final state X. The colliding hadrons have centre-of-mass energy √ s, and are treated as massless particles The observed final state F consists of a generic system of non-QCD partons composed of one or more colour singlet particles (such as vector bosons, photons, Higgs bosons, Drell-Yan (DY) lepton pairs and so forth) with momenta q µ i (i = 3, 4, 5, . . . ). The total momentum of the system F is denoted by and the kinematics of the system can be expressed in terms of the total invariant mass M , the transverse momentum q T with respect to the direction of the colliding hadrons (omitting the azimuthal dependence), and the rapidity in the centre-of-mass system of the hadronic collision, Y , The fully differential hadronic cross section can therefore be written as where dσ ab is the differential partonic cross section, ξ 1 , ξ 2 are the partonic momentum fractions and f c/h the distribution function for finding parton c in hadron h. Since F is colourless, the LO partonic cross section can be either initiated by qq annihilation, as in the case of the Drell-Yan process, or by gluon-gluon fusion, as in the case of Higgs boson production. In the case of the Born cross section, the kinematics of the colour-neutral system F is fully constrained such that In order to explain the basic idea of the subtraction formalism, we first notice that at LO, the transverse momentum q T of the final state system F is identically zero. Therefore, as long as q T > 0, the N n LO QCD contributions (with n ≥ 1) are given by the N n−1 LO QCD contributions to the F + jet(s) final state. Consequently, if q T > 0 we have: where the notation N n LO stands for: N 0 LO =LO, N 1 LO =NLO, N 2 LO =NNLO and so forth. Equation (4) implies that if q T > 0 the infrared (IR) divergences that appear in the computation of dσ F N n LO | q T >0 are those already present in dσ F +jet(s) N n−1 LO . Therefore, provided that the IR singularities involved in dσ F +jet(s) N n−1 LO can be handled and cancelled with the available subtraction methods at N n−1 LO, the only remaining singularities at N n LO are associated with the limit q T → 0 and we treat them with the q T subtraction method. Since the small-q T behaviour of the transverse momentum distribution is well known through the resummation program [17] of logarithmicallyenhanced contributions to transverse-momentum distributions, we can (in principle) exploit this knowledge to construct the necessary N n LO counterterms (CT) to subtract the remaining singularity at q T = 0, thereby promoting the q T subtraction method proposed in Refs. [1] to N n LO.
The generic form of the q T subtraction method [1] for the N n LO cross section is where the symbol "⊗" denotes convolutions over the momentum fractions and the flavour indices of the incoming partons and is explicitly defined as The counterterm dσ F CT N n LO constitutes the contribution to the N n LO cross section which cancels the divergences of dσ F +jet(s) N n−1 LO in the limit q T → 0 and renders the term in square brackets finite for all values of q T . The n-th order counterterm can be written as where we note that the dependence of the function Σ F N n LO (q T ) on the transverse momentum q T is not kinematically related to the Born-level process.
The functions Σ F N n LO (q T ) and H F N n LO correspond to the n-th order truncation of the perturbative series in α s of the functions where the labels a and b stand for the partonic channels of the N n LO correction that are mapped to that the Born cross section. The function Σ F (q T ) embodies all the terms of the form log(q 2 T /M 2 ) that are divergent in the limit q T → 0 and reproduces the logarithmically singular behaviour of dσ F +jet(s) in the small-q T limit. Terms proportional to δ(q 2 T ) as well as IR finite terms are absorbed in the perturbative factor H F . The hard coefficient function H F N n LO thus encodes all the IR finite terms of the n-loop contributions.
According to the transverse momentum resummation formula [2] and using the Fourier transformation between the conjugate variables q T and the impact parameter b, the perturbative hard function H F and the corresponding counterterm are obtained by the fixed-order truncation of the identity The large logarithmic corrections are exponentiated in the Sudakov form factor S c (M, b) of the quark (c = q,q) or of the gluon (c = g), which has the following expression: where the functions A and B permit a perturbative expansion in α s : Explicit expressions for the coefficients A (n) g and B (n) g that are relevant for Higgs production are collected in Appendix A up to n = 3. In particular, we also give the B The analytical form of the function Σ F ; (3) in Eq. (8) can be obtained by expanding Eq. (10) to the corresponding matching order. The full analytical formula for Σ F is resummation scheme independent order by order in the strong coupling constant. Therefore, the logarithmic singular behaviour for Σ F at q T → 0 at each given order in α s does not depend on the resummation scheme, and can be validated against the behaviour of the fixed-order results at small q T . To fully account for the logarithmically enhanced terms at a given order requires a sufficient depth in the resummation accuracy prior to its fixed-order expansion in Eq. (8). Specifically, the LO Higgs boson q T distribution receives singular contributions from up to NLL (next-to-leading-logarithm) resummation [33,34], the NLO Higgs boson q T distribution requires the expansion of NNLL resummation [24,35,36,37], and the NNLO Higgs boson q T distribution has been recently validated against the singular contributions from N 3 LL resummation [38,39].
The structure of the symbolic factor denoted by H F C 1 C 2 cc;ab in Eq. (10), depends on the initial-state channel of the Born subprocess and is explained in detail in Refs. [19,20]. Here we limit ourselves to the case in which the final state system F is composed of a single Higgs boson, F ≡ H, in which case, where H H g is the hard-virtual function and respectively C g a and G g a the gluonic helicity-preserving and helicity-flipping hard-collinear coefficient functions.
The gluonic hard-collinear coefficient function C g a (z; α s ) (a = q,q, g) has the following perturbative expansion In contrast, the perturbative expansion of the helicity flip hard-collinear coefficient function G ga , which is specific to gluon-initiated processes, starts only at O(α s ), and can be expanded as [19,20] G g a (z; α s ) = ∞ n=1 α s π n G (n) g a (z) .
The IR finite contribution of the n-loop correction terms to the Born subprocess is contained in the hard-virtual function (which does not depend on z 1 or z 2 ), Using Eqs. (10) and (13), then, after integration over b and dropping the renormalisation group predictable terms that are produced by evolving α s to a common scale (i.e. setting µ F = µ R = M ), we obtain the resummation scheme independent Note that in the literature, it is often the rapidity-integrated variant H H gg←ab (z) that is quoted which is related to H H gg←ab (z 1 , z 2 ) via the convolution The H H function in Eq. (17) can be expanded perturbatively without approximation to any order in the strong coupling constant α s . Inserting the expansions of the hard functions into Eq. (17), then, Explicit expressions for the known fixed-order coefficients are collected in Appendix A.
The new third-order contribution is given by The second-order helicity-flip functions G g a (z), the third-order collinear functions C g a (z) and the third-order hard-virtual coefficient H H;(3) g are only known in parts or not at all, thereby presenting an obstacle to applying the q T subtraction formalism at N 3 LO. Nevertheless, within the q T subtraction formalism, all these resummation coefficients can be inferred for any hard scattering process whose corresponding total cross section is known at N 3 LO. This point is discussed in detail in Sect. 3.
Although the hard-virtual coefficient H is currently not known in analytical form, parts of it can be inferred from known results in threshold resummation. This relies on the knowledge of the general structure of H F c (to all orders), which relates H to the finite part of the n-loop virtual Matrix Element [20]. To this end, we split H H;(3) g into two pieces, where H can be computed using the corresponding hard-virtual factor C th (3) gg→H [43] from threshold resummation (in the large-m t limit) and the exponential equation that relates hard-virtual coefficients in threshold-and q T -resummation (Eq. (81) of Ref. [20]). We find, with L t = ln(M 2 /m 2 t ) and ζ n denoting the Riemann zeta-function for integer values n (ζ 2 = π 2 /6, ζ 3 = 1.202 . . . , ζ 4 = π 4 /90). Note that we neglect all the third-order terms in the exponent of Eq. (81) in Ref. [20], considering the entire O(α 3 s ) correction (in the exponent) as unknown. However, the full top-mass dependence of H (2) ) represents a single coefficient (of soft origin) belonging to the finite part of the structure of the IR singularities contained in the third-order virtual amplitude of the corresponding partonic subprocess gg → H.
As a consequence, the only missing ingredients to H H; (3) are the functions G (2) ) . The details on their numerical extraction will be discussed in the following section.

The Higgs boson total cross section at N LO
We start this section by reviewing some properties of the hard-scattering function H F cc←ab . This function is resummation-scheme independent, but it depends on the specific hard-scattering subprocess c +c → F . The coefficients H F ;(n) cc←ab of the perturbative expansion in Eq. (9) can be determined by performing a perturbative calculation of the q T distribution in the limit q T → 0. In the right-hand side of Eq. (10), the function H F controls the strict perturbative normalization of the corresponding total cross section (i.e. the integral of the total q T distribution). This unitarity-related property can be exploited to determine the coefficients H F ;(n) cc←ab from the perturbative calculation of the inclusive cross section. In particular, the integral of the full q T spectrum in Eq. (5) must reproduce the inclusive cross section σ F (tot.) , Since the hard-scattering function H F cc←ab is accompanied by δ(q 2 T ), we evaluate the q T spectrum on right-hand side of Eq. (5) according to the following decomposition [2] σ F (tot.) where dσ F (fin.) is directly related to the quantity in square brackets in the right-hand side of Eq. (5) The relation in Eq. (25) is valid order-by-order in QCD perturbation theory [2]. If the perturbative coefficients of the fixed-order expansion of σ F (tot.) , H F and dσ F (fin.) /dq 2 T are all known, the relation (25) has to be regarded as an identity, which can be explicitly checked. Since the fixed-order truncation of dσ F (fin.) /dq 2 T is free of any contribution proportional to δ(q 2 T ), its NLO contribution does not contain the coefficient H F ; (1) , and so forth. Therefore, H F ;(3) can be isolated from the the N 3 LO term in Eq. (25): where α s = α s (µ 2 R ). If all the components on the left-hand side of Eq. (27) are known analytically (as it was the case at NNLO in Refs. [25,26]), the function H F can be extracted exactly in analytical form. At NLO the extraction of the function H F ;(1) is straightforward for Drell-Yan and Higgs boson production. The function H F ;(2) at NNLO (for Higgs (F = H) boson production [25] and Drell-Yan (F = DY ) [26]) can be obtained with a dedicated analytical computation using the analogue of Eq. (27) at NNLO. Since the transverse momentum distributions for H+jet and DY+jet at NNLO are not known analytically, Eq. (27) can be used only numerically to compute H F ; (3) .
As was elaborated on at the end of the previous section, the general structure of the coefficient H F ; (3) is not known in analytic form for any hard-scattering process. Nonetheless, within the q T subtraction formalism, H F ;(3) can be reliably approximated for any hard-scattering process whose corresponding total cross section is known at N 3 LO. As identified in Eqs. (21) and (22), the only missing ingredients to H F ;(3) are the functions G (2) ) . Their contribution to Eq. (27) can be approximated as follows: where the third-order coefficient C N 3 embodies the numerical extraction of the hard-virtual coefficient H (2) ) plus the approximation of the z i -dependent functions by a numerical constant proportional to δ(1 − z i ). The resulting coefficient H ). In other words, the approximation that is made in Eq. (28) is related only to the functions G (2) g a (z) and C g a (z), whose functional dependence on the variable z goes beyond terms proportional to δ(1 − z), and which involves not only gluon-togluon transitions (a = g), but also contributions from other parton species (a = q,q). The latter are not explicitly distinguished in the above approximation, which fully attributes their numerical contribution to the gluon-induced processes.
The method outlined in Eq. (28) to approximate the unknown terms in the hard-virtual function H H gg←ab numerically is not new. It was first used in Ref. [2] in order to compute the second order function H H; (2) gg←ab numerically at NNLO, providing a reasonable estimate of the exact result to better than 1% accuracy. Notice that Eq. (28) ensures that one recovers the total cross section (at N 3 LO in this case) with no approximation. After integration over the transverse momentum q T , Eq. (24) provides the same total integral (numerically in this case) as in the fully analytical case. Even more, for IR-safe observables (at fixed order) where the back-to-back kinematical configuration (q T = 0) is located at a single phase space point (e.g. the q T distribution, the angular separation ∆ϕ γγ between the two photons for a Higgs boson decaying into diphotons, etc.), the fixed order result is also exact, i.e. the integral of the analytical unknown terms in Eq. (28) (which all have q T = 0) is located in one single point of the exclusive differential distributions.
The previous considerations about the approximation underpinning Eq. (28) were regarding the total cross section or differential distributions in which the Born-like configurations belong to one single phase space point. In order to quantify the quality of the approximation proposed in Eq. (28) at the differential level when the Born differential cross section populates the entire differential range, we perform a detailed numerical study of the Higgs boson rapidity Y ≡ y H distribution in Sec. 4.1 at NNLO. Anticipating these results, we find that in the rapidity range 0 ≤ y H ≤ 4 the approximated NNLO result differs by less than 0.2% from the exact NNLO Higgs boson rapidity distribution.

Implementation and setup of the numerical calculations
To extract the value of C N 3 , we first introduce the numerical tools and the calculational setup in this section. We use the same setup for the inclusive and differential predictions presented in Sections 3.2, 4.1, 4.2 and 4.3.
We consider Higgs boson production in proton-proton collisions at a centre-of-mass energy of √ s = 13 TeV. In our computation, we set the Higgs boson mass to M ≡ M H = 125 GeV and the vacuum expectation value to v = 246.2 GeV. The Born process is initiated via gluon-gluon fusion mediated through a top-quark loop, which can be integrated out in the large-m t limit (m t → ∞). In this limit, the production of the Higgs boson is described through an effective gluon-gluon-Higgs boson vertex [45]. The mass of the top quark is taken as m t = 173.2 GeV, which enters in the contributions that have a residual m t dependence (e.g. Eqs. (42) and (23) and effective vertex coefficient corrections at N 3 LO). With the top quark loop replaced by an effective vertex, we consider a five-flavour scheme QCD with all light quarks being massless. We use the central set of the PDF4LHC15 PDFs [46] as implemented in the LHAPDF framework [47] and the associated strong coupling constant with α s (M Z ) = 0.118. Note that we systematically employ the same order in the PDFs (in particular the set PDF4LHC15_nnlo_mc) for the LO, NLO, NNLO and N 3 LO results presented in this paper. The central factorization and renormalization scale is chosen as µ ≡ µ R = µ F = M H /2. The theoretical uncertainty is estimated by varying the default scale choice independently for µ R and µ F by factors of {1/2, 2} while omitting combinations with µ R /µ F = 4 or 1/4, resulting in the common seven-point variation of scale combinations.
As stated in Sec. 3 and in Ref. [1], the computation of the total cross section or differential distributions with the q T subtraction formalism can be separated into two main parts by inserting Eq. (26) into Eq. (25): The contribution dσ F +jet(s) in Eq. (29) is computed with the parton-level event generator NNLOJET which provides the necessary infrastructure for the antenna subtraction method up to NNLO [9]. Processes at NNLO with the structure of dσ F +jet(s) implemented in NNLOJET are: F = H [48], F = γ * , Z [49,50] and F = W ± [51]. In this paper we focus on Higgs production F = H, where the relevant matrix elements in NNLOJET are: (H + 1)-parton production at two loops [52], (H + 2)parton production at one loop [53,54,55] and (H + 3)-parton production at tree-level [56,57,58].
The subtraction formalism that we are applying to Higgs boson production could be easily extended to Z and W ± production [59]. N 3 LO ) are also required. We use the analytical coefficient function for the total Higgs boson cross section that was recently calculated in Ref. [13] and which is available within the public program ihixs 2 [61]. This program is further used to compute any of the analytical total cross-section ingredients required to extract the coefficient C N 3 .
The numerical computation of the integral of the difference dσ Once cancellations between the terms on the right-hand side of Eq. (29) take place, the numerically calculated total cross sections and differential distributions have to be q cut T independent (within the statistical errors) over some range of q cut T . At the lower end of this range, numerical instabilities in dσ F +jet(s) NNLO (arising from the large dynamical range in this calculation) will limit the accuracy of the result, while at the higher end of the range, missing non-logarithmic terms in dσ F CT N 3 LO will start to become significant. The numerical stability of dσ F +jet(s) NNLO at small q T using NNLOJET has been systematically validated for Higgs boson production (with q cut T = 0.7 GeV in Ref. [38]) and Drell-Yan production (with q cut T = 2 GeV in Ref. [39]) at the LHC. In Sections 3.2, 4.1, 4.2 and 4.3, we document numerical results obtained with the q T subtraction formalism using q cut T = (2 ± 1) GeV.

The numerical extraction of C N 3
In the following, we describe the numerical results regarding the extraction of the C N 3 coefficient and the corresponding N 3 LO total cross section. In Fig. 1 Table 1. The solid red central line in Fig. 2, and the associated red band are obtained using this single value.
The numerically extracted C N 3 coefficient allows the total cross section to be computed at N 3 LO using the q T subtraction method, which serves as a closure test of the approach and the approximations used, and allows the impact of uncertainties associated with the numerical evaluation of the ingredients to be quantified. In Fig. 3 we compare the fully analytical N 3 LO Higgs boson total cross section [13] (red dot) and our estimation (red dot with error bar) for three central scales, using q cut T = 2 GeV. The yellow dots with error bar represent our best approximation without the use of the C N 3 coefficient (i.e. C N 3 = 0), that can be considered as the prediction of the q T subtraction method in the case in which the total cross section is unknown (e.g. for Drell-Yan at N 3 LO). The uncertainty bars in the q T subtraction prediction correspond to the statistical errors of the numerical computations and are mainly due to the finite contribution in Eq. (26) at N 3 LO-only. The green crosses and purple squares correspond to our N 3 LO prediction using q cut T = 1 GeV and 3 GeV respectively. Notice that the q cut T variation is performed at N 3 LO-only, while the NNLO cross section is evaluated at fixed q cut T parameter. The NNLO cross section is also shown in Fig. 3 (blue dots) in order to put the size of the N 3 LO corrections in relation to the previous perturbative order. The total cross sections shown in Fig. 3 are reported in Table 2.   Table 1. The error bars for each particular C N 3 point are obtained propagating the statistical uncertainties of the different terms involved in the computation. The red band corresponds to our best estimation for C N 3 obtained with the central scale µ = M H /2 at q cut T = 1 GeV, as detailed in the text.  Table 1: Extracted values of the C N 3 coefficients as a function of the q cut T as shown in Fig. 2 for each scale choice. In bold typeface the C N 3 coefficient (for the case q cut T =1 GeV) which constitutes our best estimation. The uncertainty for each one of the C N 3 coefficients is determined with the customary propagations of the uncertainties. The first column is used to label each particular scale choice used in Fig. 2. [pb]

The rapidity distribution of the Higgs boson
In this section we use the C N 3 coefficient (extracted in Sec. 3.2) to produce differential predictions at N 3 LO. In particular, we present differential results for the rapidity distribution of the Higgs boson. In Sec. 4.1 we first estimate at NNLO the uncertainties introduced in the rapidity distribution by the procedure proposed in Eq. (28). In Sec. 4.2 we present the rapidity distribution at N 3 LO with the estimation of the uncertainties associated to the variation of the q cut T and C N 3 parameters.

The NNLO rapidity distribution
In this section we aim to quantify the uncertainty in the approximation used in Eq. (28). This approximation was first proposed in Ref. [2] for Higgs production at NNLO. Since all the ingredients of the q T subtraction formalism at NNLO are known in analytical form [25], it is possible to quantify the difference induced by the approximation compared to the exact result. This analysis further allows to assess the potential impact of the approximation that could be present at N 3 LO in Sec. 4.2 and 4.3 below. For this quantitative study we consider the collinear functions C where H H;(2) g (δ q T (1) ) is considered as unknown for the present NNLO study. These so-called unknown functions (for this exercise) which depend on the variables z i in Eq. (20) are approximated with a single numerical coefficient C N 2 proportional to δ(1−z 1 )δ(1−z 2 ) (the C N 2 here was labeled as C N in Ref. [2]) in direct analogy to Eq. (28): In Fig. 4 we show the rapidity distribution of the Higgs boson at NNLO computed with the exact q T subtraction (blue hatched band) and the NNLO prediction using the C N 2 coefficient (dot, cross and square points). For this particular example at NNLO, we employ the three-point scale variation: Repeating the analysis performed for Table 1 and Fig. 2, we obtain: C N 2 = 28 ± 1. The numerical value of the C N 2 parameter corresponds to a specific H H;(2) g hard coefficient: which is obtained with the same method that was used to arrive at Eq. (23). Using this C N 2 parameter we can produce differential predictions which are obtained mimicking the strategy that we intend to apply at N 3 LO.  In the lower panel of Fig. 4 we show the ratio to the exact NNLO result, i.e. we present the ratio for each scale. As expected, the approximation presents its best behaviour at central rapidity and the deviation from the exact results is at per mille level throughout the considered rapidity range of |y H | ≤ 4.
We performed at NNLO variations of the q cut T value between 0.1 GeV and 3 GeV, and the NNLO cross sections (and differential distributions) present deviations within a range of size 0.26% (the largest deviation is always observed for the scale choice µ = M H /4). We consider q cut T = 1 GeV enough to proceed at NNLO (and as our reference value), as we can understand from Table 2 at NNLO.

Numerical stability of the N 3 LO rapidity distribution
In this section, we quantify the numerical stability (as well as the involved intrinsic uncertainties) of the Higgs boson rapidity distribution at N 3 LO concerning the q cut T and C N 3 parameters and the statistical uncertainties introduced by dσ H (fin.) /dy H at N 3 LO-only.
In Fig. 5 we show the rapidity distribution at N 3 LO obtained with the q T subtraction method using the C N 3 coefficient determined in Sec. 3.2 (C N 3 = −943 ± 222). The NNLO prediction is always computed with q cut T = 1 GeV. The red band in Fig. 5 shows the size of the seven-point scale variation for q cut T = 2 GeV.  Figure 5: Rapidity distribution of the Higgs boson as computed using the q T subtraction formalism at N 3 LO. All bands include the seven-point scale variation as detailed in Table 1. The red band constitutes our result with q cut T = 2 GeV using the central value for the C N 3 coefficient (C N 3 = −943). The pale yellow band is obtained as the envelope between the prediction at q cut T = 1 GeV and 2 GeV using C N 3 = −943. The black band is computed at fixed q cut T = 2 GeV taking the two extremal values of the C N 3 coefficient according to the uncertainty (C N 3 = −943 ± 222), and performing seven-point scale variation as described in the text.

HN3LO + NNLOJET
The pale yellow band is calculated as the envelope of the scale variation bands for two different values of q cut T : 1 GeV and 2 GeV. Therefore, the pale yellow band in Fig. 5 can be taken as an estimate of the uncertainty due to the variation of the q cut T parameters at N 3 LO. In Fig. 3 (and Table 2), we observed that the total cross section (for the three central scales) is rather stable as a function of the q cut T value. The variations of the N 3 LO cross sections were at the per mille level of accuracy if we consider q cut T = 2 ± 1 GeV, which is far better than the associated statistical uncertainty (see Table 2). The uncertainty estimate due to the q cut T variation performed in Fig. 5, which is differential in the Higgs-boson rapidity, confirms the stability of the total cross section reported in Table 2. The rapidity distribution is almost insensitive to the change in the q cut T parameter in the region where the bulk of the cross section is concentrated (|y H | ≤ 3.6). At large rapidities (|y H | ∼ 4), where the overall contribution to the total cross section is less than 0.5%, we found the largest deviations. Such deviations are mainly related to the numerical uncertainties from dσ H (fin.) /dy H at N 3 LO-only.
Finally, we consider the uncertainty introduced by the statistical errors of the C N 3 coefficient. The black band in Fig. 5 is obtained as the envelope of the seven-point scale variation at q cut T = 2 GeV now considering for each scale the two extremal C N 3 coefficients corresponding to its maximum and minimum statistical deviations: C N 3 = {−1165, −721}. The envelope is therefore taken from a total of 14 rapidity distributions (two extremal predictions for each one of the seven scales). The net effect of this C N 3 variation result in an overall enlargement of the red band at q cut T = 2 GeV. Our final uncertainty estimate in the rapidity of the Higgs boson at N 3 LO is computed as the envelope of three bands: seven-point scale variation only, combined with q cut T variation, and combined with C N 3 variation.

The rapidity distribution of the Higgs boson at N 3 LO
In this section we present our predictions for the Higgs boson rapidity distributions at the LHC, applying the N 3 LO q T subtraction method presented in Sec. 2. The setup of the calculation is summarised in Sec. 3.2.  the band additionally includes the uncertainties due to q cut T and C N 3 as described in Sec. 4.2. Going from LO to NNLO, the scale µ = M H /2 is always at the center of the respective scale variation band in Fig. 6. The central prediction at N 3 LO, on the other hand, almost coincides with the upper edge of the band, as was already observed for the total cross section [12,13], see Table 2 and Fig. 3. Figures 3 and 6 respectively show a substantial reduction in the size of the scale variation band at N 3 LO both in the total cross section and in differential distributions.
In the central rapidity region of |y H | ≤ 3.6, the impact of the N 3 LO corrections on the NNLO result is almost independent of y H with a flat K-factor about 1.034 for the central scale choice. The combined theoretical uncertainty at N 3 LO is at most of ±5% level with respect to the central scale choice. The uncertainty on the y H distribution is reduced by more than a factor of 1/2 by going from NNLO to N 3 LO. The N 3 LO uncertainty band lies fully within the scale variation band at NNLO, exhibiting a stable perturbative behaviour. The only exception is the very high rapidity region, where the q cut T uncertainty becomes the dominant source for the size of the N 3 LO band as shown in Fig. 5.
The N 3 LO corrections to the Higgs boson rapidity distribution have been investigated recently in [15], where an analytic calculation of the leading terms of the threshold expansion of the rapiditydifferential coefficient function was presented. Based on the behaviour of the threshold expansion for the total cross section [12], it is anticipated that the currently known terms in the threshold expansion for the rapidity distribution [15] are not yet sufficient to provide a proper description. Comparing Fig. 6 with the results obtained in Ref. [15], although with different choices of PDFs and scale-variation prescriptions, we observe that the central scale N 3 LO values for the rapidity region y H < 0.5 agree well between the two calculations. Both calculations display a considerable reduction of scale uncertainties going from NNLO to N 3 LO in this central rapidity region. For the rapidity region y H > 1, however, larger differences are observed between the two calculations, where the results using the q T subtraction formalism generally yield smaller N 3 LO corrections (within the NNLO scale uncertainty band).

Conclusions and outlook
In this paper we have performed a detailed study of Higgs boson production at the LHC using the q T subtraction formalism at N 3 LO. We systematically describe the q T subtraction formalism for a generic colourless and massive system F ({q i }) produced at hadron colliders. Fully differential cross sections for this type of final state system are separated into δ(q T ) and q T = 0 contributions. The contribution for q T = 0 is calculated, using a phase space cut-off q cut T , as the difference between F ({q i }) + jet(s) production and q T counterterms. Specifically, we use the NNLOJET package to compute NNLO Higgs-plus-jet production and expand the Sudakov from factor in the hard resummation scheme to the matching order for the corresponding q T counterterms. The contribution at δ(q 2 T ) is further factorized into convolutions of the Sudakov form factor, the hard-virtual function, the helicity-flip coefficient function, the hard-collinear coefficient function as well as the PDFs (Sec. 2). The factorization guarantees that all the process-dependent contributions proportional to a form factor are included in the hard-virtual function, which depends on both initial-and final-state particles. All other factorized contributions only depend on the initial states. Some of the factorized ingredients contributing at δ(q 2 T ) are not known analytically at N 3 LO for the moment. We collect all analytically available contributions and approximate the unknown pieces by a constant coefficient C N 3 which is scale-and process-independent (Sec. 3). Using the available inclusive total cross section for N 3 LO Higgs production and the known pieces from the q T subtraction formalism, we numerically extract the value of C N 3 . By comparing the numerical values for C N 3 using different scales and q cut T setups in the extraction, we conclude from mutually consistent results that C N 3 is independent of the scale choice with a value obtained for µ = M H /2 and q cut T = 1 GeV of C N 3 = −943 ± 222 (Sec. 3.2). As a proof-of-concept implementation of the q T subtraction method at N 3 LO, we calculate the total cross section and rapidity distributions for Higgs boson production at LHC using a new Monte Carlo generator HN3LO [60]. Using the extracted value of C N 3 , we perform a closure test for the inclusive total cross section for three different scale choices and find excellent agreement with the exact results (from ihixs 2 [61]) at the 0.2% level. For the differential rapidity distribution of the Higgs boson, we first study the systematic error from the C N 3 approximation by considering the NNLO calculation and introducing an approximate C N 2 . The NNLO y H distribution exhibit per-mille level agreement between the C N 2 approximation and the exact result, supporting the reliability of the procedure. We calculate the y H distribution at N 3 LO employing a seven-point scale variation and carefully assess systematic errors arising form different q cut T and C N 3 values. Compared to the NNLO y H distributions, we observe a large reduction of theory uncertainties by more than 50% at N 3 LO. The scale variation band at N 3 LO stays within the NNLO band with a flat K-factor of about 1.034 in the central rapidity region (|y H | ≤ 3.6). Both the systematic error analysis and the phenomenological predictions confirm that our calculations at N 3 LO using q T subtraction formalism are well under control. The approximation related to the C N 3 coefficient in our approach can be easily replaced by the full analytical results once available.
With the upcoming larger data set and more accurate measurements of Higgs properties at the LHC, we prepare precise theoretical tools that could match the frontier accuracy of experimental results. More differential properties at N 3 LO involving the Higgs boson and its decay products can be studied using the same framework established in this paper. The current N 3 LO calculation, using the approximation of large top quark mass, attains a level of precision that several other contributions will need to be taken into account for a full study of precision phenomenology [62]: finite top quark mass effects, heavy-light quark interference contributions and electroweak corrections.
where N f is the number of massless QCD flavours and the SU (N c ) colour factors are C A = N c and C F = (N 2 c − 1)/(2N c ). Throughout this paper we always use the hard resummation scheme [20] to report explicit expressions for the perturbative expansion of these individual coefficients. The hard resummation scheme states that all the contributions proportional to δ(1 − z) are associated with the hardvirtual functions H F c . This directly implies that H F c is process dependent. The collinear C ab and G ab functions and the resummation coefficients A c and B c are independent of the final state system F .
The truncation of Eq. (10) at a given fixed order requires the explicit knowledge of resummation coefficients and hard collinear coefficient functions. For F = H at NLO, the knowledge of the coefficients A (1) ga (a = q,q, g) and H H;(1) g are sufficient to compute the inclusive total cross section and differential distributions. Assuming that the Higgs boson couples to a single heavy quark of mass m Q , the first-order coefficient H H;(1) g in the hard resummation scheme is [20] H H;(1) The function c H (m Q ), which depends on the NLO virtual corrections of the Born subprocess, is given in Eq. (B.2) of Ref. [40]. In the limit m Q → ∞, the function c H becomes Therefore, the complete set of coefficients necessary to compute Higgs boson production (in the limit in which the mass of the top quark Q = t is larger than any other scale involved in the process) at NLO are The coefficients A (1) g and B (1) g are process and resummation scheme independent. The collinear functions C (1) ga (a = q,q, g) are process independent, while H H;(1) g depends on the final-state system (F = H). Together, they depend on the resummation scheme in such a way to ensure the resummation scheme independence of Eq. (10) at NLO. In Ref. [24] was shown that the NLO hard-virtual coefficient H F ;(1) c is explicitly related to dσ F LO and to the IR finite part of the NLO virtual correction to the Born cross section.
At NNLO, the coefficients A (2) g and B (2) g are needed [2,20,24], where γ (1) g is the coefficient of the δ(1 − z) term in the NLO gluon splitting function [29,30], which reads The coefficient A g does not depend on the resummation scheme whereas B (2) g in Eq. (40) is valid in the hard resummation scheme and both coefficients are process independent.
The general structure of the hard-virtual coefficients H F c has been established in Ref. [20]. Although H F c is in principle process dependent, Ref. [20] showed it can be directly related in a universal way to the IR finite part of the all-order virtual amplitude of the corresponding partonic subprocess cc → F . The relationship between H F c and the all-order virtual correction to the partonic subprocess cc → F has been made explicit up to NNLO and is based on the definition of universal subtraction operators that cancel the IR divergences of the two-loop (NNLO) virtual corrections to the Born cross section [32]. These universal second-order operators contain an IR finite term of soft origin (δ (1) q T ) that only depends on the initial-state partons [20]. In the case of Higgs boson production, the hard-virtual factor H F =H;(2) g in the large-m t limit (in the hard resummation scheme) is given by [25] where L t = ln(M 2 /m 2 t ). The two-loop scattering amplitude [31] used in the computation of H F =H;(2) g includes corrections to the large-m t approximation. Due to the large size of the expressions for C (2) ab (z), we refrain from explicitly quoting them here and instead refer to Eqs. (37)-(40) of Ref. [20] using the full results of Refs. [25,26]. These collinear coefficients C (2) ab have been independently computed in Refs. [27,28]. At NNLO, in Eq. (13) the first order G (1) ga helicity-flip functions are required which read [19] where C q;q = C F and C g = C A . The first-order functions G (1) ga are resummation-scheme independent and do not depend on the final-state system F .
At N 3 LO, the numerical implementation of Eq. (10) requires the following ingredients: A ga (a = q,q, g) and H The explicit expression of the B c (a = q, g) coefficients in the hard scheme can be computed from Refs. [41,42]. In the particular case of the gluon channel then in the hard resummation scheme, we obtain

B Convolutions at N 3 LO
The numerical implementation of Eq. (10) requires the computation of several convolutions between splitting functions, collinear and helicity-flip functions. In principle, taking the N -moments of the functions involved in the calculation, one can avoid the use of convolutions, since in N -space they correspond to simple products. However, the numerical implementation of Eq. (10) in the Monte Carlo code HN3LO was carried out in the z-space (e.g. as in the codes HNNLO [1], DYNNLO [63], 2γNNLO [64], etc.), and therefore the new third order convolutions have to be calculated as well.
The convolutions in Eqs. (20), (21), (28) and (31) between two functions (f (z) and g(z)) of the the variable z are defined through the following integral (f ⊗ g) (z) ≡ 1 z dy y f z y g(y) .
In the case of processes initiated by gluon fusion, the complete list of third order convolutions to be calculated can be found in Table 3. All the remaining convolutions in Eq. (10) at N 3 LO already contributed to the previous orders and they are regarded as known.
The symbol γ (n) ab in Table 3 denotes the usual splitting functions of n-th order and they contribute to Eq. (10)) since the PDFs have to be evolved from the scale b 2 0 /b 2 to the factorization scale    1, 2). The repeated subindices a and b imply a sum over the parton flavors q,q, g. The first and last subindices denote the partonic channel in which they are contributing, i.e. the convolutions in the first column are used in the gg partonic channel whereas the second (and last) column is for the qg and gq partonic channels. µ F . The first three rows in Eq. (3) were calculated in Ref. [66] and cross-checked with a dedicated computation for the results presented in this paper. The public Mathematica package MT [65] is used to calculate the necessary convolutions (i)-(vi) in Ref. [66], which can be further expressed in terms of harmonic polylogarithms (HPLs) [69] using the Mathematica package HPL [68]. The remaining convolutions in Eqs. (vii)-(xii) of Table 3 were computed for this work. The MT [65] package is not able to solve all the convolutions of weight 3 and 4 that are needed in (vii)-(xii). For instance, the MT package cannot handle convolutions in which their result has to be expressed in terms of multiple polylogarithms (or Goncharov polylogarithms GPLs) [70,67,71] as it is the case when the collinear functions C (2) gj are involved. For those, we have computed the convolutions (vii)-(xii) with a newly developed code Convo, which is able to provide results in terms of GPLs and also can handle terms that are individually divergent, but finite after addition.
where the plus distribution D 0 [1 − z] is defined as usual After performing all the convolutions listed in Table 3, their final expressions (each one of the convolutions) are finite in the domain z ∈ (0, 1). Even more, convolutions evaluated in the domain z ∈ (0, 1) produce results in R. It is possible to write the expressions in Table 3 (after simplifying) in terms of twelve GPLs that are not reducible to polylogarithmic functions of type Li n (z), and cannot be combined (e.g. through the shuffle algebra) with other GPLs in order to produce simpler results. The list of the irreducible GPLs is presented in Table 4. All remaining GPLs appearing in the convolutions of Table 3 can be related to the set given in Table 4 using the results of Refs. [73,68,72] and performing the customary shuffle algebra. The numerical implementation of the GPLs in Table 4 was made using the package GiNaC [74,75]. The basis of GPLs in Table 4 is not unique, but sufficient for numerical evaluation. An example of a third order convolution is the following integral