Two-loop gg → Hg amplitude mediated by a nearly massless quark

We analytically compute the two-loop scattering amplitude gg → Hg assuming that the mass of the quark, that mediates the ggH interaction, is vanishingly small. Our computation provides an important ingredient required to improve the theoretical description of the top-bottom interference effect in Higgs boson production in gluon fusion, and to elucidate its impact on the Higgs boson transverse momentum distribution.


Introduction
The recent discovery of the Higgs boson at the LHC and the close proximity of its properties to the Higgs boson of the Standard Model, strongly suggest that electroweak symmetry breaking is triggered by an elementary scalar field.However, since the Higgs sector of the Standard Model is not very compelling and since there is a large number of theoretically appealing alternatives, experimental exploration of the Higgs boson properties has a very high priority for elementary particle physics.It is expected that the results of this research program will allow us to sharply contrast the observed properties of the Higgs boson with the Standard Model expectations.Because the Standard Model is a renormalizable theory, we can predict the expected properties of the Higgs boson with a precision that is only limited by our ability to perform the required computations and by the lack of understanding of non-perturbative phenomena that affect the outcomes of hadron collisions.The latter issue will probably prevent us from doing sub-percent precision physics for the Higgs couplings, but it is irrelevant for studying the proximity of the Higgs couplings to their Standard Model values with a few percent precision.Reaching this, the few percent, precision in theoretical predictions for Higgs physics is non-trivial; it requires many ingredients including improved perturbative predictions for major Higgs boson production processes.Providing such predictions is the main motivation behind the computation reported in this paper.
The major Higgs boson production mechanism at the LHC is the gluon fusion gg → H.The interaction of the Higgs boson with gluons is mediated by quarks; since the Higgs Yukawa couplings are proportional to quark masses, the top quark provides the dominant contribution to the ggH interaction vertex.Significant theoretical advances in describing this contribution enabled the prediction of the strength of the ggH interaction in the limit of a very heavy top quark m t → ∞ with a residual uncertainty of about 4% [1].At this level of precision, many other effects have to be taken into account; a detailed discussion can be found in Ref. [1].One such effect, that we focus on in this paper, is the modification of the gg → H interaction strength by loops of bottom quarks. 1t may sound surprising that we need to care about the bottom quark loop contribution.Indeed, simple power counting indicates that the bottom quark contribution is suppressed relative to the top quark contribution by m2 b /m 2 H ∼ 10 −3 .However, a more careful analysis of the bottom quark contributions reveals that it is enhanced by two powers of a logarithm ln(m 2 H /m 2 b ) ∼ 6.5.As the result, the relative suppression of the bottom quark loop relative to the top quark loop becomes much weaker, m 2 b /m 2 H ln 2 (m 2 H /m 2 b ) ∼ 10 −1 , making a detailed understanding of the bottom quark contribution quite relevant at the few-percent precision level.Since the NLO QCD corrections to gg → H are known to be significant, it is gratifying that these corrections are available for an arbitrary relation between the quark mass and the mass of the Higgs boson [2,3].
The situation becomes more complex if we consider less inclusive quantities, for example the transverse momentum distribution of the Higgs boson or the cross section for the production of the Higgs boson in association with a jet.In this case, the double logarithmic enhancement becomes p ⊥ -dependent, i.e. some of the large ln(m 2 H /m 2 b ) logarithms turn into ln(p 2 ⊥ /m 2 b ).These p ⊥ -dependent logarithms represent a serious problem for p ⊥ -resummations since their origin and their structure in high orders of QCD are not understood. 2In the absence of a clear understanding of how to resum these terms, the extent to which these p ⊥ -dependent logarithms affect the Higgs transverse momentum distribution was studied empirically [6][7][8][9].The results of these studies indicated a few percent differences in predicted transverse momentum distribution of the Higgs boson, depending on how these non-canonical log(p ⊥ /m b ) terms are treated in the resummed calculations.
It is interesting to put these studies on a more solid ground.We believe that a good starting point is the computation of the scattering amplitude for gg → Hg process in the approximation where the mass of the quark that mediates the ggH interaction is treated as the smallest kinematic parameter in the problem.Indeed, such a computation will give us a solid perturbative result that can be used directly to improve the prediction for the transverse momentum distribution of the Higgs boson in the region m b ≪ p ⊥ and, at the same time, an interesting data point for attempting the resummation of the Sudakov-like logarithms described above.The goal of this paper is to compute the two-loop gg → Hg amplitude in the approximation m b → 0.
We remark that the computation of the gg → Hg scattering amplitude for a nearly massless internal quark is an interesting theoretical challenge.Indeed, in contrast to the limit of a large internal quark mass, there is no algorithmic procedure to expand the Feynman diagrams that contribute to gg → Hg around the vanishing quark mass.It is possible to get around this problem by delaying the expansion in the small quark mass until it becomes clear how it can be performed.If one thinks about this problem keeping in mind the established technology for higher-order computations that includes i) generation of Feynman diagrams; ii) projection of scattering amplitude on Lorentz-invariant form factors; iii) reduction of contributing Feynman integrals to master integrals and iv) derivation and solution of the differential equations for master integrals, it is easy to realize that the safest point to start the expansion in the small quark mass is when the differential equations for the master integrals are about to be solved.Indeed, since the differential equations contain all the information about the singular behaviour of the master integrals in the limit of the vanishing quark mass, it should be possible to solve these equations by expanding the solutions around this singular point without any assumptions about the behavior of the integrals in the limit m b → 0. Proceeding along these lines, we can find the master integrals by consistently neglecting all the terms that are power-suppressed in the m b → 0 limit and, at the same time, keeping all the non-analytic O(log m b ) terms.This procedure was recently discussed in Ref. [10] in the context of the inclusive Higgs production in gluon fusion.In this paper we describe how to generalize it to the case of gg → Hg. 3However, we would like to stress that performing the expansion at the level of the differential equations is not optimal since the derivation of the differential equations requires the reduction to master integrals for an arbitrary relation between the quark mass and other kinematic parameters.As we explain in Section 3, these reductions are so demanding in terms of computing resources, that their successful completion should not be taken for granted.It is clear that, for our purposes, the full reduction is an overkill since we are interested in the limit m b → 0; nevertheless, it is non-trivial to take this limit consistently at the time of the reduction.It will be important to develop a computational method that will allow us to do that and we leave this interesting problem for the future.
The paper is organized as follows.In Section 2 we introduce our notation, de-scribe the parametrization of the scattering amplitude in terms of invariant form factors and explain how to apply the renormalization procedure to get the finite result.In Section 3 we describe how the invariant form factors are obtained from Feynman diagrams and how the contributing integrals are expressed in terms of master integrals.The master integrals are computed by solving the differential equations, as we discuss in Section 4. By solving the differential equations, we determine the integrals up to the integration constants.Some of these constants can be obtained by imposing regularity conditions on the solutions, but some can not and have to be computed separately.We discuss a few examples in Section 5. We present results for the helicity amplitudes in Section 6, discuss analytic continuation in Section 7 and conclude in Section 8.

The scattering amplitude
We consider the process mediated by a bottom quark loop in a theory that includes gluons, N f massless quarks and the bottom quark.The masses of the bottom quark and the Higgs boson are denoted by m b and m h , respectively.We define the Mandelstam variables subject to the constraint s + t + u = m 2 h , and note that for the process Eq.(2.1) all Mandelstam invariants are positive s > 0, t > 0 and u > 0. Following Ref. [12], we define the dimensionless variables The scattering amplitude is a function of x, y, z, κ with an overall multiplicative factor (−m 2 h ) −ǫ per loop.We would like to compute the scattering amplitude in the Euclidean kinematics.This is achieved by taking m 2 h < 0, s < 0, u < 0, t < 0 and keeping m 2 b > 0. In terms of x, y, z, this implies that for the scattering amplitude is explicitly real.We denote the scattering amplitude for the process (2.1) as where f a 1 a 2 a 3 are the SU(3) structure constants and ǫ j (a j ) is the polarization vector (color label) of a gluon with momentum p j , j = 1, 2, 3. Using Lorentz symmetry and gauge invariance, one can show that the scattering amplitude A is given by a linear combination of just four form factors.In particular, using the transversality conditions ǫ i • p i = 0, i = 1, 2, 3, and imposing a cyclic gauge fixing condition we can write the amplitude tensor in the following way (2.7) The four form factors F 1,..,4 (s, t, u, m b ) are Lorentz scalars; they admit a perturbative expansion in the QCD coupling constant.The expansion of the unrenormalized form factor reads where α 0 is the bare QCD coupling constant.
To perform the ultraviolet (UV) renormalization of the form factors we proceed as follows.First, we express the bare coupling constant and the bare bottom quark mass through their renormalized values.Note that this also applies to the renormalization of the Yukawa coupling since g Y = m b /v, where v is the vacuum expectation value of the Higgs field.Second, each of the form factors is multiplied by the gluon field renormalization constant raised to an appropriate power.This is sufficient to perform the ultraviolet renormalization.
In practice, we subtract the contribution of the massless quarks to the bare coupling constant in the MS scheme and the contribution of the b-quark to the bare coupling constant at zero momentum transfer.This implies the following relation between the bare coupling constant α 0 and the renormalized one at the scale µ R , where is the number of massless quark species employed in the computation and The quark mass renormalization is performed by replacing the bare quark mass with the on-shell renormalized quark mass.Technically, this amounts to making the following substitution in the form factors where and expanding them to the appropriate order in the strong coupling constant.We find Finally, the gluon wave function renormalization is performed by multiplying every form factor by with δ w defined after Eq.(2.9), and expanding in α s .With these notations, the UV-renormalized form factors become with (2.13) The UV-renormalized form factors still contain poles in ǫ that are of soft and collinear origin.For a generic NNLO QCD scattering amplitude, it was shown in Ref. [13] that all such poles can be written in terms of tree-and one-loop amplitudes of a given process.Although we perform a two-loop computation, the one-loop amplitude for gg → Hg is the leading term in the perturbative expansion; for this reason, it is sufficient to use the NLO QCD results of Ref. [13], to isolate the infra-red and collinear poles.We therefore write ,fin j = F (2),UV j where in the case of three external gluons the I 1 (ǫ) operator reads (2.15) We perform the UV renormalization at the scale µ 2 R = m 2 h .We verified explicitly that the IR poles in the form factors are removed by the subtraction in Eq. (2.14).

Calculation of the form factors
The direct computation of the decay amplitude A, using the standard methods for multi-loop computations, is difficult because the amplitude depends on the polarization vectors of the external gluons.We can get around this problem by computing the form factors instead.To this end, we design projection operators to extract contributions of different Feynman diagrams to the four form factors.We define four projection operators P µνρ j by requiring that they satisfy the following equation where, for consistency with Eq.(2.6), sums over polarizations of external gluons are taken to be We stress at this point that all Lorenz indices in Eq.(3.1) have to be understood as d-dimensional.The explicit form of the projection operators can be found by making an Ansatz in terms of the same linearly independent tensors as in Eq.(2.7) where j ∈ {1, 2, 3, 4}.The scalar functions c (j) i are unknown a priori; they are found by requiring that Eq.(3.1) is satisfied.We obtain With these results at hand, we can compute each of the form factors separately.Since the form factors are independent of the external polarization vectors, all the standard techniques employed for multi-loop computations can be applied.In practice, we proceed as follows.We generate the relevant one-and two-loop Feynman  (3.2, 3.3, 3.4).Once this step is completed, each contributing diagram is written in terms of integrals that depend on the scalar products of the loop momenta between themselves and the scalar products of the loop momenta with the external momenta.We can assign all Feynman integrals that contribute to the scattering amplitude to three integral families, two planar and one non-planar.These integral families are given by where top ∈ {PL1, PL2, NPL} is the topology label and the propagators [1], [2], ..., [9] for each topology are shown in Table 1.The integration measure is defined as We note that the loop momenta shifts required to map contributing Feynman diagrams on to the integral families are obtained using the shift finder implemented in Reduze2 [15].All algebraic manipulations needed at different stages of the computation are performed using FORM [16].Once the amplitude is written in terms of scalar integrals, we simplify them using all possible loop momenta shifts with a unit Jacobian; this can also be done using the momentum shift finder of Reduze2.When the contributions of all diagrams to the form factors are summed up, significant simplifications occur; for example, only integrals with up to three scalar products are left, although some individual diagrams receive contributions from integrals with up to four scalar products.
Having determined all scalar integrals that contribute to the amplitude, we need to reduce them to master integrals.The reduction procedure relies on a systematic Prop.Topology PL1 Topology PL2 Topology NPL Table 1: Feynman propagators of the three integral families.
application of the so-called integration by parts identities (IBPs) [17,18] to the integrals that belong to the three topologies defined in Table 1.This procedure is automated so that, as a matter of principle, one can use the publicly available programs Reduze2 [15,[19][20][21], FIRE5 [22,23] and LiteRed [24] to perform the reduction.However, in practice, the reduction to master integrals appears to be very challenging, due to the presence of a mass parameter in some of the propagators.We stress that, although we will eventually obtain the result for the amplitude assuming that the mass parameter is small, we retain the full mass dependence at the intermediate stages of the computation, including the reduction to master integrals.We have found that the publicly available reduction programs and, in some cases, also their private versions, were unable to successfully complete the reduction of the most complicated non-planar 7-propagator integrals.In order to reduce those integrals, we wrote a FORM program that produces and solves the IBPs for the 7-propagator non-planar integrals, thereby reducing them to 6-propagator integrals.We found that this step by itself is relatively straightforward and not too timeconsuming; however, once the reduction of the produced 6-propagator integrals is attempted, the reduction procedure stalls.In order to simplify this step as much as possible, we perform it only at the level of combinations of integrals that are actually required for computing the amplitude or the differential equations for the master integrals, but not for all contributing integrals individually.Moreover, we want to work with as compact expressions as possible and we do this by choosing wisely the basis of the 7-propagator master integrals.The main criterion that we impose is that, upon reduction, all 7-propagator integrals are written in terms of master integrals, whose coefficients do not contain any non-factorizable unphysical poles which mix the Mandelstam variables (s, t, u) and the space-time dimensions d.We searched for the right basis empirically, by fixing the Mandelstam variables s, t  and u to numerical values and performing the reductions with our code. 4On one hand, using numerical values for the Mandelstam variables makes the reduction very fast; on the other hand, it does not change the dependence of the final result on the space-time dimension.We found many different bases which fulfill the d-factorization requirement and we have chosen the basis that leads to the differential equations with the nicest properties, as explained in detail in Section 4. The steps described above allowed us to express all the integrals that contribute to the scattering amplitudes in terms of master integrals and to derive the differential equations for master integrals retaining full dependence on the internal quark mass.We will now explain how these differential equations are used to compute the master integrals.

Solving the differential equations for the master integrals
Following the procedure outlined in the previous Section, we write the form factors as linear combinations of the master integrals.Examples of master integrals with seven propagators that need to be computed are shown in Figs. 2, 3, 4. Master integrals with six or less propagators are obtained from the ones with seven by removing some of the internal lines.
To compute these master integrals, we consider their derivatives with respect to the kinematic variables that they depend upon.These derivatives are given by Feynman integrals that belong to the integral families that we discussed in the previous Section; for this reason they can be expressed through the master integrals.Following these steps, we obtain a system of partial differential equations that the master integrals satisfy.
Derivation of the differential equations is facilitated by the fact that derivatives with respect to kinematic invariants can be written as linear combinations of derivatives with respect to the four-momenta of external particles; the latter derivatives can be easily computed if we use Eq.(3.7) to write down the master integrals.Specifically, treating s, t and u as independent variables, we obtain where Partial derivatives with respect to y = t/m 2 h and z = u/m 2 h at fixed m 2 h are then related to the partial derivatives in Eq.(4.1) in a straightforward manner The partial derivative with respect to the the b-quark mass is trivially related to the The differential equations are computed by applying Eq.(4.2) and Eq.(4.1) to the master integrals and using the integration-by-parts identities to reduce the resulting integrals to master integrals.In this way, coupled systems of differential equations in κ, y and z are found for the list of master integrals that we denote by {I i } throughout this section.
We can also compute the derivatives of the master integrals with respect to m 2 h .However, these differential equations are not useful since, when the integrals are expressed in terms of the variables κ, y, z and m 2 h , the m 2 h -differential equations trivialize and only provide the (already known) information on the canonical mass dimensions of the master integrals.For this reason, we set m 2 h = 1 at the beginning and re-introduce it back at the very end of the computation.
The differential equations in κ, y and z take the following form Matrices A k are rational functions of κ, y, z and ǫ.It is essential that these matrices are sparse and, to a large extent, triangular.This allows us to organize the process of solving the differential equations by starting from the simplest integrals with smaller number of propagators and gradually moving to more complex ones.The two-loop tadpole integral is computed independently and used as an input for the differential equations for master integrals with three or more propagators.We are interested in solving the differential equations as an expansion in the normalized b-quark mass squared κ, around κ = 0. From the structure of the differential equations it follows that κ = 0 is a singular point; as the consequence, we have to look for the solutions of the differential equations using the following Ansatz I i (κ, y, z, ǫ) = j,k∈Z,n∈N c i,j,k,n (y, z, ǫ) κ j−kǫ log n (κ). (4.4) In practice we observed that for the computation of the master integrals that appear in the form factors, a simpler Ansatz is sufficient5 As indicated in Eq.(4.5), the strongest κ-singularity that we encountered in the master integrals is κ −1 ; this is related to the fact that the master integrals that we have chosen have at most one propagator raised to the second power.As a rule, integrals with higher powers of propagators have stronger κ-singularities.In principle, since we are interested in the computation of the scattering amplitude in the limit κ → 0 and since the scattering amplitude has at most logarithmic κ-singularities in this limit, it would have seemed natural to choose master integrals with similar or weaker singularities.However, if this is done, it becomes more difficult to solve the system of differential equations.This can be understood by a closer proximity of integrals that have similar singularity structure in κ, compared to integrals whose structure of singularities is very different.
The coefficients c i,j,k,n are functions of y, z, ǫ.To determine them, we start with the differential equations with respect to κ.We use the Ansatz Eq. (4.5) in κdifferential equations and require that the coefficients of the κ j−kǫ log n κ terms vanish independently of each other.This gives a system of linear algebraic equations for the coefficients c i,j,k,n that can be solved straightforwardly.
Suppose that a particular sector has N coupled master integrals.Upon solving the κ-differential equations in this sector, we are with N unknown integration "constants" that are, in fact, functions of y and z.If in the massless case 6 there are N 0 master integrals in this sector, then there are N 0 integration constants that can be determined by matching the massive results to the massless ones.The massless limit of the integrals that we study in this paper was computed in Refs.[25,26] and can be borrowed from there.
The remaining N −N 0 coefficient functions have to be determined by considering the differential equations in y and z.We use the Ansatz Eq. (4.5) in the (y, z)differential equations and again demand that the coefficients of the κ j−kǫ log n κ terms vanish; this gives the required N − N 0 (y, z)-differential equations for the coefficient functions.We solve these differential equations order by order in ǫ.We find that, similar to the massless case [25,26], the coefficient functions can be expressed in terms of Goncharov polylogarithms (GPL).The GPL's are defined through the iterative formula subject to additional constraints The denominators that appear in the recursive integrands of the GPL's in Eq.(4.6) are the ones that appear in the matrices A k after expanding around κ = 0; they assume the following values It is easy to see from the definitions of GPL's that these denominators lead to branch points at x = 1 − y − z, y, z = 0, y = 1, z = 1 and y = −z.Physically, only the first three singular points are allowed while the last three are not.This implies that the corresponding GPL's can appear in the results for the master integrals only in such combinations where these unphysical singularities cancel.As we explain below, this feature allows us to simplify calculation of master integrals in certain cases.We expand the coefficient functions c i,j,k,n in ǫ through the weight four. 7We also adjust the expansion of the master integrals in κ in such a way that the leading O(κ) contribution to the amplitude can be computed.Note that this requires expanding some of the integrals to relatively high order in κ since some master integrals appear in the differential equations with coefficients that scale as κ −n , n > 0, in the κ → 0 limit.
Upon solving the differential equations, we write the solutions for each of the master integrals in the following way (4.9) The lower limit of the ǫ expansion, r , is bounded below by −4.We include ancillary files with the paper that contain all master integrals required for computing the form factors, expanded in κ and ǫ as in Eq. (4.9).Note that this form of the solution gives access to different κ-branches, i.e. terms O(κ −2ǫ , κ −ǫ , κ 0 ).Each of these branches should correspond to contributions of distinct "regions", using the language of the "strategy-of-regions" [27,28], or "modes", in the language of effective field theories.Knowing the results for each of the κ-branches separately should be useful for understanding how to resum the log(κ) terms.Note also that individual branches have stronger ǫ-singularities than the complete integral.
Finally, we note that after solving the (y, z)-differential equations, we can only determine the master integrals up to the constants of integration, that need to be computed separately.We explain how to do this in the next Section.However, once this is done, we have the complete expression for the master integral and we can check it by comparing numerically the expansions in κ and ǫ of all the master integrals in various kinematic points in the Euclidean region with FIESTA [29].For all integrals required for the calculation of the gg → Hg amplitude, we found a perfect agreement between the analytical results obtained in this paper and the numerical results obtained with FIESTA.

Boundary conditions
By solving the differential equations, we can only obtain the master integrals up to the integration constants.These constants have to be determined separately.To this end, it suffices to compute the master integrals at an arbitrary kinematic point and then match the results to the solutions of the differential equations.However, this procedure is our last resort since, usually, there are other, simpler, ways to obtain the required integration constants.
We have already mentioned some of them in the previous Section.For example, if we are able to determine a master integral from the κ-differential equation, the integration constant is the massless branch of a particular integral known from Refs.[25,26].
Another way to determine the integration constants arises if the homogeneous parts of the (y, z)-differential equations exhibit unphysical singularities in x, y and z variables.Since, upon integration, singularities of differential equations become singularities of master integrals and since only certain, physical, singularities in x, y and z can appear in master integrals, we determine some of the integration constants by requiring that the unphysical singularities of the differential equations do not appear in the master integrals.
When none of the above applies, the integration constants have to be computed by matching the value of an integral to the solution of the differential equations for some values of y and z.It is difficult to describe how this has been done since there is no method that covers all the cases.In practice, we have used different techniques such as integration over Feynman parameters, expansion-by-regions, Mellin-Barnes integration and fitting numerical values of the integrals to linear combinations of the fixed-weight irrational numbers, to determine integration constants for master integrals.We give a few examples below to illustrate how these techniques are applied.
Consider the integral I PL2 (1, 1, 0, 0, 1, 0, 1, 0, 1).It is given by the following expression . (5.1) In the limit of small κ = −m 2 b /m 2 h , the integral behaves as 2) The singularity at x = 1 − y − z = 0 is allowed and the two constants of integration C 1,2 can not be determined from the differential equations.To find these constants, we need to extract the non-analytic terms that arise in the limit κ → 0. We do this by first re-writing the integral over the loop momenta through an integral over Feynman parameters.We obtain (5.4) Inspecting the above equations, we conclude that the two branches, κ −ǫ and κ −2ǫ , appear due to the integration over two distinct regions To project the integrand on one of the two branches, one should Taylor expand the integrand in a variable that is small for a particular branch and extend the integration over this variable to the positive half of the real axis.Upon applying this procedure, we arrive at the following expression for the branch κ −2ǫ and for the constant C 2 (5.7) The integration over ξ and α can be easily performed and we obtain (5.8) Upon writing we split the integral into two parts; the first integral can be evaluated exactly and the second can be evaluated upon expanding in ǫ.We find (5.10) To obtain the second constant, we need to understand how the branch O(κ −ǫ ) arises in the integral.As we already pointed out, this branch corresponds to the region where 1 − y ∼ κ.We therefore change variables y → 1 − f , Taylor expand the integrand in f assuming that f ∼ κ ≪ α, µ, ξ and extend the integration over f to the positive half of the real axis.We find Integrations over f, ξ and µ are straightforward and we obtain (5.12) As the second, more complicated example, we consider the evaluation of one of the non-planar, seven-propagator integrals.The integral reads (5.13) The constants of integration that need to be determined are in the branch κ −ǫ , including the coefficient of the logarithmic term κ −ǫ log κ.To determine these constants, we compute the integral at a particular kinematic point z → 1, y → 1.This implies that x = 1 − y − z = −1, so that the integral receives an imaginary part.
To compute the integral Eq.( 5.13) at the kinematic point described above in the limit κ → 0, we write it as an integral over suitably chosen Feynman parameters.In particular, we start by combining two pairs of propagators into single propagators and then shift the loop momenta l → l + k − ηp 1 and k → k + (1 − η)p 1 + p 3 .The integral becomes where p η 1 = (1 − η)p 1 and p H = −p 1 − p 2 − p 3 .Next, we integrate over l and k, substitute z = 1, y = 1 and arrive at the following representation for the integral where We notice that the integrand is linear in η and ξ.We also notice that the coefficient of O(η) term in ∆ is proportional to 1 − ξ; this means that upon integration over η, it cancels the (1 −ξ) factor in the integrand Eq.(5.16).As the result, the integrations over η and ξ can be performed exactly.Integrating over η and ξ, we arrive at the following represention of the non-planar integral at the kinematic point z = 1, y = 1 where To determine the coefficient of the κ −ǫ−1 branch, that arises in the limit κ → 0, we need to consider two integration regions Region 1 : As always, we perform the Taylor expansion of the integrand in all the variables that are parametrically small and then extend the integration region over the small variables to the positive half of the real axis.An important comment worth making is that, in order to be able to treat all the different terms in Eq.(5.18) separately, one has to introduce an additional regulator; in particular, the prefactors (1 − ρ(1 − t)) and (1 − ρ) in Eq.(5.18) must be modified in the following way (5.21) Constructing the expansion of Eq.(5.18) along the lines sketched above, we find that the κ −ǫ−1 branch of lim κ→0 I NPL (2, 1, 1, 1, 0, 1, 1, 1, 0) is obtained as the sum of five integrals.We write All of these integrals involve integration over some "small" variables that can be easily performed since the corresponding integration region extends from zero to infinity.We present the expressions for the five integrals after these simple integrations are carried out.We find (5.23) We compute these integrals directly.We note that through weight three, it is straightforward to calculate them analytically.However, at weight four, we compute some of the integrals numerically and then fit them to a linear combination of the irrational constants of weight four that include log 2 to an appropriate power and Li 4 (1/2).Working through O(ν 0 , ǫ 1 ), we obtain8 (5.24) Using these results in Eq. (5.22), we obtain the coefficient of the κ −1−ǫ branch of the integral lim κ→0 I NPL (2, 1, 1, 1, 0, 1, 1, 1, 0) that arises in the κ → 0 limit, at the point z = 1, y = 1.This result is then matched to the solution found for this branch from the (y, z)-differential equation and the constant of integration is determined.Other branches of this integral, as well as other integrals, that require determination of the integration constants can be studied along similar lines.

Helicity amplitudes
In this Section, we present the results of the computation of the H → ggg scattering amplitude.It is convenient to write the result fixing gluon helicities.We obtain the helicity amplitudes from the general expression for the H → ggg amplitude given in Eqs.(2.5, 2.7).We write where λ 1,2,3 are helicity labels of gluons with momenta p 1,2,3 respectively, and the dependence of the amplitude on the b-quark mass is suppressed.Since each gluon has two independent polarizations, there are eight helicity amplitudes in total but only two of them are independent.The remaining six amplitudes can be obtained by permuting the external gluons and applying parity and charge conjugation.We choose the two amplitudes A ++± (s, t, u) as independent and write them using the spinor-helicity notations. 9The polarization vectors for external outgoing gluons with momentum p j and the reference vector q j read We note that the reference vectors q j must be chosen for each gluon in accord with Eq. (2.6).Using these notations, we write the two independent helicity amplitudes as The helicity coefficients Ω j are linear combinations of the form factors defined in Section 2; they read Similar to the form factors, the helicity coefficients Ω j can be written as an expansion in the strong coupling constant.We write where m b is the b-quark pole mass and α s is the strong coupling constant at the scale µ = m h defined in the theory with N f active flavors.The two-loop helicity coefficients are not finite, but their infra-red divergent parts are described by the Catani's formula Ω (2l) where the I 1 (ǫ) operator is defined in Eq.(2.15).We note that, originally, we defined the renormalized coupling constant in a scheme with N f active flavors since the contribution of the b-quark loop was subtracted at zero momentum.This is not optimal since we are interested in a kinematic situation where all momenta transfers are large compared to m b .In this case, the appropriate strong coupling to use is the one defined in the scheme with N f +1 active flavors.The relation between the two couplings at the scale µ = m h reads α We use this relation to re-write the helicity amplitudes using the strong coupling constant defined in a theory with N f + 1 active flavors.Since the relation between the two coupling constants is finite, it induces changes in the finite parts of the twoloop helicity amplitudes.We denote the helicity coefficients written with the strong coupling in the theory with N f + 1 active flavors as Ω.We obtain We also note that it is far from obvious that the on-shell renormalization scheme for the b-quark mass is, indeed, physically appropriate.The helicity amplitudes are proportional to m 2 b , where one power comes from the Yukawa coupling constant and the other from the helicity flip in the quark loop.It is most likely that the appropriate choice of the mass parameter related to the Yukawa coupling is the MS mass defined at the scale m h .However, the proper choice of the mass parameter responsible for the helicity flip is much less clear.It will be very interesting to understand the scale and scheme choices in the virtual amplitude that minimize the magnitude of the NLO QCD corrections to physical observables, for example to the Higgs transverse momentum distribution.We plan to return to this question in the future work.
Full results for helicity coefficients Ω ++± can be found in the ancillary files provided with this submission.Although, on the scale of known two-loop helicity amplitudes, our Ω's are quite compact, they are nevertheless complex.It is therefore instructive to study them in a few interesting limits, where they simplify dramatically.
One such limit is the soft limit.It describes a situation where the energy of one of the gluons in the H → ggg amplitude becomes small.We take the gluon with momentum p 3 to be soft; this implies the following hierarchical relations between the kinematic variables m 2 b ≪ t ∼ u ≪ m 2 h ∼ s.For the dimensionless variables introduced earlier, the soft limit implies κ ≪ y ∼ z ≪ 1.
It follows from Eq. ( 6.3) that in the soft limit the helicity amplitudes diverge as 1/ √ yz ∼ 1/ p 2 ⊥ , where we introduced the transverse momentum of the Higgs boson p 2 ⊥ = ut/s.This is the standard soft divergence present in any scattering amplitude.
The helicity coefficients Ω ++± develop logarithmic singularities in the soft limit.It is convenient to write these helicity coefficients in terms of the logarithms of the ratio of the bottom quark mass and the Higgs boson mass log κ = log (−m 2 b /m 2 h ), logarithms of the ratio of the transverse momentum p ⊥ divided by the bottom quark mass log (y z/κ) = log (−p 2 ⊥ /m 2 b ) and the logarithms of the ratio of two small parameters u and t, log(y/z) = log(t/u).To simplify the notation, we define In the soft limit, L ≫ 1, while τ and ξ, defined above, are quantities of order one.
Expanding the helicity coefficients in the soft limit and substituting N c = 3, we find Ω (1l),fin coefficients of leading and subleading logarithms.We find Ω (2l),fin +++ = L 4 5η 4 72 + 5 36 These results show that the structure of large logarithmic corrections in the collinear limit is complicated; it does not appear that any one-loop structures get iterated even at the level of the leading logarithms.

The analytic continuation
Up to this point, we studied the scattering amplitude in the decay kinematics.In this case the imaginary part is generated by an overall prefactor (−m 2 h − i δ) −n ǫ , where n is the number of loops.However, this is not sufficient; indeed, to study Higgs boson production in gluon fusion, we need to perform an analytic continuation of the scattering amplitude.The analytic continuation of the gg → Hg amplitude for the point-like ggH interaction vertex was described in Ref. [31] and we closely follow that paper in our discussion below.Following Ref. [31] The region (1a) + is the decay kinematic region that we considered so far in this paper; the region (2a) + is the kinematic region required to describe Higgs production in gluon fusion and the region (4a) + is needed to determine all helicity amplitudes for the Higgs production process from the two independent ones computed in the previous Section.In Ref. [31] it was shown how to start from the helicity amplitudes defined in region (1a) + and continue them to any other kinematic configuration, including regions (2a) + and (4a) + .While the analytic continuation of the spinor structures is straightforward, some care has to be taken in continuing the helicity coefficients Ω ++± , which are written in terms of harmonic and two-dimensional harmonic polylogarithms defined in [25].The analytic continuation of the polylogarithms can be achieved with a proper change of variables.
In the region (2a) + the Mandelstam variables are m 2 h > 0 , s > 0 , t, u < 0 , (7.4) and the analytic continuation from the decay kinematics is performed by providing small and positive imaginary parts to all positive Mandelstam variables.In particular, we require s → s + i δ.We define two auxiliary variables which fulfill the following condition In the region (4a) + , instead, we have m 2 h > 0 , s < 0 , t < 0 u > 0, with u → u + i δ.Similarly to Eq.(7.5), we define the following auxiliary variables One can show that, by changing the arguments of the polylogarithms from (y, z) to (u 1 , v 1 ) and (u 2 , v 2 ) in regions (2a) + and (4a) + respectively, and rewriting the result appropriately, one can extract all imaginary parts explicitly.In this way, the result for the helicity amplitudes can be rewritten in terms of explicitly real one-and two-dimensional harmonic polylogarithms of the new arguments u 1 , v 1 and u 2 , v 2 , respectively.The one-and two-loop helicity coefficients in the decay kinematics Eq.(7.1) and analytically continued to the scattering kinematics Eqs.(7.2,7.3) are available as ancillary files with the arXiv submission of this paper.

Conclusions
We computed the two-loop helicity amplitudes for the process H → ggg, mediated by a quark loop, in the approximation when the mass of the quark is small compared to other kinematic parameters of the process.The expansion in the small quark mass is used to solve the differential equations for the master integrals, while in all other steps of the computation no approximations are made.
The methodological results of this paper establish a framework that allows one to expand the scattering amplitudes around the limit where all or some of the particles in the quantum loops can be treated as nearly massless.There are many potential applications of this approach, including production of the Higgs boson at high transverse momentum and the electroweak corrections in the Sudakov regime.
On the phenomenological side, there are several ways in which the results of the computation reported in this paper can be used.The two-loop helicity amplitudes can be combined with the real emission process gg → Hgg to study modifications of the Higgs boson transverse momentum distribution due to the interference of top quark and bottom quark loops in the process pp → H + j in higher orders in QCD.Such a study seems especially worthwhile given the recent proposals to constrain charm and bottom Yukawa couplings from the kinematic distributions of Higgs bosons produced in proton collisions [32,33].
Furthermore, it has been recognized long ago that conventional transverse momentum resummation framework can not be applied to the production of the Higgs boson through the loop of light quarks since large non-universal double logarithmic corrections are generated by the virtual amplitude itself.Although these double logarithms are of the Sudakov origin, they are subtle as the process requires the helicity flip on the soft quark lines.It would be interesting and instructive to understand if the resummation of these Sudakov-like logarithms can be achieved.We expect that the computation of the helicity amplitudes reported in this paper and the limits of these helicity amplitudes described in the previous Section, will contribute towards this goal.

Figure 1 :
Figure 1: Examples of two-loop Feynman diagrams that contribute to the process gg → Hg. diagrams with QGRAF [14].A few examples of the two-loop Feynman diagrams that contribute to the gg → Hg amplitude are shown in Fig. 1.The projection operators are applied diagram by diagram and the polarization sums are computed following Eqs.(3.2, 3.3, 3.4).Once this step is completed, each contributing diagram is written in terms of integrals that depend on the scalar products of the loop momenta between themselves and the scalar products of the loop momenta with the external momenta.We can assign all Feynman integrals that contribute to the scattering amplitude to three integral families, two planar and one non-planar.These integral families are given by

Figure 3 :
Figure 3: Same as in Figure 2, but for the integrals of the PL2 family.

2 Figure 4 :
Figure4: Same as in Figure2but for the integrals of the NPL family.Note that in this case all the three integrals are irreducible.The leftmost integral does not contribute to the form factors because of the color structure of the corresponding Feynman diagrams.For this reason we do not compute it.
Figure 2: Seven-propagator Feynman integrals of the PL1 family, that appear in the form factors and require IBP reduction.The two integrals at the top are irreducible and correspond to two sectors in the family PL1 that contain master integrals with seven propagators.The integrals at the bottom correspond to reducible integrals.All momenta are incoming.