Higgs Boson Production at Hadron Colliders at N3LO in QCD

We present the Higgs boson production cross section at Hadron colliders in the gluon fusion production mode through N3LO in perturbative QCD. Specifically, we work in an effective theory where the top quark is assumed to be infinitely heavy and all other quarks are considered to be massless. Our result is the first exact formula for a partonic hadron collider cross section at N3LO in perturbative QCD. Furthermore, this result represents the first analytic computation of a hadron collider cross section involving elliptic integrals. We derive numerical predictions for the Higgs boson cross section at the LHC. Previously this result was approximated by an expansion of the cross section around the production threshold of the Higgs boson and we compare our findings. Finally, we study the impact of our new result on the state of the art prediction for the Higgs boson cross section at the LHC.


Introduction
With the discovery of the Higgs boson [1,2] at the Large Hadron Collider (LHC) at CERN we have entered a new era of particle physics phenomenology. With conclusive evidence for the existence of the Higgs boson the Standard Model (SM) of particle physics has become a self consistent theory. It explains the mechanism of electro-weak symmetry breaking, the origin of elementary particle masses and it allows to derive concise predictions to energies far beyond current experimental reach. The SM is however a phenomenologically incomplete theory and needs to be extended to obtain a satisfying description of all known physics. Higgs boson measurements will provide a unique window to deepen our understanding of fundamental interactions and to stringently test possible extensions of our current knowledge.
The inclusive cross section for the production of a Higgs boson represents a prototypical example of experimental and theoretical synergy. Its role in the extraction of fundamental coupling constants is key and it provides an invaluable tool to discover potential deviations from the SM. Experimentally it can be determined at the LHC to astounding precision. In order to exploit the full potential of LHC phenomenology experimental precision must be matched by equally precise theoretical prediction.
The dominant production mechanism of a Higgs boson at the LHC is gluon fusion. In comparison with other processes perturbative QCD corrections to the gluon fusion cross section are large.
In order to match current and future experimental precision this simple fact demands computation of this process to very high order in perturbation theory. Next-to-leading oder (NLO) [3][4][5] corrections to this process are available since more than two decades. Corrections at next-to-next-to leading order (NNLO) were computed in refs [6][7][8] in an effective theory (EFT) of QCD where the top quark is considered to have infinite mass and all other quarks are massless [9][10][11][12]. In ref. [13] next-to-next-to-next-to leading order (N 3 LO) corrections were computed in terms of an expansion around the production threshold of the Higgs boson. This result marked the first computation of a Hadron collider observable to this order in perturbation theory. At the desired level of precision the inclusion of many sub-dominant effects, such as electro-weak corrections and quark mass effects, in a prediction for the hadron collider observable are essential. Furthermore, a critical assessment of all sources of uncertainties is required. A comprehensive study achieving this goal was performed in ref. [14] and resulted in the state of the art prediction for LHC measurements (see also refs. [15,16]).
In this article we go beyond the previous approximation of the N 3 LO corrections to the Higgs boson gluon fusion cross section in the EFT in terms of a threshold expansion and compute it exactly. Our calculation strongly relies on various ingredients already entering the computation of ref. [13]. Specifically, we require matrix elements integrated over phase space for the production of the Higgs boson in association with up to three partons and involving up to three loops. Purely virtual corrections were computed in ref. [17,18]. Contributions with one parton in the final state and two loops were calculated in refs. [19][20][21][22][23]. Matrix elements involving two final state partons and one loop (RRV) or tree level matrix elements with three final state partons (RRR) were computed for the purposes of refs. [13,[24][25][26] in terms of a threshold expansion. Furthermore, our result relies on infrared subtraction terms formed out of convolutions [27,28] of splitting functions [29,30] and an ultraviolet counter term based on lower loop amplitudes [31]. Both were already computed for the purpose of ref. [13].
In order to obtain our result we compute N 3 LO corrections to the partonic cross section due to RRV and RRR matrix elements. The integration over the loop and final state momenta involves complicated, high-dimensional integrals. In order to facilitate our computation we employ the framework of reverse unitarity [6,[32][33][34][35] that allows to relate phase space integrals to cuts of loop integrals. Subsequently, we employ powerful loop integration techniques to actually compute our phase space integrals. In particular, we make use of integration-by-part identities [36,37] in order to express our integrated matrix elements in terms of a limited set of master integrals. We then proceed to compute these master integrals using the framework of differential equations [38][39][40]. The solution of differential equations requires the calculation of one boundary condition per master integral. To obtain these boundary conditions we perform an expansion of every master integral in terms of a threshold expansion. We then match the individual terms in the expansion to so-called soft master integrals that were explicitly computed in refs. [41,42].
When solving differential equations for RRR master integrals we encounter an obstruction in the form of 2×2 systems of differential equations that cannot be solved by conventional means. The solution to these systems is given in terms of elliptic integrals. The appearance of elliptic integrals in the computation of Feynman integrals is well established [43][44][45][46][47][48][49][50] but still poses a considerable challenge. The majority of known analytic results for Feynman integrals can be expressed in terms of iterated integrals referred to as generalised poly logarithms [51]. A profound understanding of their analytic properties [51][52][53][54][55] has been key to the success of higher order perturbation theory. The quest for a similar understanding of iterated integrals involving elliptic functions is subject of ongoing research and has already produced vast literature [56][57][58][59][60][61][62][63][64][65][66][67][68][69][70][71][72]. In particular, methods to find solutions for differential equations, the understanding of functional relations among such integrals and their analytic continuation from one kinematic regime to another are of importance. In this article we present our pragmatic solution to the problem at hand and produce for the first time a hadron collider cross section that involves the analytic treatment of an elliptic integral.
Having obtained analytic results for all required matrix elements with different parton multiplicity in the final states we combine them to form the exact correction to the partonic Higgs boson production cross section at N 3 LO. We then convolute our newly obtained result and all required lower order cross section with parton distribution function in order to derive physical predictions for hadron collider cross sections. We study in detail the deviations of our results from the previous approximation of the N 3 LO cross section [13,14]. Our computation allows us to remove one source of uncertainty due to the truncation of the threshold expansion from the state of the art prediction for the Higgs boson production cross section [14] and we update the previous result.
This article is structured as follows: In section 2 we setup the notation for our computation of the inclusive Higgs boson production cross section. Next, we discuss in detail the analytic computation of the missing RRV and RRR coefficient functions in section 3. We outline the general computational framework in section 3.1. We discuss the treatment of elliptic integrals found when solving differential equations in section 3.2. In section 3.3 we introduce a class of iterated integrals that serve as the main building blocks for our final result. Next, we describe the structure of our analytic results in section 3.4. In section 4 we present numerical results for the EFT Higgs boson cross section through N 3 LO in QCD perturbation theory. We compare our new results to previous predictions obtained with a threshold expansion in section 5. Finally, we draw our conclusions in section 6.

Set-Up
In this article we consider scattering processes of two protons that produce at least a Higgs boson.
(2.1) P 1 and P 2 are the momenta of the colliding protons and p h the momentum of the Higgs boson. The master formula for the inclusive Higgs boson production cross section is 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 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 . We define The sum over i and j ranges over all contributing partons. Furthermore, we define the variablē z = 1 − z. The partonic Higgs cross section is given byσ ij (z, m 2 h ). In this article we compute the partonic cross section through N 3 LO in perturbative QCD in an effective theory where the top quark is infinitely heavy and has been integrated out [9][10][11][12]. In this theory the Higgs boson is coupled directly to gluons via an effective operator of dimension five [73][74][75][76], where H is the Higgs field, G a µν is the gluon field strength tensor and L SM,5 denotes the SM Lagrangian with n f = 5 massless quark flavours. The Wilson coefficient C 0 is obtained by matching the effective theory to the full SM in the limit where the top quark is infinitely heavy.
Within the effective theory, we can write the partonic cross section as Dividing out the Born cross section,σ we can write the bare partonic coefficient functions as, (2.7) The initial state dependent prefactors N ij are given by Here, g, q,q and Q indicate that the initial state parton is a gluon, quark, anti-quark or quark of different flavour than q 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 (n) ij→H+m is the coefficient of α n S in the coupling constant expansion of the modulus squared of all amplitudes for partons i and j producing a final state Higgs boson and m partons summed over polarizations and colors. To compute the n th order partonic coefficient functions we require all combinations l-loop matrix elements with m external particles such that m + l = n.
The occurring loop amplitudes are plagued by ultraviolet divergencies which we regulate using dimensional regularisation and work in d = 4−2ǫ space-time dimensions. We renormalise the Wilson coefficient and strong coupling constant in the MS scheme. Squared matrix elements with fixed parton multiplicity in the final state are separately infrared divergent. These infrared divergences are canceled by summing over all contributing squared matrix elements and performing a suitable redefinition of the parton distribution functions. The resulting partonic cross section is free of divergencies and we refer to the corresponding partonic coefficient function as η ij (z). Various definitions regarding renormalisation and mass factorisation can be found in appendix B. The cross section, eq. (2.2), can be written in terms of finite partonic coefficient functions and physical parton distribution functions f R i as The partonic coefficient functions can be split into two contributions The term η (n), SV ij (z) is comprised of distributions that act on parton distribution functions. The super-script SV signifies that this term represents so-called soft-virtual contributions that arise from kinematic configurations where any parton produced in conjunction with Higgs boson is soft. The coefficient η (3), SV ij (z) was computed in ref. [25] and confirmed by ref. [26]. The coefficient functions η (3), reg. ij (z) represent the so-called regular contributions. Their functional form was approximated with a power series in 1 − z in refs. [13,24,35]. The main result of this article is the complete computation of the coefficient functions η (3), reg. ij (z). We supply this result in a machine readable format in an ancillary file together with the arXiv submission of this article.

Calculation of Coefficient Functions
In order to obtain the partonic coefficient functions η (3) ij (z) we require contributions arising from matrix elements with up to three loops (l ≤ 3) and up to three partons (m ≤ 3) in the final state such that 3 = l + m. The purely virtual matrix elements were computed in refs. [17,18]. Matrix elements with two loops and one emission were computed in refs. [19][20][21][22][23]. Matrix elements with two real emissions and one loop (RRV) and three real emissions (RRR) are so-far publicly only available in terms of the first two expansion terms in the expansion around the production threshold of the Higgs bosons [41,42]. In this article we complete the computation of the N 3 LO coefficient functions. We start by outlining the strategy involved in this computation. Next, we explain the treatment of an ellitpic integral that is part of the RRR coefficient functions. We introduce a class of iterated integrals that serve as building blocks of our partonic coefficient functions. Finally, we obtain the N 3 LO coefficient functions and describe their structure.

Computation of Matrix Elements
In order to obtain RRV and RRR coefficient functions we start by generating all required Feynman diagrams with QGRAF [77]. Next, we perform spinor and colour algebra in a private c++ code based on GiNaC [78]. With this we obtain the loop and phase space integrand for our partonic coefficient functions.
Next, we want to perform the inclusive integral of our integrands over all loop momenta and final state parton momenta. The phase space measure for producing a Higgs boson and m partons is given by All final state momenta are chosen in-going such that the energy component in the above equation appears with a minus sign. In order to perform the loop and phase space integration we employ the framework of reverse unitarity [6,[32][33][34][35] that allows to treat phase space and loop integrals on equal footing. In particular, we represent the on-shell constraints in terms of cut propagators.
The subscript c serves as a reminder that this propagator is cut. Cut propagators can be differentiated just like normal propagators.
d dx They satisfy the condition We can now apply integration-by-part (IBP) identities [36,37] on our combined loop and phasespace integrands. A private c++ implementation of the Laporta algorithm [79] allows us to express our partonic coefficient functions in terms of a limited set of master integrals. To compute these master integrals we work with the method of differential equations [38][39][40]. This method allows to derive a system of partial differential equations for a vector of our master integrals I(z) of the form ∂ ∂z I(z) = A(z, ǫ) I(z). Here, I(z) is a vector of n master integrals and A(z, ǫ) is a n × n matrix with ratios of polynomials in z and ǫ as entries. In order to have a complete system of differential equations we define 550 and 362 master integrals for RRR and RRV respectively. The commonly used strategy to solve such differential equations is to find a n×n transformation matrix T such that Here, A ′ (z, ǫ) is holomorphic in ǫ as ǫ → 0. Having obtained such a form the solution for our master integrals can be easily expressed in terms of a Laurent series in ǫ by Here, I ′ 0 represents a vector of boundary conditions that has to be determined by other means. For the RRV and RRR master integrals such a boundary condition is easily obtained by matching the full solution obtained in eq. 3.8 to an expansion of the required integrals I(z) around the point z = 1 that can be performed by means presented in ref. [41,42].
The art in solving differential equations rests in finding an adequate transformation matrix T . For certain differential equations in a single parameter an algorithmic solution exists [80][81][82][83] and was nicely formulated in ref. [83]. This method applies if a transformation matrix can be found that is comprised of ratios of polynomials in the parameters z and ǫ. For a large subset of integrals in our vector of master integrals I such transformations can be found and we rely on a private implementation of the algorithm outlined in ref. [83] to do so.
For another large class of master integrals it is necessary to find a transformation matrix that contains square roots of polynomials of our parameter z. For these cases we can find the desired transformation by finding suitable algebraic variable transformations that rationalises the square roots involved. Once the roots are rationalised we can again employ the aforementioned algorithm. We point out that this procedure is not particularly algorithmic but leads to a desired solution fairly easily.
We encounter a further obstruction when solving differential equations for the system of RRR master integrals. This obstruction involves the presence of elliptic integrals and we elaborate on our solution in the following section. When solving differential equations for master integrals contributing to the triple real coefficient functions of Higgs boson production at N 3 LO we encounter two coupled 4 × 4 systems of differential equations that we could not decouple order by order in the dimensional regulator by conventional means. In this section we discuss the differential equations in question and present our solution.

An Elliptic Integral in Higgs Production
In figure 1 we display two scalar phase space integrals. Let us choose four master integrals with the same propagators as the scalar integral in figure 1b.
We choose . (3.10) These four integrals satisfy a system of differential equations of the form The vector y(z) represents the inhomogeneous part of the differential equation. The matrix A 1 (z, ǫ) in the homogeneous part of the differential equation is holomorphic in ǫ as ǫ → 0. The homogeneous part of the differential equation that does not decouple as ǫ → 0 is given by the matrix As we can see, for ǫ = 0 two of the master integrals decouple and we are left with a 2 × 2 system for the homogeneous solution of the differential equation.
∂ ∂z In order to decouple our original system of eq. (3.11) we want to find a transformation matrix T E such that We show in appendix A that the functions t ij (z) can be written in terms of complete elliptic integrals and pre-factors. However, this solution is quite unwieldy and we choose another approach here. For all practical purposes it is sufficient to simply define the functions t ij (z) to be the solution to the differential equation eq. (3.14). The homogeneous differential equations for the master integrals E ′ 1 and E ′ 4 , defined by are decoupled as we send ǫ → 0. The inhomogeneity can then be decoupled order by order by in ǫ by standard techniques. A general solution for the differential equations can subsequently be found as illustrated by eq. (3.8).
The second set of master integrals that have the same propagators as the scalar integral depicted in figure 1a can be chosen in such a way that the homogeneous part of their differential equations takes identically the same form as the one already discussed. Therefor we can apply the same transformation matrix to decouple the system order by order in ǫ. With this we found a transformation matrix T that allows us to express the differential equations for all master integrals required for RRV and RRR contributions to Higgs production at N 3 LO in the desired form, eq. (3.7).
In order to derive numerical results for the functions t ij we can solve the differential equations eq. (3.14) in terms of a generalised power series ansatz using the Frobenius method. Consider for example an ansatz for the solution of the system of differential equations as an expansion around z = 0 and z = 1. (3.16) We derived the required structure of our ansatz by regarding the asymptotic solution of the differential equations around the considered expansion points.
Here, t 0 ij and t 1 ij are some numerical boundary constants. Inserting the ansaetze into the system of differential equations we find the following recurrence relations. b (n+2) 11 By comparison to the asymptotic solution given in eq. (3.17) we find all starting conditions for the solution to the recurrence relations. Specifically, we find the conditions b Any choice of boundary conditions will lead to a transformation matrix that satisfies the differential equations eq. (3.14). The only restriction we impose is that the transformation has to be invertible, i.e. that det(T E ) = 0. In accordance with this criterium we make the following choice for the asymptotic solution of the differential equation: t 1 11 = t 1 22 = 0 and t 1 12 = t 1 21 = 1. We find that with this choice the determinant of the transformation matrix is given to all orders in z by Fixing the asymptotic behaviour of the functions t ij (z) in one limit automatically determines their behaviour at any other point. Computing the explicit values for t 0 ij explicitly given the choice we made for t 1 ij is however a non-trivial task. At this point it is useful to reflect on the practical aim of our computation. We desire a solution that is numerically sufficiently precise to determine the complete N 3 LO coefficient functions for values of z in the interval [0, 1] as required for cross section predictions. In this light our solution for the t ij (z) should allow for the desired precision and should be improvable if necessary. This can be achieved by computing an approximation based on a truncated power series.
The regular singular points of our 2 × 2 system of differential equations (3.14) are located at the following values.
Consequently, the power series of the functions t ij (z) around the point z = 1 has a radius of convergence r 1 = 1. Similarly, the power series around the point z = 0 has radius of convergence The domains of convergence for the two power series overlap on the interval z ∈ (0, | 1 2 11 − 5 √ 5 |). In order to determine the boundary constants t 0 ij in terms of the t 1 ij we first compute the truncated power series around both limits under consideration for each t ij (z). Next, we evaluate both series for each t ij (z) at a point within the interval z ∈ (0, | 1 2 11 − 5 √ 5 |). Equating the results allows us to establish a relation among the constants t 0 ij and t 1 ij up to a small, numerical remainder. The remainder can be systematically improved upon by increasing the truncation order of the power series.
Let us briefly introduce a simple method of estimating the size of the remainder of the truncated series. Suppose a function f (x) is given by the convergent series If we truncate the series before order N its remainder would be given by Suppose that asymptotically the ratio of to consecutive coefficients remains constant.
Under this assumption we can estimate the modulus of the remainder to be bounded by Note, that the series converges for |r N x| < 1.
In order to obtain sufficiently high precision for our coefficient functions we perform an expansion of the functions t ij around the expansion points z = 0, z = 1 and z = 1 2 . For each expansion we compute several hundred terms and match the boundary conditions within the overlaps of the domains of convergence. Estimating the remainder of the power series expansion at our matching points suggests that we can easily determine the boundary values with a relative accuracy of 10 −42 or better if needed. In addition to estimating the remainder as described above we evaluate the different power series for the same t ij for several points in the intervals where all series converge and only observe relative deviations at levels smaller than 10 −42 .
In order to further study the convergence of our power series approximation we may regard the asymptotic behaviour of the recurrence relations given in eq. (3.18) and eq. (3.19) as n → ∞. (3.26) We see that asymptotically b Here, the c i are some numerical constants. The numbers 9 22 ± 5 is larger than one. From this we again draw the conclusion that the power series around the expansion point z = 1 is convergent everywhere within the unit interval. The power series around z = 0 is convergent if z < 1/ − 11 2 − 5 . This asymptotic analysis also supports the validity of the procedure to estimate the remainder of the power series truncated at order N defined in eq. (3.25).

Iterated Integrals
In this section we briefly introduce a class of iterated integrals [84] that is particularly convenient to express the solution of differential equations as in eq. (3.8). We define with J(z) = 1. We refer to one ω i (z) as a letter and to an ordered set of letters, {ω n (z), . . . , ω 1 (z)} that defines an iterated integral as a word. Many well known classes of iterated integrals, such as harmonic poly logarithms (HPLs) [85] or generalised poly logarithms (GPLs) [51], that are widely used in particle physics, are sub-classes of this type of iterated integrals. For example the GPLs are given by G(a n , . . . , a 1 , z) = J 1 z − a n , . . . , The presence of the elliptic integrals t ij (z) in the solution of our differential equations does not allow for a solution purely in terms of GPLs. For this reason it becomes necessary to define an extension of GPLs in this article. Already several generalisations of GPLs to accommodate elliptic functions exist in the literature (see for example [60,71,72,[86][87][88]). In the following we review several properties of iterated integrals (see for example [53][54][55]89]). Iterated integrals form a so-called shuffle algebra. The last line of the above equation is then regulated as discussed above. If several right-most letters have poles at the lower end point of the integration we simply iterate the above procedure. We want to be able to rewrite an elliptic integral with argument z in terms of iterated integrals with argumentz = 1 − z or w = 1 2 − z. Let us illustrate how this can be achieved by regarding a transformation from z toz. (3.37) The last line in the above equation is a numerical constant. In order to write the integral in the penultimate line in terms of an iterated integral with upper integration boundz we have to first rewrite the iterated integral in the integrand with an upper integration boundz ′ . To do this we simply apply the above equation iteratively to the integrand. Notice, that the above procedure may be ill defined if the integrand we are considering is divergent at any of the end points. This case is easily avoided by shuffle regulating both end points prior to applying eq. (3.37). Let us demonstrate this step with a well known example. Consider the iterated integral In the above equation we employed a shuffle identity such that right most letter of the function is regular at the new lower integration point z = 1 and that the left most letter is regular at the new end point z = 0. We now may write Combining the the results of eq. (3.38) and eq. (3.39) we find the famous di-Logarithm identity. (3.40) In this example it was possible to determine the integration constant to be π 2 6 analytically. If this is not possible the constant can also be determined numerically with finite precision by simply evaluating the function under consideration before and after variable transformation numerically for any value of z.
The iterated integral representation of eq. (3.28) allows to easily compute truncated power series expansions for the iterated integrals. For example By proceeding iteratively we can easily compute the power series inz for any iterated integral to arbitrary power. In order to obtain compact expressions for our analytic results it is of importance to be able to derive functional relations among our iterated integrals. One of the big advantages of GPLs is that their functional relations are well studied (see for example [51,[53][54][55]90]). The case of generic iterated integrals is not understood at the same level. In ref. [60] it was outlined how relations among iterated integrals involving elliptic functions can be found using IBP identities. Here, we proceed differently.
First, note that our final analytic result will be a linear combination of iterated integrals and pre-factors a i (z), i a i (z)J( ω i ,z).
can be satisfied for some c i = 0 for arbitrary values ofz. The coefficients a i (z) and corresponding iterated integrals J( ω i ,z) are understood to be identical to those appearing in our final result. In order to determine the unknown coefficients c i we expand eq. (3.43) inz. Every coefficient of every power inz has to vanish separately in order for the equation to be satisfied. This allows us to build a system of equations that is large enough to solve for the unknown coefficients c i . If we find a certain linear combination of iterated integrals and coefficients that cannot be constrained with this procedure we found a relation of functions. Let us illustrate the procedure with a trivial example. Consider the simple shuffle relation and let us pretend we do not know already know the coefficients c i . After expanding inz we find We can now create a system of equations by regarding each coefficient inz separately. Technically, we want to find the kernel of the system of equations. We find that the kernel for our example is spanned by the vector {c 1 , c 1 , −c 1 } T . This means we found the shuffle identity Of course this procedure only guarantees that the so-found relations are satisfied up to the order in z at which we truncate our power series. However, we may convince ourselves that the relations are correct by computing as many higher order terms as are to our liking. A more involved example of such an identity is given by (3.48)

Analytic Solution for Partonic Coefficient Functions
In the previous sections we described how we derive differential equations for all master integrals required for RRR and RRV partonic coefficient functions. Furthermore, we outlined how we find a suitable transformation matrix that transforms the differential equations into the form of eq. (3.7). Once, this form is obtained the solution to the differential equations can be conveniently written as in equation eq. (3.8). Iterated integrals as given in eq. (3.28) are particularly suited to represent this solution. Once we calculated all master integrals and computed all boundary conditions we simply insert the master integrals into our IBP reduced matrix elements and obtain the desired result for the partonic coefficient functions. In this section we describe the structure of our final result for the partonic coefficient functions. The set of all letters, the so-called alphabet, that appear in the iterated integrals that constitute the Higgs boson cross section at N 3 LO is given by Note, that the alphabet required to describe all our master integrals individually contains additional letters that drop out in the final expression.
The partonic coefficient functions are comprised of iterated integrals with up to five letters. Typically we find that there are several thousand different iterated integrals in each partonic coefficient function. Applying the procedure outlined in the previous section to find functional identities among these integrals we find that we can express them in terms of only 365 different iterated integrals that cannot be re-written as GPLs in a straight forward fashion. Out of those 188 have letters containing elliptic integrals t ij . For the remaining ones a representation in terms of GPLs may exist.
Having derived moderately compact expressions for our coefficient functions we want to find a method to evaluate them numerically. The conceptually simplest way to evaluate the iterated integrals is to perform every integral numerically. The fact that all our integrals are real valued and are finite renders this approach straight forward. Integrating 5 dimensional integrals numerically is however not particularly fast if a certain level of precision is desired. As an alternative, we want to represent the entire partonic coefficient functions in terms of power series expansions.
Let us first investigate for which values of z we can perform a convergent series expansion. In order to extract this information we regard all singularities and branch points that occur in our alphabet and the algebraic factors of our coefficient functions. We find that they are located at values of z of Here, we included the regular singular points of the differential equations of the elliptic sector, eq. (3.13). In order to evaluate our functions to high precision within the physical interval, z ∈ [0, 1], we decide to perform a power series expansion around the points z = 1, z = 1 2 and z = 0. The associated radii of convergence are then r 1 = 1, r 1 2 = 1 2 and r 0 = 1 2 11 − 5 √ 5 . To obtain a series expansion around our three different expansion points we perform an expansion of all iterated integrals as outlined in the previous section. As the default upper bound for our iterated integrals is the parameterz the expansion around the point z = 1 can be carried out simply by expanding the iterated integrals at the integrand level and integrating subsequently as demonstrated in eq. (3.41). In order to obtain an expansion around z = 0 and z = 1 2 we first re-express our iterated integrals in terms of iterated integrals with upper integration bound z and 1/2 − z respectively. As outlined in section 3.3 this procedure requires us to determine certain integration constants which we obtain numerically by matching series expansions around different expansion points. To ensure that the numerical error introduced by truncating series expansions is sufficiently small we estimate it as explained in eq. (3.25).
We expand the coefficient of every iterated integral in the partonic coefficient function separately around each of our three expansion points and combine the result with the expansion for the iterated integrals. In order to obtain numerical values for the coefficient functions within the unit interval of z we evaluate the expansion around z = 1 in the interval z ∈ [0.75, 1], the expansion around z = 1 2 in the interval z ∈ 1 13 , 0.75 and the expansion around z = 0 in the interval z ∈ 0, 1 13 . We truncate the expansion around z = 1 at O((1 − z) 50 ), the expansion around z = 1 2 at O ( 1 2 − z) 200 and the expansion around z = 0 at O z 100 . Using the estimator introduced in eq. (3.25) we find that this approximates the coefficient functions at any point in the unit interval to a relative numerical precision of 10 −10 or better. This is supported by evaluating the different expansions for several points within the overlaps of their respective domains of convergence and calculating their difference. The numerical precision may of course be improved arbitrarily by simply including more terms in the respective series expansions.

Results
In the previous sections we calculated analytic results for the partonic coefficient functions η (3) ij (z). Our analytic results agree with the power series around z = 1 for the same functions obtained in refs. [13,14]. The leading behaviour of the coefficient functions as z → 0 was correctly predicted in ref. [91]. The coefficient function η (3) qQ was calculated already in ref. [92] and agrees with our result. We derived a representation of the coefficient functions in terms of power series expansions that is particularly useful for numerical evaluation. In this section we present numerical results for the Higgs boson production cross section through N 3 LO.
Let us start by regarding the functional dependence of our coefficient functions. In figure 2 we display the shape of the regular coefficient functions for each distinct partonic initial state. The quark -gluon and and gluon-gluon initial state coefficient functions behave as ∼ log 5 (1 − z) as we approach the value z = 1. The coefficient functions with two quarks in the initial state are tending towards zero in this limit. The limit z → 0 is characterised by a power divergence and all coefficient functions behave as ∼ log 5 (z) z .
In order to derive physical predictions for hadron collider phenomenology we need to convolute our partonic coefficient functions with parton distribution functions (PDF). Throughout this article we will use the PDF sets PDF4LHC15 [93]. We choose a Higgs boson mass of 125 GeV and a top quark mass of m t (m t ) = 162.7 GeV. We choose a value for the strong coupling constant of α S (m Z = 91.1876 GeV) = 0.118. If not stated otherwise we derive numerical predictions for proton-proton collider with a center of mass energy of 13 TeV. We use a private c++ code to perform the numerical convolutions of PDFs and partonic coefficient functions.
In figure 3 we display the impact of the N 3 LO corrections on the hadronic cross section for different initial states as a function of the perturbative scale µ. The gluon-gluon (red) and quarkgluon (blue) initial state contributions were rescaled by a factor of 10 −2 in order to fit nicely. We observe that the numerical impact of these two channels is clearly dominant over all other initial state configurations. The nominally smallest corrections for each channel can be found in an interval of µ ∈ [40,90] GeV.
In figure 4 we combine the contribution from all partonic coefficient functions and evaluate their contribution to the hadronic cross section including lower orders in perturbation theory as a function of the perturbative scale µ. We show LO, NLO, NNLO and N 3 LO predictions in green, orange, blue and red respectively. We observe that the dependence on the perturbative scale is greatly reduced at N 3 LO compared to lower orders. Furthermore, NNLO and N 3 LO predictions overlap within the interval of µ ∈ m h 4 , m h . To derive a concrete numerical prediction we choose the value of the cross section at µ = m h 2 . We vary the perturbative scale in the interval m h 4 , m h in order to estimate the effect of missing higher order corrections at N 4 LO and beyond. As can be seen from figure 4 this procedure is not conservative enough at leading and next-to-leading order. Regarding the progression of the series from NLO onward we observe convergent behaviour. The nominal size of the corrections is greatly reduced at each successive order. Uncertainty estimates based on scale variation overlap at NNLO and N 3 LO.
Our prediction for the Higgs boson production cross section at the LHC based on a computation in perturbative QCD in the large top quark mass limit through N 3 LO of expansion around z = 1 exploits a kinematic enhancement of the gluon luminosity in the collision of protons for lower values of partonic center of mass energy to yield reliable predictions. The point z = 1 represents the production threshold for a Higgs boson, i.e. the lowest possible amount of energy required to produce a Higgs boson. In ref. [14] seven additional terms in the power series were added. The quality of a threshold expansion for N 3 LO corrections was furthermore studied in refs. [24,25,94]. Having now the complete coefficient functions at our disposal we want to reflect on previous estimates and compare our exact analytical findings to the approximate results.
Using the same set-up as in the previous section to derive numerical predictions we find that the hadronic cross section through N 3 LO in perturbative QCD in the infinite top quark mass limit based on thirty terms in the threshold limit is given by We observe a difference of 0.11 pb with respect to our new prediction, eq. (4.1). The scale variation interval in eq. (4.1) is slightly larger. In ref. [14] it was estimated the effect of missing higher order terms in the threshold expansion are less than 0.18 pb. We now see that this estimate was sufficiently conservative.
In the remainder of this section we want to study the behaviour of N 3 LO corrections as a function of the order where the threshold expansion is truncated. In particular we want to investigate its performance for contributions arising from different partonic initial states. In figure 5 we show the N 3 LO correction due to different initial sate partons based on a threshold expansion (red) as a function of the order at which the expansion is truncated. In blue we also display our new result to all orders in the threshold expansion as a reference. We observe that the first four terms show particularly large changes in the derived prediction. Starting from the fifth term we observe slow asymptotic improvement towards the full result. The nominally largest gluon-gluon  [pb] (c)   and quark-gluon channels are approximated better than their purely quark initiated counter parts. The sum of all channels can be seen in figure 5f. In order to see more clearly the quality of the threshold expansion for each channel we show in figure 6 the impact of N 3 LO corrections on the hadronic cross section due to different partonic initial states. The predictions in red are now based on a threshold expansion normalised to the respective all order result. The x-axis shows the order at which the threshold expansion is truncated. The line in blue at one serves as a reference. We observe that contributions originating from the gluon-gluon channel are approximated within several per-mille including only a few terms in the expansion.

1.005
Truncation Order     Similarly the quark-gluon initiated contributions are approximated reasonably well below a level of ten percent. All other contributions are considerably different from the exact result and receive corrections of the order of 100 % even with thirty terms in the expansion. Their nominal effect on the inclusive cross section is however negligible. The fact that the threshold expansion works best for gluonic initial states can be explained by the fact that the probability to extract gluons from a proton is peaked towards lower momentum fractions, i.e. closer to the production threshold. For quarks this enhancement is not as large. The relatively slow improvement towards the exact result of the predictions as more and more terms in the threshold expansion are included can be understood from the high energy behaviour of the partonic coefficient functions. As we displayed in fig. 2 the coefficient functions have a power like divergence ∼ log 5 (z) z as z → 0. While the threshold expansion is formerly convergent within the entire physical interval a relatively slow convergence to capture the high energy behaviour can be expected.

Conclusions
In this article we present an exact computation of the Higgs boson production cross section at hadron colliders through N 3 LO in perturbative QCD in the infinite top quark mass limit. The main result of this article are analytic formulae for N 3 LO corrections to the regular partonic coefficient functions. We provide these functions in an ancillary file together with the arXiv submission of this article.
To obtain our result we compute matrix elements for the production of a Higgs boson in association with three partons at tree level and with two partons at the one-loop level. In order to perform required phase space integrals we employ the framework of reverse unitarity and make use of loop integration techniques such as IBP identities and master integrals. We compute all required master integrals using the framework of differential equations. When solving differential equations we encounter elliptic integrals in the solution for triple real radiation master integrals. We find that an analytic solution for our master integrals can be easily found by embedding the solution of our differential equations in a fairly general class of iterated integrals. We discuss in detail how we find relations among iterated integrals involving elliptic functions and how we evaluate them efficiently numerically.
Having obtained analytic expressions for all required partonic cross sections we embed them in a numerical code and derive predictions for hadron collider cross sections. We find that N 3 LO corrections are small compared to corrections at previous orders and that the dependence on the perturbative scale is greatly reduced. We perform a detailed comparison with a previous approximation of N 3 LO corrections based on an expansion around the production threshold of the Higgs boson including 37 terms [13,14]. We observe that our new results are in excellent agreement with this approximation. Dominant contributions due to gluon initiated partonic cross sections are approximated rather well by the threshold expansion. Quark initiated contributions on the other hand are approximated rather poorly. The estimate of missing higher orders in the threshold expansion in refs. [13,14] was sufficiently conservative to cover the difference to the exact result.  Table 1: Cross sections and uncertainties as function of the collider center of mass energy.
To derive precise predictions for hadron collider phenomenology many effects beyond the effective theory cross section considered in this article have to be take into account. The finiteness of quark masses and neglected electro-weak effects play an important role. It is particularly important to critically asses all non-negligible sources of uncertainty. A detailed study of the inclusive production cross section for the Higgs boson considering all such effects was conducted in ref. [14]. Repeating this analysis is beyond the scope of this article. However, it easily possible to modify the final predictions for hadron collider cross sections of ref. [14] such that the results of this article are taken into account. Specifically, we include the exact contributions to the cross section at N 3 LO in the EFT and remove uncertainties due to the truncation of the threshold expansion. Otherwise, we can simply use the results of ref. [14] that are neatly combined in a new numerical code iHixs 2 [95]. In table 1 we show updated predictions for the gluon fusion Higgs boson production cross section at the LHC as in ref. [95].

Acknowledgements
I would like to thank Babis Anastasiou, Claude Duhr, Falko Dulat and Franz Herzog for many invaluable discussions. I would like to thank Claude Duhr and Babis Anastasiou for useful comments on the manuscript. I am grateful to Lorenzo Tancredi and Stefan Weinzierl for useful discussions about elliptic integrals appearing in this computation. Furthermore, I would like to thank Achilleas Lazopoulos and Falko Dulat for comparisons of numerical results using the latest version of iHixs 2 [95]. My research was supported by the European Comission through the ERC Consolidator Grant HICCUP (No. 614577).

A The Elliptic Integral
In section 3.2 we discuss a coupled system of two differential equations that describes the homogeneous solution to master integrals appearing in triple real radiation matrix elements when integrated over phase space. The particular system is given by Equivalently, we can say that E 0 4 satisfies a second order differential equation.
First, a solution to this differential equation was found by Stefan Weinzierl in terms of an elliptic integral. The homogeneous part of a differential equation for a Feynman integral has to be satisfied by the maximum cut of the corresponding Feynman integral. In ref. [40] it was proposed that it is sufficient to normalise the leading singularities of Feynman integrals to constants in order to decouple their differential equations order by order in the dimensional regulator. For this to hold true the physical linear combinations of leading singularities themselves must satisfy the homogeneous differential equation for ǫ = 0. Computing the leading singularity of E 4 we find We can rewrite the quartic polynomial under the square root as Following the prescription of ref. [57] we define two integrals (r 4 − r 2 )(r 3 − r 1 ) K(1 − m).
Here, K(m) is the complete elliptic integral of the first kind. We find that both integrals I 1 and I 2 are solutions to our second order differential equation eq. A.2. In principle we could now follow a procedure outlined in ref. [57] to construct a transformation matrix T E that allows us to decouple the system of differential equations order by order in ǫ. Specifically, we find that the functions t ij (z) defined in section 3.2 are given by linear combinations t ij (z) = c 1 I 1 + c 2 I 2 + c 3 z ∂ ∂z I 1 + c 4 z ∂ ∂z I 2 , c i ∈ C. (A.5) The derivatives of the functions I 1 and I 2 with respect to z yield a sum of elliptic integrals of first and second kind with algebraic pre-factors. We can determine the coefficients c i analytically by equating the power series expansions of the above equation with the results obtained in section 3.2. However, any of these analytic expressions is quite unwieldy.

B Various Ingredients for Higgs Boson Production
In this appendix we summarise various standard ingredients for the perturbative calculation of the inclusive Higgs boson production cross section. In order to perform renormalisation in the MS scheme we substitute the bare coupling and Wilson coefficient as The renormalisation factors for the strong coupling constant and Wilson coefficient required for a computation through N 3 LO [18] are given by (B. 2) The coefficients at the various orders in the coupling constant β i are given by the QCD beta function [96][97][98][99].
In order to obtain infrared finite cross sections we are required to perform a suitable redefinition of our parton distribution functions.
The infrared counter term Γ consists of convolutions [28] of splitting functions P (n) ij [29,30] and can be derived from the DGLAP equation. Its perturbative expansion required for an N 3 LO accurate calculation of the differential Higgs boson production cross section is given by