Two-loop amplitudes for processes $g g \to H g, q g \to H q$ and $q \bar{q} \to H g$ at large Higgs transverse momentum

We compute the two-loop QCD corrections to amplitudes for processes $g g \to H g, q g \to H q$ and $q \bar{q} \to H g$ in the limit when the Higgs transverse momentum is larger than the top quark mass, $p_\perp \gg m_t$. These amplitudes are important ingredients for understanding higher-order QCD effects on Higgs transverse momentum distribution at large $p_\perp$.


Introduction
The discovery of the Standard Model (SM)-like Higgs boson at the LHC five years ago got rapidly transformed into an active experimental exploration of this new particle. Indeed, a detailed knowledge of Higgs boson properties and its coupling to other particles is essential for understanding its role in the electroweak symmetry breaking and for early clues about physics beyond the Standard Model. Since in the SM the Higgs couplings to gauge bosons and matter particles can be computed theoretically to a very high precision, the existence of equally precise measurement program is crucial to search for differences between measurements and predictions that may then be interpreted as signals of physics beyond the Standard Model (BSM).
Unfortunately, most recent results from the Run II of the LHC show that the Higgs boson fits very well the expected profile of the SM Higgs particle and no signs of New Physics have been seen so far. These conclusions are so far limited by statistical and systematic errors that, on average, are in the O(15 − 20) percent range but can be much larger for certain couplings and cross sections. It is expected that during the Run II and the high-luminosity phase of the LHC, the precision of Higgs couplings measurements will significantly increase, reaching eventually a few percent accuracy.
This accuracy has to be matched on the theory side and we have seen quite very impressive accomplishments in refining predictions for major Higgs production and decay processes in recent years. For example, the inclusive Higgs boson production in gluon fusion is now known to an impressive next-to-next-to-next-to-leading order (N 3 LO) QCD in the infinite top quark mass limit [1] and the H+jets cross section has been computed through next-to-next-to-leading order (NNLO) QCD in the same approximation [2][3][4][5].
The approximation of an infinitely heavy top quark is justified as long as typical values of kinematic parameters relevant for particular cross sections are smaller than O(2m t ). Although this criterion is satisfied for the majority of events selected for both inclusive and H + j cross sections, there are good reasons to look at regions of phase-space where this condition is explicitly violated. For example, with the dramatic increase of statistics promised by the highluminosity run at the LHC, we will have access to Higgs transverse momentum distribution at high values of p ⊥ ≥ m t . This is a very interesting regime since, as a matter of principle, it allows us to disentangle two terms in the effective SM Lagrangian -the point-like Higgs coupling to gluons and the modification of the Higgs-top Yukawa coupling [6][7][8][9][10]. 1 Amazingly, first experimental attempts to explore Higgs boson production at high-p ⊥ have recently been undertaken [12].
To fully benefit from this opportunity, it is important to have as precise predictions for Higgs p ⊥ -distribution at large transverse momenta as possible. Since for computations at high p ⊥ ≥ m t , the Higgs coupling to gluons cannot be treated as point-like, all existing higher-order computations, including most recent NNLO QCD predictions for H + j production [2][3][4][5] are of little use. In fact, when mass effects are accounted for, the p ⊥ -distribution appears to be known only at leading order which, in this case, is determined by one-loop diagrams. Since NLO QCD corrections for processes with gluons in initial state are known to be large [13][14][15], it is quite conceivable that large corrections to Higgs transverse momentum distribution at high p ⊥ are to be found as well. Computing two-loop contributions to relevant amplitudes and setting up the stage for a full NLO QCD computation of the Higgs boson transverse momentum distribution at high p ⊥ is the main goal of this paper.
We note that the relevant two-loop amplitudes for a NLO computation of Higgs plus jet production mediated via a massive quark-loop were considered recently in Refs. [16,17]. However, in those papers the limit of a small quark mass m q ≪ m H ∼ p ⊥ was considered. This limit is relevant for the bottom quark contribution to effective ggH interaction vertex but it is not the right limit to describe high-p ⊥ regime of the Higgs boson production.
To address the high-p ⊥ case we impose the following hierarchy between kinematic variables and particle masses m 2 h ≪ m 2 t ≪ s, t, u. This result is then applicable to the case where the Higgs boson is produced via a top quark loop at high p ⊥ . 2 To compute the scattering amplitude in that limit, we will follow an approach developed in Refs. [16,17,19] and expand the relevant Feynman integrals in small parameters, namely in m 2 h /m 2 t and m 2 t /s, using the differential equations that these Feynman integrals satisfy. We note that the computation of relevant integrals for arbitrary Higgs and quark masses is still ongoing; planar master integrals have recently been computed in [20].
The remainder of the paper is organized as follows. In Section 2 we explain the notation, introduce the relevant amplitudes, explain their decomposition into invariant form factors and describe the renormalization. In Section 3 we discuss how form factors are computed. We explain how to calculate the master integrals with the differential equation method in Section 3.1. In Section 3.2 we provide an example of how integration constants for differential equations can be computed. The final results for helicity amplitudes are presented in Section 4. The amplitudes are originally computed in the kinematic region where t > 0, s, u < 0; in Section 4.1 we describe the analytic continuation to other relevant scattering regions. We conclude in Section 5. We include ancillary files with this submission that contain analytic results for all relevant amplitudes in the different kinematic regions.

The scattering amplitudes
Production of the Higgs boson in association with a jet at a hadron collider can occur in several different ways; the relevant partonic processes can be found by crossing the Higgs decay processes to the production kinematics. We consider all quarks in Eq.(2.1) as massless. The Higgs boson interaction with gluons and massless quarks is facilitated by loops of top quarks; this is the only quark that we consider massive in this article. Some examples of Feynman diagrams that contribute to (crossed versions of) processes shown in Eq.(2.1) are presented in Fig. 1. The goal of this paper is to compute two-loop contributions to scattering amplitudes for processes g g H g g g H g g q H q Figure 1: The one-loop Feynman diagrams that contribute to the quark-loop induced processes gg → Hg and qg → Hq.
in Eq.(2.1) assuming that the Higgs boson mass and the top quark mass are smaller than all other kinematic invariants. We start by defining the Mandelstam variables We trade four dimensionful Mandelstam variables for a dimensionful variable s and three dimensionless variables In the large transverse momentum region and in the limit of a small Higgs mass the following hierarchy of scales applies For the top quark and Higgs boson with masses m t ∼ 173 GeV and m h ∼ 125 GeV respectively, |η| ∼ 0.13 and can be treated as a small parameter.
A Euclidean region where all Mandelstam variables s, t or u are negative does not exist since |m 2 h | = |s + t + u| ≪ |s|, |t|, |u| in the kinematic region that is of interest to us. At least one of the Mandelstam variables has to be positive and without loss of generality we choose t to be positive and s, u negative. Furthermore we will compute our amplitudes initially in the region where m 2 h < 0 and m 2 t > 0, in other words the parameters will satisfy If we analytically continue to the region where m 2 h > 0, our results will represent the physical scattering processes All other production channels can be found from crossing and analytic continuation of the computed amplitudes in the region specified in Eq. (2.5), as we will describe in Section 4. Note that because the Euclidean region does not exist, all the amplitudes have explicit imaginary parts.
We follow Refs. [16,17] and define the partonic amplitudes corresponding to the processes shown in Eq. (2.1) as The color structure of the amplitudes is completely factorized and captured by the SU(3) structure constants f a 1 a 2 a 3 and the usual Gell-Mann matrices T a jk for the gluon and quark channels respectively. The color indices are denoted by a 1,2,3 and j, k for gluons and quarks, respectively. The gluon polarization vectors are transversal ǫ i · p i = 0, i = 1, 2, 3 and the spinors satisfy the massless Dirac equations / p 1 u(p 1 ) = / p 2 v(p 2 ) = 0.
To understand the Lorentz structure of the amplitude, we write it as a sum of parity conserving Lorentz tensors of relevant ranks. The amplitudes must furthermore satisfy the Ward identity which implies that an on-shell amplitude must vanish after replacing any of the gluon polarization vectors with their momenta. After imposing these constraints, the H → ggg and H → qqg amplitudes can be written as a sum of four (two) tensors, respectively. They read (2.9) The above decomposition corresponds to the cyclic gauge fixing condition for the gluon polarization vectors ǫ 1 · p 2 = ǫ 2 · p 3 = ǫ 3 · p 1 = 0. (2.10) The form factors F q,g j are scalar functions of the Mandelstam variables and the quark mass. In the following we will drop the upper index q and g for simplicity, unless they need to be explicitly specified.
The unrenormalized form factors F j can be expanded in the bare QCD coupling constant constant α 0 as The LO contribution F (1) j with the full dependence on the quark mass was calculated in Refs. [21,22]. In this paper, we will compute the two-loop contributions to form factors F (2) j assuming that the Higgs boson transverse momentum is large and the Higgs boson mass is parametrically smaller than the mass of the top quark. Some examples of two-loop diagrams that contribute to Higgs boson production in association with a jet are shown in Figs. 2 and 3.
The unrenormalized form factors that we compute have poles in ǫ = (d − 4)/2; these poles are of ultraviolet (UV) and/or infrared (IR) origin. We perform the subtraction of these poles in two steps. First we UV renormalize the above bare form factors F un j in Eq. (2.11) (2.12) We express the bare strong coupling constant and the top quark mass parameter in F un j in terms of renormalized parameters and we include for each external gluon the wave-function renormalization factor. The strong coupling constant gets renormalized in the mixed scheme; this implies that contributions of N f massless quarks are renormalized in the MS-scheme whereas top quark contributions are subtracted at zero momentum. The top quark mass is renormalized in an on-shell scheme. The corresponding formulas read (2.14) the number of colors. The wave-function and mass renormalization constants are Renormalization of the gluon wave-function is taken into account by multiplying the form factors with for each of the external gluons.
Following the described procedure, we express the UV-renormalized form factors in terms of bare ones. We find where i = q, g denotes the H → ggg and H → qqg form factors respectively. Unfortunately, even after the UV renormalization is performed, the form factors still exhibit poles in ǫ. These are the infra-red and collinear poles that appear in the virtual amplitude; they disappear once elastic and inelastic partonic processes are combined to compute physical cross sections. Since the structure of IR-singularities is universal [23] and since they, as we said, will eventually get cancelled against real emission corrections, it is useful to separate them in the two-loop amplitude. We write where again i = q, g and I q,g 1 (ǫ) are the so-called Catani operators Our final results for the form factors, {F fin j } are finite in the limit ǫ → 0. Note that in order to perform the final IR subtraction we require the one-loop amplitudes to order ǫ 2 .

Computing the form factors
The bare form factors are expressed in terms of Feynman diagrams that we produce with QGRAF [24] and independently with FeynArts [25]. We allow for massless external quarks and both massive and massless internal quark loops. Some examples of Feynman diagrams that one has to consider are shown in Figs. 2 and 3. We follow procedures outlined in Refs. [16,17] to express the form factors in terms of scalar integrals by applying projection operators as follows (3.1) Explicit expressions for projection operators can be found in Refs. [16,17]. Both FORM [26] and FormCalc [27] have been independently used to implement the algebraic manipulations related to the projection in d dimensions. The resulting form factors are expressed as linear combinations of scalar integrals Prop. Topology PL1 where the integration measure is chosen to be Scalar integrals that appear in the form factors belong to one of the three integral families that we refer to as {PL1, PL2, NPL}. Sets of propagators that define each topology are shown in Table 1.
After an amplitude is projected on a form factor, all scalar integrals are reduced to a set of master integrals (MI) using the integration by parts identities (IBP) [28,29]. The reduction has been previously performed in Refs. [16,17] using public versions of FIRE5 [30,31], Reduze2 [32][33][34][35] and an in-house routine written in FORM [26]. 3 The MIs are computed by solving differential equations in kinematic variables; the differential equations are solved perturbatively, expanding in the small parameters κ and η, as will be explained in Section 3.1.
We note that MIs contain logarithmic singularities ∝ log(m 2 t ) ∼ log (κ) as κ → 0. These are mass singularities that are expected to be present in the high-p ⊥ kinematics. In addition, there are Feynman integrals that develop logarithmic singularities ∝ log (η) ∼ log(m 2 h ) as η → 0; This happens whenever all the massless external partons couple directly to massless internal propagators, such as for example is the case for the top center diagram in Figure 2. The resulting MI which appear after pinching this diagram also contains logarithmic singularities ∝ log (η). For these MI it is possible to cut massless propagators in their corresponding diagrams such that the squared momentum flowing into the cut equals m 2 h and therefore we expect a singular behavior as m 2 h → 0. Note that the top loop itself always gets screened by the top mass and therefore the log (η) singularities are attributed to a specific scaling of the loop momenta running through massless propagators.
Since the Higgs boson always couple to top quarks, we expect that all the log(m 2 h ) singularities are the artifacts of computational procedure and that they should cancel in the final result for form factor. We have confirmed this expectation by an explicit computation.
Another interesting point is that three sectors of non-planar MI (one sector with six propagators and two top sectors with seven propagators) have integrals whose expansion around κ → 0 starts with non-integer powers of κ, i.e. I ∼ κ −1/2 ∼ m −1 t . This non-analytic behavior indicates contributions to scalar integrals beyond the standard soft-collinear paradigm. It is interesting to see, however, that none of these non-analytic terms survives in final results for physical amplitudes.
To conclude, after the results for MIs are used to calculate the unrenormalized form factors, the form factors are written as an expansion in κ and η lim m h m t →0,mt→0 The Yukawa coupling and the helicity flip in one of the quark lines contribute each a factor of m t , which results in the overall factor of κ in the above result. In Eq. (3.4) we retain terms that are leading powers in the squared Higgs mass η and up to next-to-leading power in the squared top quark mass. Since there are no logarithms in η in the final result, we could have put η → 0 from the beginning. However, in our computation we did not do this and kept m 2 h ∼ η = 0 throughout the calculation, 4 only cutting off the expansion of form factors to leading power after inserting the MIs. It was argued in e.g. Ref. [38] for the quark channels that an expansion to leading power in η gives a good approximation to the full amplitude with non-zero Higgs mass. We have checked this statement explicitly by comparing expanded and un-expanded one-loop amplitudes for both quark and gluon channels. We conclude that Eq. (3.4) is expected to provide a reasonable description of the form factors with non-zero Higgs and top mass. We will next describe how to compute the MIs using the method of differential equations.

Solving for the two-loop master integrals
The master integrals with seven propagators correspond to Feynman diagrams shown in Figs. 2; all other MIs that contain six or even less propagators can be obtained from the highest-level ones by pinching. We note that all the master integrals for H+jet production were recently computed in an approximation m 2 q=b ≪ s ∼ t ∼ u ∼ m 2 h , in Refs. [16,17]. In this paper we are instead interested in computing master integrals for high energies and transverse momenta m 2 q=t ≪ s, t, u in a situation when the quark mass is larger than the Higgs mass, m 2 h ≪ m 2 q=t . To derive differential equations, we start by taking derivatives of the integrals with respect to the kinematic invariants m 2 t , s, t, u. The derivatives with respect to the Mandelstam variables can be expressed in terms of linear combinations of derivatives with respect to the four-momenta of the external particles The planar master integrals corresponding to the case where m 2 h = η = 0 have been computed in [36,37].
Here we use the notation p i · ∂ p j = p i,µ ∂/∂p µ j . The derivatives with respect to dimensionless variables defined in Eq. (2.3) are related to above differential operators through the following equations We apply the derivatives in Eqs. (3.2,3.1) to the set of master integrals and use integrationby-parts identities to reduce all the integrals back to master integrals. This procedure leads to a linear system of coupled partial DE for all the MIs that we will denote in this Section by {I i }. After expressing the MIs in terms of the chosen variables, the derivative with respect to the Mandelstam variable s becomes trivial and provides the mass dimension of the MIs. Therefore, it suffices to solve the MIs for the case s = 1 and re-introduce it back at the end of the calculation. The DEs take the following form The matrices A κ,η,z are sparse and can be put in a triangular form. We may then solve the system starting from the simplest integrals, which then serve as inhomogeneous contributions to the DEs of integrals with more propagators. The integrals which depend on a single scale, e.g. the two-loop tadpole integrals, are computed independently and serve as an input for the DEs.
The three matrices A k are rational functions of η, κ, z and ǫ. The MIs have been chosen such that the dependence on space-time dimensionality d does not mix with the kinematic variables inside the denominators that appear in A k . The matrices have singularities at η = 0, −1/2, −1, in other words at m 2 h = 0, m 2 h = 2m 2 t and m 2 h = 4m 2 t respectively. The pole at m 2 h = 2m 2 t is expected to be spurious and can be avoided by taking canonical combinations of the MI as in [20]. At the point η = m 2 h = 0 there are singularities at κ = 0, −1/4, −(1+z)/4, −z/4, which corresponds to poles at m 2 t = 0 and s, t, u = 4m 2 t respectively. The latter three poles arise from the top threshold when the invariant mass of a pair of final state particles in the processes in Eq. (2.1) is equal to 2m t . These considerations imply that the matrices can be conveniently expanded in m 2 h /(4m 2 t ) = −η and 4m 2 t /s = −4κ and the DE then solved perturbatively in small η and κ. The order of expanding in η and κ is irrelevant. Furthermore, since the DEs have singularities at both η = 0 and κ = 0, the solutions are expected to contain terms beyond a usual analytic Taylor expansion in η and κ. The structure of the differential equations implies the following ansatz I i (κ, η, z, ǫ) = j,k,l,m∈Z,n∈N c i,j,k,l,m,n (z, ǫ) η j−kǫ κ l/2−mǫ log n (κ). (

3.4)
A more detailed analysis of the differential equations shows that at two loops there are at most two powers of κ −ǫ and log(κ) and at most one power of η −ǫ . The following simpler ansatz therefore suffices c i,j,k,l,m,n (z, ǫ) η j−kǫ κ l/2−mǫ log n (κ). (3.5) The maximal value for the powers j, l of the variables η and κ, respectively, are chosen such that we can expand the form factors to leading power in η and to next-to-leading power in κ.
We note that this requires computing some of the MI to higher suppressed powers in η and κ.
As we already alluded to in the paragraph above Eq. (3.4), we need to include powers of η −ǫ in the ansatz for exactly six MI, which all appear in the planar topology PL1. A detailed study of these six MI shows that they indeed have terms that scale as η −ǫ when η → 0 but there are no 1/η singularities. For all other MI, the expansion in η → 0 correspond to a simple Taylor expansion. We conclude that none of the MI have singularities in 1/η, which fixes the lower bound of the index j in the sum of Eq. The coefficient functions c i,j,k,l,m,n depend on z and ǫ. We determine them by substituting the ansatz for integrals in Eq.(3.5) into the differential equations and equating terms with the same powers of η, κ and log(κ) on both sides of the relevant equations. This procedure relates the c i,j,k,l,m,n coefficients to each other via a system of linear algebraic equations. We note that the DEs allow powers of η −1 in our complete ansatz in Eq. (3.4). Therefore, requiring that solutions to DEs do not contain poles in η provides additional relations between coefficient functions.
Some of the coefficient functions c i,j,k,l,m,n remain undetermined after solving the differential equations in η and κ. However, we can solve the DEs in such a way that these undetermined coefficient functions appear in the leading power expansion of η, i.e. in terms that correspond to j = 0 in our ansatz Eq. (3.5). The "massless" coefficients c i,0,0,0,0,0 correspond to a completely massless "version" of the MIs which is obtained by setting m 2 h and m 2 t to zero at the integrand level. These integrals are well-known and serve as an input in our calculation. Indeed, all the needed planar massless master integrals have been computed in Refs. [39,40]. 5 The non-planar massless master integrals have been taken from Refs. [42][43][44].
After fixing the massless coefficients to the known computed massless MI, we are left with undetermined coefficients c i,0,k,l,m,n . To find them, we use the ansatz in Eq. (3.5) in the z differential equation and again equate terms with matching powers of η, κ and log(κ) on both sides of the differential equation. The DEs in z are relatively simple and can be solved 5 Note that in the massless limit, the integral I P L2 1,1,1,−1,1,0,1,1,1 can be written as 1 ǫ (I P L2 1,1,1,0,1,0,1,1,1 − uI P L2 1,1,1,0,1,0,1,1,2 ) plus lower sub-topologies. The order ǫ pieces of the two planar master integrals are of weight five, but in the difference only terms with weight four survive. These terms were computed in Ref. [41]. order by order in an expansion in ǫ. Similar to the case of massless MIs, the solutions are expressed in terms of Harmonic Polylogarithms (HPLs) which form a subset of the Goncharov polylogarithms G(l 1 , · · · , l n weight n The letters that we encounter in the z DEs are very simple; the alphabet reads The first letter corresponds to branch points at s = 0 or u = 0 when η = m 2 h = 0, the second to t = 0. After solving the equations in z, we expand the solutions in ǫ keeping all the terms up to weight four 6 c i,j,k,l,m,n (z, ǫ) = r (i,j,k,l,m,n) 0 +4 r=r (i,j,k,l,m,n) 0 ǫ r c (r) i,j,k,l,m,n (z). (3.8) The powers of ǫ in the expansion are bounded below by r = −4. Typically, individual coefficient functions have higher singularities in ǫ than the expanded solution. This feature is understandable since massive internal particles screen infra-red and collinear singularities; for this reason, full results for master integrals should typically be less singular in the ǫ → 0 limit than their massless branches. After solving the DEs in z we are left with unknown integration constants that need to be determined. For the MIs in the planar topologies PL1 and PL2 we could fix many of the constants by requiring that the unphysical singularities at z = −1 cancel. We are allowed to do this since the corresponding planar diagrams do not have any cuts in the t-channel, but only in s and u. After requiring that these unphysical branch points at t = 0 vanish, all of the constants in topology PL2 become fixed. We are left with one constant in the family PL1 and six in the family NPL that we need to determine in some other way. In the next Section we will explain how we computed these constants using the Mellin-Barnes representation of the relevant integrals.
We note that in order to compute the amplitude to order O(ǫ 0 ), we are required to compute coefficient functions of some integrals to weight five and a few even to weight six. By using the DEs in η and κ, we could find many connections between contributions of weights five and six to the coefficient functions of the MIs. After substituting MIs into the amplitude, most of the unknown weight five and all of the weight six contributions cancel amongst each other. The few weight five pieces that are left, appear only in the planar families PL1 and PL2 and for these we needed to integrate the DEs in z to weight five. However the integration constants of these weight-five contributions cancel in the final result for the amplitude and therefore they did not have to be computed.  1, 1, 1, 2, 0, 1, 1, 0).
The expansions of the MIs in κ and ǫ have been, whenever possible, numerically compared with FIESTA [45] at the point m 2 h = η = 0 and an agreement was found within the integration errors of FIESTA. We include with this paper ancillary files that contain our solutions for all MIs in the form of the ansatz in Eq. (3.5), expanded in η and κ to orders that are sufficient to compute the amplitude to leading and next-to-leading power in η and κ, respectively.

Integration constants and numerical checks using Mellin-Barnes
As we already mentioned several times, by solving differential equations we determine master integrals up to integration constants that have to be determined in a different way. Many integration constants can be fixed by requiring that integrals have regular limits at certain singular points of the differential equations, for example at η → 0. However, there are seven integration constants that are left undetermined by these considerations and we have to compute them explicitly. To accomplish that, we use the Mellin-Barnes representation to calculate the relevant master integrals at certain kinematic points and then match the results to solutions of differential equations.
The Mellin-Barnes representation has been used before to compute the massless coefficient functions of some planar [42] and non-planar master integrals [46]. Since we relate the coefficient functions corresponding to higher powers in expansion in η to coefficient functions at leading power in the η expansion, all undetermined integration constants appear in the coefficient functions that can be computed by setting η to zero. In other words we have to keep the non-vanishing top quark mass 7 but we may set m 2 h = 0 in our Mellin-Barnes computation from the very beginning. The Mellin-Barnes representation is ideal for organizing the computation as an expansion in a small parameter and isolating different κ-branches since different powers of κ appear naturally after residues are computed in Mellin-Barnes integrals.
We are interested in computing an integration constant of the coefficient function that corresponds to a factor κ −1−2ǫ . We choose the kinematic point s = u = −1, t = 2, m 2 h = 0. The integral is well defined and regulated by dimensional regularization. However, the region integrals that represent the term ∝ κ −1−2ǫ are not regulated by the dimensional parameter ǫ. This can be already seen from the solution of the DE that predicts terms ∝ κ −1−2ǫ log 1,2 (κ). Such non-analytic behavior is typically cured by an introduction of analytic regulators in the context of asymptotic expansion of Feynman integrals [47,48]. The parameter δ introduced in Eq. (3.2.9) is the analytic regulator that makes expansion of all branches the integral I top 011120110 well-defined.
To proceed further, we introduce Feynman parameters and integrate over two loop momenta. The integral can be expressed as powers of the Symanzik U and F polynomials The sum inside the delta function can be chosen to be any combination of the Feynman parameters [49,50]. We have chosen the delta function as δ(1 − x 1 − x 2 ). The integration over the Feynman parameters are nontrivial but may be performed by using the method of Mellin-Barnes. Namely, we may split up terms inside the brackets by introducing Mellin-Barnes integration parameters The contour runs parallel to the imaginary axis in the complex z-plane and is chosen such that the singularities of Γ(−z) and Γ(λ+z) are to the right (left), respectively of the integration contour. After we integrate over Feynman parameters, we are left with the following Mellin-Barnes integral to perform (3.2.11) The Mellin-Barnes integrations can be performed with the help of packages collectively known as MBTools [51]. For example, the contours of the z 1..4 integrals can be systematically deformed [52] in a way that allows one to take the limit δ → 0. Indeed, because δ is an analytic regulator, we need to take δ → 0 at fixed ǫ and then deform the contour further to extract poles in ǫ and, eventually, arrive at the ǫ expansion. We note that, as follows from Eq. (3.2.11), poles in z 1 correspond to different powers of κ; for our purposes we require the pole at z 1 = −1 − 2ǫ. After extracting the κ −1−2ǫ branch and expanding the result in ǫ, we use the Barnes lemma to perform the Mellin-Barnes integrations. In most cases, these integrations are straightforward. However, we also obtain a contribution which requires to deal with the integrand that contain polygamma functions. A typical integral reads The integration over z 4 is performed using the method of residues. The integration contour runs along the imaginary axis with Re(z 4 ) small and negative. We may close the integration contour to the left as the integrand will vanish fast enough along the half circle in the left complex plane with infinite radius. By Cauchy's theorem we pick up the ladder of residues in z 4 in the left complex plane, Re(z 4 ) < 0. Application of this procedure leads to the following representation 1 n 4 (n + 1) 2 π 2 n 2 + 6 (n + 1) 2 + 3n 2 3n 2 − 3 n 2 + n 2 ψ (1) (n + 1) − 1 ψ (0) (n + 1) +3n(n + 1) 2 ψ (0) (n + 1) 2 − n(n + 1)((n − 3)ψ (1) (n + 1) + 2n(n + 1)ψ (2) (n + 1)) .
These sums can be performed with e.g. the XSummer [53]. The final result for the O(κ −1−2ǫ ) branch of the I top 011120110 at the kinematic point s = u = −1, t = 2, m 2 h = 0 reads (3.2.13) We then match the solution of the differential equation to this result and determine the integration constant.
In addition to determination of constants, we also used the Mellin-Barnes representation for numerical checks of our solutions for master integrals. Namely for non-planar MIs that, in the κ → 0 limit develop power-like singularities with half-integer exponents, we were unable to use FIESTA for numerical checks. For such integrals we compared individual coefficient functions in κ with the Mellin-Barnes representation and found perfect agreement for all of them. In particular, the massless contributions as well as other coefficient functions that are completely fixed by the DE, all agree with the Mellin-Barnes result. Note that this is a nontrivial check on both the solution for the differential equation and the Mellin-Barnes representation that we used to extract integration constants for certain branches. One example of the coefficient functions that are completely fixed by the DEs are those corresponding to the κ −1−2ǫ log(κ) and κ −1−2ǫ log 2 (κ) terms that appear in our solution of the above integral I top 011120110 , which we have checked to agree exactly with the corresponding logarithms in κ in Eq. (3.2.13) for the chosen kinematic point.

Helicity amplitudes
Once the master integrals are computed, we use them to derive the form factors and calculate the analytic expressions for helicity amplitudes. We define positive and negative helicity spinors for massless external quarks and gluons in the standard way (see e.g. [54]) Here q is an arbitrary light-like reference vector. For our computation, the reference vectors are fixed by gauge conditions outlined in Eq.(2.10).
The helicity amplitudes are defined as Eight helicity configurations are needed to describe the H → ggg amplitude. However, only two of them are independent since the other six may be related to them by the use of charge and parity conjugation. For the H → qqg amplitude there are four possible helicity configurations in total, since QCD interactions cannot change the helicity of the massless quarks and therefore the helicity of the outgoing quark must be opposite to that of the outgoing anti-quark in Eq. (4.4). We have chosen to treat the following amplitudes as independent

6)
A q −++ (s, t, u, m t ) = 1 √ 2 [23] 2 [12] s Ω q −++ (s, t, u, m t ). (4.7) The amplitudes are dimensionless and the helicity coefficients Ω i (s, t, u, m t ) have a mass dimension one. We may obtain the other helicity assignments for the amplitude by complex conjugation and by permuting the external legs as follows A g ++− (p 1 , p 2 , p 3 ) = A g +−+ (p 1 , p 3 , p 2 ) , (4.8) The complex conjugation should only be applied to spinor-helicity structures and not to the helicity coefficients Ω i (s, t, u, m t ). The helicity coefficients can be expressed in terms of the form factors introduced in Eq.(2.9) as follows (4.12) We expand the helicity coefficients in the strong coupling constant and extract an overall coefficient m 2 t /v in order to have dimensionless one-and two-loop helicity coefficients Once the form factors have been renormalized and IR-subtracted, the resulting helicity coefficients will also be finite as seen from Eq. (4.12). We are interested in a kinematic region where all Mandelstam variables are much larger than the top mass m t . Therefore we prefer to define the amplitude in terms of a strong coupling constant that runs with N f + 1 active flavors. The relation between the coupling constants defined in the two schemes reads (4.14) This change in the strong coupling constant leaves the one-loop coefficients unchanged, but the two-loop finite remainder of the helicity amplitude changes as follows according to Eq. (4.14) The helicity coefficients Ω correspond to using a strong coupling constant α (N f +1) s (µ 2 ) that evolves with N f + 1 active flavors.
Unfortunately, analytic results for helicity amplitudes are too long to be presented here. Instead, we provide ancillary files that contain finite remainders of the relevant helicity amplitudes Ω defined in Catani's subtraction scheme Eq. (2.17) with the submission of this paper.

Analytic continuation
Our goal is to compute the two-loop amplitudes that are needed to describe production of the Higgs boson with high transverse momentum at the LHC. The relevant production channels are gg → Hg, qq → Hg, qg → Hq andqg → Hq. Amplitudes for these processes are obtained from amplitudes for H → ggg and H → qqg processes that we have computed and crossing some final state particles to initial states. Our H → ggg and H → qqg amplitudes have been computed in the region t > 0, s, u < 0, while the physical scattering processes are defined in the kinematic regions where the invariant mass of the two initial partons is positive instead. For this reason we are interested in computing the amplitudes in the regions s > 0, t, u < 0 and u > 0, s, t < 0 as well, which we will refer to as (2a) + and (4a) + respectively. The amplitude in these regions can be found by analytically continuing our result from the region t > 0, s, u < 0 which we refer to as (3a) + . The three scattering regions are defined in terms of the Mandelstam invariants as (2a) + : s > 0 , t, u < 0 , (4.1.1) (3a) + : t > 0 , s, u < 0 , (4.1.2) (4a) + : u > 0 , s, t < 0 .  The method to perform the analytic continuation from the region (3a) + , where our computation has been performed, to the other two regions was explained in Ref. [43] and we refer to this paper for details. The spinor products are left unchanged during the analytic continuation but Harmonic Polylogarithms may receive imaginary parts when continued to regions (2a) + and (4a) + . We introduce the variable u j for the three scattering regions Our helicity amplitudes are expressed in terms of the new variables 0 ≤ u j ≤ 1 in the three corresponding regions. In this way the imaginary part of the amplitudes is explicit and all the HPL that appear in the results are real-valued with the alphabet {0, 1} in each of the scattering regions. The Harmonic Polylogarithms can be numerical evaluated with the Mathematica package HPL [55] or the Fortran code CHAPLIN [56]. The helicity amplitudes Ω q +++ , Ω g +−+ Ω q −++ in all three scattering regions are provided in the ancillary file together with the submission of this paper.

Conclusions
We computed the two-loop helicity amplitudes that are needed to describe production of the Higgs boson with large transverse momentum at the LHC. The Higgs boson interaction with gluons and massless quarks is mediated by loops of massive top quarks. However, the top quark mass is considered to be small compared to Higgs bosons transverse momentum. Clearly, in this kinematic regime the Higgs boson mass is also small compared to its transverse momentum and we effectively neglect it in our computation.
Although the dependence of the scattering amplitudes on the Higgs boson mass is simple and can be obtained by a simple Taylor expansion, the expansion of the amplitudes in the top quark mass contains non-analytic terms O(ln(m 2 t /p 2 ⊥ )) and is, therefore, non-trivial. We construct the expansion of the amplitudes using differential equations for master integrals that allow us to obtain both analytic and non-analytic terms in an expansion in a controlled way. Our final results for the amplitudes are expanded to leading power in the Higgs boson mass which, essentially, corresponds to setting the Higgs boson mass to zero, and to next-to-leading power in the top quark mass squared. We expect that the two-loop amplitudes computed in this paper will allow for a robust estimate of the number of Higgs bosons that are expected to be produced at the LHC with very large transverse momentum, and the comparison of this prediction with the experimental result [12].