On the computation of finite bottom-quark mass effects in Higgs boson production

We present analytic results for the partonic cross-sections contributing to the top-bottom interference in Higgs production via gluon fusion at hadron colliders at NLO accuracy in QCD. We develop a method of expansion in small bottom-mass for master integrals and combine it with the usual infinite top-mass effective theory. Our method of expansion admits a simple algorithmic description and can be easily generalized to any small parameter. These results for the integrated cross-sections will be needed in the computation of the renormalization counter-terms entering the computation of finite bottom-quark mass effects at NNLO.


Introduction
There have been recent improvements in the determination of the inclusive Higgs boson cross section with the computation of the three-loop corrections to the production of a Higgs boson via gluon fusion in the infinite top-mass limit [1]. The uncertainty due to scale variation has been shown to be less than 3% in this case. With such a level of precision arises the question of the uncertainty coming from other contributions. In particular, corrections coming from a finite bottom-quark mass account for around 6% of the total Higgs boson production cross section at leading order and around 4% at NLO [2,3]. This means that one can expect the uncertainty coming from the unknown finite bottom-quark mass effects at N 2 LO to contribute substantially to the total uncertainty on the inclusive Higgs cross section at the LHC and the computation of such effects is desirable.
The dominant production mechanism of the Higgs boson within the Standard Model at the LHC is gluon fusion through a heavy-quark loop. Because of the large value of the Yukawa coupling y t of the top quark, the corresponding cross-section is dominated by contributions involving a top-quark loop [4]. Since the Higgs boson is light compared to the top quark, it is justified to work in the limit of an infinite top quark mass. This induces an effective interaction of the Higgs boson with gluons described an effective scalar operator where c H is the Wilson coefficient of the effective interaction, F a µν is the gluon field-strength tensor, H is the physical Higgs field, and v is the vacuum expectation value of the Higgs field. The Wilson coefficient can be found in the literature, see for example [5][6][7].
Despite its much smaller mass, the bottom-quark also contributes substantially to the inclusive Higgs cross-section. In particular, the top-bottom interference, i.e. contributions where diagrams containing a top-quark loop are interfered with diagrams containing a bottom-quark loop, has to be taken into account. These contributions are suppressed by a factor of y b /y t = m b /m t 0.025 compared to the ones with two effective vertices but become important well above the bottom-quark threshold, where m 2 b s ∼ m 2 H , due to large logarithms. Note that contributions where two effective operators are replaced by b-quark loops are further suppressed by a factor of y b /y t and will not be considered here. These large logarithms can be resummed to all orders in the strong coupling constant and have been found to have only a mild impact on exclusive observables [8]. Here we concentrate on the fixed-order inclusive cross-section with complete m b dependence whose computation presents two main difficulties: Firstly because the bottom-quark mass is much smaller than the physical Higgs mass it is impossible to simplify the calculations using an infinite-mass effective theory for these contributions; secondly the integrals required at N 2 LO with full dependence on m 2 b are unknown in the literature and very difficult to compute because of the many mass scales involved. We suggest that a more amenable strategy is to compute these finite bottom-mass effects as an expansion in small m 2 b . We explore the difficulties associated with such a double expansion -large top mass and small bottom mass -and obtain analytic expressions for the first order of the small m 2 b expansion of the top-bottom interference contributions to the inclusive Higgs cross section via gluon fusion at NLO QCD. This shall serve both as a proof for the feasibility of such an approach and as a preparation for the computation of the corresponding N 2 LO QCD corrections. In particular we develop tools and understanding that will be valuable for the N 2 LO computation: (i) we introduce a systematic method of expansion of Feynman integrals based on differential equations and (ii) we explore how the phase-space integration create non-trivial analytic behaviour of the inclusive integrals as m b → 0 and we were able to obtain the correct small m 2 b expansion. Note that the two-loop amplitudes needed at NLO along with the corresponding master integrals are already known in the literature with full m b dependence [4,[9][10][11][12] and that the real radiation matrix elements have been obtained in references [13][14][15]. While the known two-loop virtual master integrals will merely serve as a check of our method of expansion, we present analytic results for the cut two-loop real master integrals for the first time. We also present analytic expressions for the contribution of the top-bottom interference to the cross section for the production of the Higgs boson via gluon fusion at NLO at order . This will serve for the computation of counter terms for the renormalization at NNLO.
This article is organized as follows. In Section 2 we present our notation and describe the calculation. In Section 3 we present our analytic results for the first term in the small m 2 b expansion of the Higgs production cross section via gluon fusion at NLO QCD. In Section 4 we describe our method of expansion of the master integrals; and in section 5 we present the analytic results for the master integrals appearing in the virtual and real radiation contributions to the cross section. In Section 6 we give analytic results for the squared matrix elements in terms of master integrals.

Description of the calculation
We consider the hadronic production of a Higgs boson, denoted symbolically as P (P 1 ) + P (P 2 ) → H(p H ) + R(p 3 , . . .), where P denotes a proton and R is any additional QCD radiation. Let us denote by S = 2P 1 · P 2 the center-of-mass frame energy of the two colliding protons. Assuming the usual factorization, the total hadronic cross section can be written as where f i is the parton density function for partons of type i. The partonic cross section for the process i(p 1 ) + j(p 2 ) → R(p 3 , . . .) + H(p H ) is given by where p 1 = x 1 P 1 and p 2 = x 2 P 2 are the momenta of the incoming partons, s = 2p 1 · p 2 denotes the center-of-mass frame energy of the colliding partons and |M ij→RH | 2 is the corresponding squared matrix element. We can rewrite the production phase-space as dΦ(p 1 , p 2 ; p 3 , . . . , p H ) = dΦ(p 1 , p 2 ; p 3 , . . . , p H ) s dz δ(zs − m 2 H ), where the four vector p H has the mass-shell condition p 2 H = sz on the right-hand side and we have z = m 2 H /s ∈ [0, 1]. We can then use the mass-shell delta function δ(zs − m 2 H ) to put constraint on s rather than on z and use the latter as an integration variable parametrizing the soft limit of the integral. This amounts to rewriting the total hadronic cross section (2.1) as where τ = m 2 H /S is the production threshold, σ ij→RH (z) is the partonic cross section with p 2 H = zs, and we defined L ij (z) as This representation allows one to consider the inclusive cross-section as a distribution with respect to the integration variable z and allows the implementation of the soft subtraction without any explicit reference to the luminosity. Note that such a representation requires that we choose z and m 2 H as independent variables such that we must set s = m 2 H /z. We can now define more precisely the different pieces of this computation. The partonic cross-section can be decomposed as where 'int.' stands for interference and 'eff.' for effective theory. The cross section σ eff. ij→H is the cross section for the partonic process ij → RH in the infinite top-mass effective theory alone (with y b = 0 and the effective interaction given by (1.1)) and starts at order y 2 t ∝ m 2 t . Such contributions can be found in the literature up to α 5 s [16]. The cross-section σ int.
gg→H is the top-bottom interference and contains the first finite bottom-mass effects. It can be further decomposed in Born, virtual, and real contributions as where n ij is the averaging factor and we made the sum over external polarizations and spins implicit. We denote by A (n) ij→RH the n-th order QCD correction to the matrix element for the production of RH mediated by a b-quark loop (with y t = 0), and by B (n) ij→RH the nth order QCD correction to the matrix element for the production of RH in the infinite top-mass effective theory (with y b = 0). Note that there is no single diagram involving both a top and a bottom loop at this order and that the Born matrix element A (0) ij→RH contains one b-quark loop. Typical diagrams contributing to these cross sections are shown in figure 1.
The Born and virtual contributions are known analytically in the literature [12]. The integrated real contribution has only been obtained numerically and we present here its first fully analytic computation as an expansion in m b . The different channels contributing to σ int.
ij→H are: • Born and virtual: Only the gg-initiated channels contributes to the Born and virtual cross sections, such that we will only consider σ B;int. gg→H and σ V ;int. gg→H .
• Real: There are new channels contributing to the real cross sections, namely the intial states gg, qg,qg, gq, gq, qq, andqq. All these contributions can be easily obtained from σ R;int. gg→gH , σ R;int. qg→qH , and σ R;int. qq→gH , and only these will be presented here.
Our strategy for the computation of the inclusive cross-section is based on a reduction to master integrals using integration-by-parts (IBP) identities [17][18][19][20] and on the computation of these masters integrals by the method of differential equations [19,[21][22][23]. The Feynman diagrams are generated using the program FeynArts [24] and we use the program FeynCalc [25] to evaluate the color and spin traces. We use our own Mathematica code in order to compute the square of the matrix elements. The calculation is performed in Feynman gauge and we therefore add contributions with external Faddeev-Popov ghosts. The reduction to master integrals has been performed in two independent ways -using both the program FIRE [26,27] and a private C++ code 1 -and found to agree. Differential equations with respect to the b-quark mass are trivial to obtain because m b enters the integrals only in internal propagators. The differential equations are then solved as an expansion with respect to m 2 b using the method presented in section 4. Note that contributions from the b-quark loop vanish for m b = 0 such that the cross-sections σ int.
The inclusive real contributions can be computed in a similar way thank to the method of reverse unitarity [28][29][30][31][32]. This method is based on the usual relationship between cuts and discontinuities of Feynman propagators [33], which is given by In the reverse unitarity approach one takes this identity for granted and treat the phasespace integrals as loop integrals, i.e. the phase-space integrals are subject to the same IBP identities as their loop counterparts. The phase-space integrals can then be reduced to master integrals and differential equations can be obtained in the usual way. The only difference with the normal reduction is that integrals that have vanishing or negative power for any of the cut propagators must vanish. Note that powers bigger than one are allowed and correspond to derivatives of the delta function in equation (2.4), in the sense of distributions. The method has been successfully applied to a variety of cut integrals [29][30][31][32]34] and we refer the reader to these references for further discussion.

Results
In this section we present analytic results for the first order in the small-m 2 b expansion of all the relevant partonic cross sections contributing to σ int.
ij→H up to order α 3 s , as described in the previous section. Our method of expansion allows one to generate arbitrary high orders in the m 2 b expansion but we will display here only the leading terms and higher-order terms can be easily generated using the method presented in section 4. Since we are not performing any detailed phenomenological study here, we will only give expressions for the unrenormalized partonic cross sections, while singling out contributions cancelled by the different renormalization counter terms. We refer the reader to references [12] for a more detailed discussion.
We start by presenting our results for the gg-initiated channel. The inclusive Born and virtual contributions decouple completely from the phase-space integration, which contributes as an overall pre-factor The exact result for these contributions is known analytically [12] with full dependence in m 2 b and can be easily expanded 2 for small values of the bottom quark mass. We performed an independent computation and found agreement. This provides a strong check of both our method of expansion and our evaluation of the boundary conditions. The leading order of the small bottom mass expansion of the Born contribution reads where we defined where S = exp { (log 4π − γ E )} and µ is the 't Hooft scale. Note that because of the δ-function in equation (3.1), we must have that L = log m 2 b /m 2 H here. In contrast, the real contributions do depend explicitly on z as the emission of a supplementary final-state parton allows the initial-state partons to have a center-of-mass frame energy different from the Higgs boson mass m 2 H . Since the real matrix element becomes singular in the soft limit z → 1, where the energy of the supplementary finalstate parton vanishes, it is necessary to implement some form of soft subtraction. We find that the expansion in small m 2 b completely preserves the structure of the usual NLO soft subtraction 3 and we find no special difficulty there. Namely, we write whereσ R; int. gg→H = lim z→1 σ R; int. gg→H (z)/(1 − z) −1−2 is the soft limit of the real contributions and has the expected form [37]. The term in parenthesis in (3.3) is regular as z → 1 and can be safely expanded in while the second term is expanded in plus-distributions using the identity where the plus-distributions D n are defined as for an arbitrary function f . After the soft subtraction has been performed, the total next-to-leading order contribution can be written as where β 0 is the first order of the QCD beta function when expanded with respect to α s /4π, that is the g → gg splitting kernel is given by and the convolution is defined for arbitrary functions f and g as To get the correct form of the PDF counter term it is important to remember that σ B; int.
gg→H is proportional to δ(1 − z) and that we are using z and m 2 H as independent variables, not s. The first term in (3.4) is cancelled by the α s and mass renormalizations while the second term is cancelled by the introduction of the PDF counter terms. The rest of the expression is what is left after the introduction of all renormalization counter terms and can be written in terms of a δ part given by a plus-distribution part given by and a regular parts given by where Li i and S i,j denote the polylogarithms and the (Nielsen) generalized polylogarithms, respectively.
The virtual contribution alone has the singularity structure predicted by [37] and can be written as The real contribution must be -finite before integration over the phase-space since the internal b-quark mass m 2 b regulates the internal bottom-quark loop completely. This means that real contribution has only a 1/ -pole, which must be proportional to the collinear splitting kernel, . Although we have the presence of β 0 in collinear and renormalization counter terms, the matrix elements do not contain it, since there is no explicit correction to the gluon propagator. In particular, this means that our result is independent of the number of light-quark flavours N f apart in the running of α s .

Systematic expansion of Feynman integrals
In this section we introduce a method for the expansion of Feynman integrals with respect to a small parameter. It is a simple extension and systematisation of the method of expansion developed recently in [1,38,39] and we show here that it can be formulated as an explicit algorithm and that a closed-form formula in the form of a Dyson expansion can be obtained. The procedure presented here is completely general and can be carried out with any family of integrals and any external parameter. It is based on the general method of differential equations [19,21,40] as well as more recent ideas presented in [41] and we refer the reader to references [42,43] for a modern introduction to the subject.
The usual strategy is to seek solutions of the differential equations in the form of an expansion in , see for example [44,45]. Many different types of functions appear in such expansions: most notably polylogarithms, harmonic polylogarithms [46] and Goncharov polylogarithms [47]. These functions are generally embedded with a powerful algebraic structure [48,49] that enables spectacular simplifications [50]. However, for integrals with many different external parameters, more complicated functions appear, such as elliptic functions [51], and no useful algebraic structure generating relations among them has been found so far. Handling such functions properly in a physical computation remains a difficult challenge. It might be easier to obtain solutions of the differential equations as expansions in terms of an external parameter. In such expansions, only powers and logarithms of the expansion parameter appear, hence removing the hassle of dealing with complicated functions 4 . The hope is that single terms in the expansion are simpler to compute than the complete result, and that by summing a sufficient number of these terms it might be possible to obtain a satisfactory approximation.
Let us consider an arbitrary Feynman integral I that depends on the small bottom-mass m 2 b . Extracting a pre-factor of (−s) 2d−ν , where ν = i ν i denotes the total inverse power of propagators, the m b -dependence of this integral can be expressed using the dimensionless variable In what follows, we will set s = −1 and consider only the variable r from now on. In general, Feynman integrals have non-trivial analytic properties when considered as functions of their external parameters and possess branch cuts whose starting points are determined by the Landau equations. If one tries to expand the Feynman integral around a point which is not a solution of the Landau equations, the integral is analytic there and a naive expansion of the integrand will give the correct Taylor expansion. However, if the expansion point r = 0 is a solution of the Landau equations then the integral is not analytic there and can not be naively expanded. In this case the expansion can be performed using the prescription of expansion by regions [52,53] and the general analytic behaviour of the integral in the vicinity of the expansion point r = 0 is found to be of the form where the functions I (n,m) have a Laurent expansion around r = 0. Only some of the I (n,m) are non-vanishing such that the sum over n and m contains only a finite number of terms. Such an expansion makes the analytic behaviour of the Feynman integral in the vicinity of the expansion point explicit since the functions I (n,m) do not posses any branch point at r = 0. We call contributions with different n as having different scalings while the term proportional to I (0,0) has a trivial analytic structure and will be denoted the hard region.
In order to solve a system of differential equations in the form of an expansion, it is necessary to resolve the analytic structure as shown explicitly in (4.2). In particular, one has to understand how the explicit logarithms, known in the literature as integrating factors, appear in such an expansion [19,22,54]. In the rest of this section, we will see that once the analytic structure of the solution of the differential equations in the vicinity of the expansion point has been exposed, one can easily generate an expansion in the form of (4.2) explicitly.
Note that although we have introduced the parameter r as (4.1) for definiteness, the algorithm developed in what follows can be applied to any small parameter and is not restricted to expansions with respect to internal masses. Such a parameter can be the square of any linear combination of the external momenta, any internal mass, or any combination thereof. For example, while computing the virtual corrections in section 5.1 we use the (related) variable x = ( √ 1 + 4r − 1)/( √ 1 + 4r + 1) as our expansion parameter.

Analytic behaviour near a singular point of the differential equations
Let us first comment briefly on the link between the structure of a set of differential equations and the analytic behaviour of its solution in the vicinity of a regular singular point. A similar discussion can be found in reference [54]. The considerations presented here are not new nor very advanced and we refer the reader to textbook literature for more detailed discussion. Consider a family of master integrals, M(r, ) = (I 1 (r, ), . . . , I m (r, )) T , that depends on an external parameter r and further assume that it satisfies a system of differential equations with a regular singular point at r = 0, that is where the dots indicate contributions less singular at r = 0 and the residue matrix R( ) is not nilpotent. As we will discuss later, if the residue matrix is nilpotent the pole can be removed altogether by a suitable change of master integrals and the solution of equation where B( ) denote some vector of boundary conditions. Matrix exponentials are best computed using the Jordan decomposition, that we recall here for completeness and to fix notation. For every square matrix R with eigenvalues λ 1 , . . . , λ n there exists a similarity transformation S[R] such that is of block-diagonal form and can be written as with a i ∈ {λ 1 , . . . , λ n } and d i ∈ N. The matrix J[R] is the Jordan normal form of R and the d i -dimensional square matrices J d i a i are its Jordan blocks. For a given eigenvalue λ i , the number of Jordan blocks with a j = λ i is called the geometric multiplicity of this eigenvalue, while the sum of the sizes of the Jordan blocks is called the algebraic multiplicity. A matrix R is diagonalizable if and only if the geometric and algebraic multiplicities are equal for each eigenvalue, that is, if all the Jordan blocks are 1 × 1 matrices. Using the very definition of the matrix exponential, we obtain that exp is an arbitrary constant, such that only the exponential of J[R] needs to be computed. The exponential acts on each Jordan block independently and the exponential of a single Jordan block can be easily obtained 5 as We are now ready to understand how solution (4.4), which we now write as develops non-trivial analytic behaviour at r = 0. The only non-trivial analytic structure comes from the exponentiation of each Jordan block and appears in two distinct forms: (i) the exponential pre-factor in (4.7), which is of the form r λ i for c = log r, and (ii) the logarithms arising from the exponentiation of non-trivial Jordan blocks. This corresponds exactly to our decomposition (4.2). Hence explicit logarithms appear by the exponentiation of non-trivial Jordan blocks of the residue matrix and the highest power of the logarithms is given by the dimension of the corresponding Jordan block.
In turn, the fact that the general structure of the asymptotic expansion of a Feynman integral must be of the form given by equation (4.2) puts constraints on the system of differential equations satisfied by the master integrals. In particular, it means that • Feynman integrals do not have essential singularities, such that the system of differential equations has a regular singular point at r = 0.
• If the system is in a form where it can be expanded as in equation (4.3), then the residue matrix R( ) must have eigenvalues of the form At first sight, the fact that the system has a regular singular point does not imply that it can be directly expanded as in equation (4.3) but there are known methods to achieve such a task, as we discuss below.

Description of the algorithm
We now present our algorithm in detail. Let us consider a family of master integrals M(r, ) = (M 1 (r, ), M 2 (r, ), . . .) satisfying a linear system of differential equations given by where A is a m × m matrix whose entries are rational functions of r and . Under a linear redefinition of the master integrals M(r, ) → M (r, ) = T(r, )M(r, ) the system of differential equation becomes where the transformed matrix is given by (4.10) The first step of our algorithm is to find a transformation such that the system can be expanded as in equation (4.3) and solved by (4.4) in the vicinity of r = 0. In general, even if the system of differential equations has a regular singular point at r = 0, the Laurent expansion of the matrix A(r, ) around r = 0 can have poles of order higher than one. In this case, the residue matrix of the highest pole must be nilpotent, hinting at a 'fake' singularity. This issue was investigated long ago by Moser [55] who showed that it is always possible to construct a transformation -that we will call T rank in what follows -such that the transformed matrix D T rank x [A] has a singularity of order one as in equation (4.3). When A is rational, such a transformation has further been shown to be rational as well by Barkatou and Pflüegel [56,57], who also provided a practical algorithm 6 for its computation. It should be noted that such a transformation is not unique. This algorithm has been recently used in reference [41] in order to reduce maximally the degree of all the singular points with respect to a given variable simultaneously, whenever it is possible. If this is achieved, the solution of the differential equation can be easily obtained in an -expansion [23,58]. However, the very existence of such a form is a tremendous mathematical problem closely related to Hilbert's 21st problem [41]. Here we completely avoid such considerations because reducing the apparent order of a single singular point is always possible and can be automatized, possibly at the price of worsening the behaviour of the matrix A at other singular points.
The second step of our algorithm is to normalize the eigenvalues of the residue matrix of the transformed system to be proportional to as first suggested in [41]. After performing the transformation T rank , the system must have the following expansion where the dots indicate less singular contributions. The eigenvalues of the residue matrix R 0 ( ) can be normalized by the following standard trick. Let S[R 0 ] be the similarity transformation of R 0 into its Jordan normal form and further denote by u i the i-th column vector of this transformation and by v i the i-th row vector of its inverse, that is These vectors are the generalized eigenvectors of the matrix R 0 ( ). It is possible to shift the i-th eigenvalue of the residue matrix from λ i to λ i + 1 using the transformation In reference [41], such a transformation is called a balance between 0 and ∞. Since all the eigenvalues of R 0 must be of the form λ i = v i − n i with v i ∈ Z and n i ∈ N, we can normalize R 0 by applying the total transformation (4.14) The residue matrix of the transformed system D Ttot.
x [A(x, )] , where T tot. = T rank T norm. , then has eigenvalues proportional to . Note that this is independent from the order by which we construct this transformation, since the matrices u i v T i commute with each other for different i.
The third and final step of our algorithm is to construct the expansion explicitly. After the two first steps have been carried out, we must have by construction that where the eigenvalues of the residue matrix R are all proportional to . Let us further define the vectorM given by where r R( ) = exp {log(r)R( )} is the fundamental matrix of the system for which only the pole is retained. Following our previous discussion, we know that the matrix exponential contains the whole non-trivial analytic structure and we expectM to have a Taylor expansion at r = 0. Indeed by using the transformation rule (4.9), we see thatM satisfies a differential equation without any singularity at r = 0, namely Because we normalized the eigenvalues of the matrix R(r, ), every entry of the matricesÃ i must be proportional to r n with n ∈ Z, such that every term in (4.16) contributes to a single order in r for = 0. This fact is of crucial importance because it allows us to generate terms of the expansion order by order in r. The sought after expansion can be obtained by means of a Dyson series as follows. Turning (4.16) into an integral equation, we get where B( ) is an arbitrary vector of boundary conditions. Substituting this very equation into itself iteratively, we get the expansioñ where the matrix coefficients of the Dyson series are given by D 0 (r, ) = 1, . . . where the matricesÃ i are defined as in (4.16). Note that every entry of D m is proportional to r m+n with n ∈ Z, such that the expansion (4.18) indeed generates a Taylor expansion for each scaling. Note however that the boundary conditions B( ) are common to all orders and have to be obtained from the first non-vanishing order of each master integral. Summarizing, the solution of the differential equation (4.8) can be written as an expansion in the small parameter r by the following steps: • Reduce the singular rank of the matrix A at the point r = 0 to its minimal value using the rational Moser algorithm [56,57]. This produces a change of basis T rank .
• Normalize the eigenvalues of the residue matrix of D T rank r [A] at r = 0 using equation (4.14). This produces a further change of basis T norm. . Note that even if the Dyson series can be obtained order by order, the transformation matrices T rank , T norm and S[A 0 ] are in general not homogeneous in r and will mix the different powers of the expansion. In practice one has to be careful and track which terms of the Dyson series have been truncated.
The computation of the Jordan normal form in r R( ) can be computationally intensive and become untractable for big systems of differential equations. In this case, a simple extension of our algorithm presented in appendix A can be used to take advantage of the fact that most systems of differential equations ecountered in practice are not completely coupled and can be solved iteratively, avoiding the hassle of dealing with large matrices R( ).

Example: Expansion of a family of factorizable two-loop integrals
We illustrate the method of expansion presented in the previous section by considering the following family I[ν 1 , ν 2 , ν 3 , ν 4 ] of factorizable two-loop integrals: This family corresponds to the closed sub-system T V 1 [ν 2 , ν 3 , ν 4 , ν 1 , 0, 0, 0] of the two-loop integrals appearing in the computation of virtual contributions that we will consider in section 5.1. We want to expand these integrals with respect to m 2 b . As explained earlier, we set s = −1 and consider only the variable r = −m 2 b /s from now on. Note that these integrals can be factorized as products of one-loop integrals such that their computation using our method of expansion is somewhat artificial since a direct computation of the one-loop integrals would be easier. However, the very fact that these integrals are easy to compute allows us to illustrate our method pedagically in details.
Differential equations with respect to r for the vector of master integrals M can be easily obtained and read We see that the matrix A has a second order pole at r = 0 and the first step of our algorithm is to find the transformation T rank (r). This can be done either by using the rational Moser algorithm or, in this simple case, by inspection. The second transformation matrix T norm (r) can then be obtained using equation (4.14), and we have Here the transformation matrices have a very simple diagonal form but in general the two first steps of the algorithm can be computationally intensive. Under the transformation T tot. (r) = T rank (r)T norm. (r), the matrix A becomes The matrix R has the following Jordan decomposition and we see that its eigenvalues are properly normalized. We are now ready to perform the expansion. The solution of the differential equation is then given by where r R( ) = exp {log r R( )} can be easily obtained from (4.23), using the explicit formula (4.7), and is given by The first order of the Dyson series can be easily computed as Higher orders can be obtained similarly using (4.19).
we obtain the solution (4.24) explicitly as (4.25) We finally have to determine the boundary conditions B i . The three two-loop integrals considered here are factorizable, i.e. they can be written as the product of two one-loop integrals, for which all order in results are known. In particular they have the following asymptotic behavior where ψ denotes the digamma function and γ E is the Euler-Mascheroni constant. We see that the structure of these asymptotic expansions matches the solution (4.25) exactly at lowest order in r, such that we can easily determine the boundary conditions to be The higher order terms in (4.25) are then completely determined as well.

Master integrals appearing in the virtual contributions
In this section we detail the computation of the two-loop scalar integrals needed in the virtual contributions σ V;int. gg→H . There are 17 master integrals falling into the three distinct topologies shown in figure 2 and defined by Although every topology is closed under IBP reduction, such that the corresponding differential equations could be treated separately, we found easier to work with all of them at the same time, avoiding the duplication of simple master integrals shared by different topologies. We follow reference [12] and choose for our master integrals the following: where we omitted master M V 12 which has a non-trivial numerator and can be written as Up to a pre-factor of (−s) 2d−ν , where ν denotes the total inverse power of propagators, these integrals can be written using the dimensionless variable Since this pre-factor can be easily reconstructed from dimensional analysis, we will set s = −1 and use only the variable x to present our results. Introducing such a variable is necessary when considering the system at arbitrary values of the mass m 2 b but is not mandatory in our case. We use it here to simplify comparison with the known full results. Analytic continuation is unambiguously defined by the causal prescription and we have such that x → x + i0 as well.
We solve the corresponding differential equations for the vector of master integrals M V = (M V 1 , . . . , M V 17 ) as an expansion in the parameter x using the methods presented in the previous section. The reduction to master integrals and the differential equations have been obtained in two different ways and found to agree. The differential equation are then normalized and expanded using our algorithm and the solution can be written as where the fundamental matrix of the system is given by is the transformation used to achieve the normalized form presented in section 4, R V is the residue matrix at x = 0 of the properly normalized differential equation matrix and B V = (B V 1 , . . . , B V 17 ) is the vector of boundary conditions. The Dyson series is defined as in (4.18). The differential equations, transformation matrix T V tot. and the corresponding boundary conditions can be found in the ancillary file attached to this article. The expansion (5.2) can then be unambiguously reconstructed using the method presented in section 4.
Since we are expanding around x = 0, the boundary conditions must be obtained from the asymptotic expansion of the integrals at x = 0, which are using the strategy of region [52,53]. Note that there are 7 factorizable integrals that can be written as a product of the following one-loop integrals with asymptotic behaviour such that the corresponding boundary conditions can be easily obtained. We were able to obtain the all the boundary conditions at all orders in and we see that they can be expressed using gamma functions and hypergeometric functions with unit argument, which can be easily expanded in using the package HypExp [35,36].
When expanded in , our solution reproduces the results of reference [12] in the small x limit and are given by

Master integrals appearing in the real contributions
In this section we describe the computation of the integrals contributing to σ R;int. gg→H , which are unknown in analytic form in the literature. They can be classified into the two topolo- gies shown in figure 3, which are defined as where dΦ 12→Hg = dΦ(p 1 , p 2 ; p 3 , p H ) denotes the 2 → 2 phase space. We parametrize it as with invariants given by Performing the reduction, we find 16 master integrals that we choose as follows: We solve the corresponding differential equations for the vector of master integrals M V = (M R 1 , . . . , M R 16 ) as an expansion in the parameter r = m 2 b /s using the methods presented in section 4. The reductions to master integrals and the differential equations have been obtained in two different ways and found to agree. The differential equation are then normalized and expanded using our algorithm and the solution can be written as M R (r, ) = F R (r, ) B R (z, ), (5.5) where the fundamental matrix of the system is given by where T R tot. = T R rank T R norm. is the transformation used to achieve the normalized form presented in section 4, R R is the residue matrix at r = 0 of the properly normalized differential equation matrix and B R = (B R 1 , . . . , B R 16 ) is the vector of boundary conditions. Note that the boundary conditions now depends on the variable z as well, making their evaluation much more difficult than in the case of the virtual master integrals. The differential equations, transformation matrix T R tot. and the corresponding boundary conditions can be found in the ancillary file attached to this article. The expansion (5.5) can then be unambiguously reconstructed using the method presented in section 4.
The main difficulty in evaluating the boundary conditions B R i (z, ) is the fact that the phase-space integral (5.4) induces a non-trivial analytic behaviour of the integrals as m b → 0. In the corresponding two-loop integrals, where no propagator is cut, one expects the integral to have scalings of the type (m 2 b ) −n , with n = 0, 1, 2. However, one-loop integrals can only develop scalings with n = 0, 1 such that the supplementary scaling with n = 2 must only appear after the phase-space integral has been performed. There is however no such prescription as the expansion by regions in the case where some of the propagators are cut. How can one then obtain the relevant boundary conditions for the (m 2 b ) −2 scalings that appear in the solution of the differential equations? We tackle this problem by introducing a Mellin-Barnes representation [59,60] and, owing to the simple structure of the one-loop integrals appearing in this computation, we are able to resolve the general asymptotic behavior of all the real master integrals.
Let us consider a generic one-loop integrals with n-external legs and whose propagators are all regulated by the same mass m b , The general parametric representation for such an integral is given by where F is the first Symanzik polynomial and ν = i ν i is the sum of the inverse powers of the propagators. Since all internal propagators have mass m 2 b , the first Symanzik polynomial must be of the form We can then introduce a Mellin-Barnes integral as where the integral is performed along the usual Barnes contour separating the poles at ζ = n from those at ζ = d/2 − ν − n, n ∈ N. Plugging this Mellin-Barnes representation in the parametric representation (5.6), one can reinterpret the integrand as being the very same one-loop integral but with m b = 0 and a shifted dimension d = 2(ζ + ν). Denoting this integral by I Note that one can not obtain a representation such as (5.8) at higher loops because in this case the second Symanzik polynomial U is not equal to one and the exponents of U and F do not allow one to recast the integrand after introducing the Mellin-Barnes representation (5.7) in a form where the original integral with m b = 0 appears. Equation (5.8) completely elucidates the structure of this class of integrals in the limit where m 2 b → 0. Applying Cauchy's theorem, the integral must be given by a sum over the residues of poles lying at the left of the contour, that is with (ζ) < 0. Each of these poles contributes to a specific power of m 2 b such that the complete expansion in m 2 b can be obtained by inspection of the poles of the integrand of (5.8) in the complex z-plane. There are two types of poles to be considered: • The poles coming from the overall gamma function Γ(ζ + ν − d/2) that are located at with n ∈ N. The corresponding residues will have an overall factor of (m 2 b ) −(ζ+ν−d/2) = (m 2 b ) n such that they all contribute to the hard region.
• The poles coming from the integral I (2(ζ+ν)) 0 itself. It is well known that Feynman integrals can develop poles only at even integer values of their dimension, such that these poles must be located at with n ∈ Z. The corresponding factor will then be (m 2 b ) −(ζ+ν−d/2) = (m 2 b ) n− , such that these poles contribute to the (m 2 b ) − scaling.
Note that only poles lying at the left of the contour need to be considered, that is we must require n > 2 − ν. This means that the small m 2 b expansion of the massive integral I can be unambiguously obtained from the -pole structure of the massless integral with shifted dimension I (2(ζ+ν)) 0 .
As an illustration, we consider the first poles with n = 0. The residue of the overall Γfunction at ζ = 2 − − ν is easily seen to be I (4−2 ) 0 , i.e. the correct hard region. Further assume that the massless integral has an -expansion of the form where β and α are the coefficients of the first and second order poles at = 0. It is then easy to compute the residue at ζ = 2 − ν and obtain the first order in the small m 2 b expansion as where ψ( ) is the digamma function. Validity of this equation can be easily checked with the one-loop triangle integrals shown in equation (5.3). Note that in order for the pole at z = 2 − ν to contribute we need ν > 2 such that the tadpole and bubble integrals (with standard powers of the propagators) do not satisfy equation (5.9).
Generalization to higher orders is straightforward, such that we obtain the small m b expansion where β i and α i are the first and second order poles of I (d ) 0 at (4−2i−d )/2 = 0, respectively. Note that I (d ) 0 can not develop poles of order higher than two, because its integrand has only two independent singular limits (soft and collinear).
The Mellin-Barnes representation (5.8) also allows one to obtain the scalings (m 2 b ) −2 created by the phase-space integral. Taking the phase-space into account, one gets integrals of the form Note that the dimension (d) of the phase-space integral given in (5.4) does not coincide with the dimension of the shifted one-loop integral (d = 2(ζ + ν)). The trick is to perform the phase-space integrals before the Mellin-Barnes integral, creating new poles in the ζ complex plane. These poles must be summed and give rise to the aforementioned scaling.
As an illustration, let us consider a simple one-loop triangle I 3 ( , m 2 b ) similar to the one shown in equation (5.3) but with s → s 13 . Introducing our Mellin-Barnes representation, we get for the corresponding phase-space integral where we have ν = 3 for the triangle and we set d = 4 − 2 = 2(ζ + ν), i.e = −1 − ζ.
We can directly compute 12) where N is the phase-space normalization factor given in (5.4)  Hence we see that the phase-space integral has created the new scaling. Plugging (5.12) in our Mellin-Barnes representation (5.11) we can collect the first order poles and obtain The two first terms corresponding to the scalings (m 2 b ) −0 and (m 2 b ) − can be obtained by first expanding I 3 ( , m 2 b ) for small m 2 b and then integrating over the phase-space term by term. The new scaling (m 2 b ) −2 corresponds to non-commuting contributions and the phasespace integral must be performed first in this case. In some sense, the poles coming from the virtual and real integrals do not overlap in the Melling-Barnes representation. Further note that result (5.13) is finite in . This is expected because all the IR divergences are regulated by the internal mass m 2 b . Such a trick has allowed us to recover all the scalings missing from a naive application of the strategy of regions. In the general case, the one-loop integral I (2(ζ+ν)) 0 under consideration might be a complicated function of the phase-space variables s 1g and s 2g and carrying the phase-space integration explicitly might be very difficult. However, since we only need to compute the coefficient of the -poles created by the phase-space integral, it is sufficient to look at the asymptotic limit of the massless integrals for small x 1 (andx 1 ), where we must have that with m, n ∈ N. This expression can then be integrated over the phase-space and it is easy to see that only the second term contributes to the −2 scaling. In conclusion, all contributions to the −2 scaling can be obtained by first taking the asymptotic expansions of the massless integral I (d ) 0 in the limits where x 1 andx 1 are small, integrating over the phase space and finally taking the corresponding residue in the ζ plane.
As previously said, equation (5.8) allows us to obtain all the boundary conditions from the corresponding massless one-loop integrals. The most difficult massless one-loop integral to consider is the massless one-loop box with one off-shell leg (p 2 4 = m 2 H ) given by with t = s 23 , u = s 13 . See for example [61] for a derivation. Note that the box appearing in the second topology shown in figure 3 can be expressed from (5.14) by simple crossing. Even if the hard region are generally simple, the integral over the phase space can lead to complicated hypergoemetric functions. We were able to obtain all the boundary conditions at all orders in z and using the hypergeometric functions 2 F 1 and 3 F 2 , which can be easily expanded by using the Mathematica package HypExp [35,36] , together with the Meijer G function with arguments: The Meijer G function is best defined via its Mellin-Barnes representation as where the integral is performed along the usual Barnes contour. This representation allows one to obtain the -expansion of such integrals easily. As an example, let us consider the following Mellin-Barnes integral where γ denotes the usual Barnes contour. Representation (5.15) is not suitable for expanding in , because the right-pole at s = 2 and the left-pole at s = 0 would pinch the contour γ as → 0. It is thus necessary to define a modified contour γ that goes to the right of the pole at s = 2 , such that (5.15) becomes The integrand of the Mellin-Barnes integral can now be safely expanded in . Integrating order by order produces harmonic sums that can be rewritten in terms of harmonic polylogarithms using standard techniques [62,63] 7 . We obtain G 3,2 3,3 A complete -expansions of the Meijer G-functions appearing in this computation can be found in the ancillary file attached to this article. Before presenting our results for the master integrals let us comment briefly on the issue of the soft limit z → 1. In principle, we want to obtain the small m 2 b expansion of the total hadronic cross section (2.2) such that the integral over z must still be taken into account.
Extending the ideas presented in this section, we can expect this last integral to create a new non-trivial analytic behaviour as m b → 0 and induce a new scaling (m 2 b ) −3 . Indeed we find that this is actually the case for single terms in the matrix elements. In particular, we looked at the master integral M R 3 (multiplied by the its coefficient from the gg → gH matrix element) and we see using the explicit Mellin-Barnes representation introduced earlier that the z integral indeed creates a non-vanishing (m 2 b ) −3 scaling. However, we do not need to worry about such subtleties when considering the full matrix element because the correct soft limit (with massive b-quarks) is recovered at all orders in by simply expanding under the z integral, see section 3 for more details. This means that even if individual masters multiplied by their coefficients do develop terms proportional to (m 2 b ) −3 when integrated over z, the physical matrix element does not have such contributions because a soft subtraction can be implemented as usual. Such a cancellation should also work at N 2 LO and in particular the correct soft limit z → 1 should also be obtained there by simply expanding under the z integral.

Matrix element appearing in the virtual contributions
In this section we present our results for the unrenormalized squared matrix element for the process gg → H in terms of master integrals needed for the computation of the virtual contributions. The amplitude consists of two pieces: the first piece is the leading order of this process in the full theory with bottom quarks running in the loop interfered with one loop correction to the ggH vertex in the effective theory, while the second piece is the two-loop correction to the to the same process in the full theory interfered with the Higgs-gluon-gluon vertex in the effective theory. See section 2 for more details.
We write the amplitude for the process gg → H generically as where P is the following helicity projector: where µ (p 1 ) and ν (p 2 ) are the polarisation vectors for incoming gluons with momenta p 1 and p 2 respectively, and a and b are the colour indices of the gluons. Then the interference of two amplitudes can be written as where we used The first piece, the leading order for the production of a Higgs boson via gluon fusion in the full theory interfered with the NLO matrix element in the effective theory is then given by  [64]. The second piece, the double virtual amplitude for the production of Higgs boson via gluon fusion in the full theory interfered with the LO matrix element in the effective theory is given by We present A V below in terms of master integrals. For convenience we have set s = 1.

Matrix elements appearing in the real contributions
The QCD corrections to the Higgs production via real emission involves partons in the final state, which means we have to take into account the processes involving the light quarks in the initial and final states.

Quark Channels
For the real emission involving light quarks we consider the following channels: q(p 1 ) +q(p 2 ) → H(p H ) + g(p 3 ) (6.2) q(p 1 ) + g(p 2 ) → H(p H ) + q(p 3 ) (6.3) q(p 2 ) + g(p 1 ) → H(p H ) + q(p 3 ) (6.4) q(p 1 ) + g(p 2 ) → H(p H ) +q(p 3 ) (6.5) q(p 2 ) + g(p 1 ) → H(p H ) +q(p 3 ) (6.6) Quarks are in the initial state in the first channel, whereas the other channels have quarks both in initial and final states. The latter channels have the same cross section since the cross section is invariant under the exchange of two initial state momenta, p 1 and p 2 . Therefore we present only one result below. Similar to the gluon channel we write the squared matrix element as create unregulated divergences. Compared to our original algorithm this amounts to a redefinition of the boundary conditions. The computation of the inverse of the fundamental matrix F −1 (r) as a series in r can be easily obtained by the method of series reversion. Up to the pre-factor of r R that can easily be inverted, the fundamental matrix has an expansion of the form where F n are matrix coefficients that do not depend on r. We are looking for a leftmultiplicative inverse G(r) of F(r) with G(r)F(r) = 1: Setting G(r) = n G n r n , Cauchy multiplication formula directly gives G 0 F 0 = 1, . . .
such that the matrix G(r) can be defined order by order as .3 can then be applied iteratively to all the subsystems of the differential equation under consideration to generate a solution of the complete differential equations.