The one-loop amplitudes for Higgs + 4 partons with full mass effects

We present compact analytic formulae for the one-loop amplitudes for Higgs + 4 parton scattering, 0 → ggggh, 0 →q¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{q} $$\end{document}qggh and 0 →q¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{q} $$\end{document}qq¯′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\overline{q}}^{\prime } $$\end{document}q′h, mediated by a loop of massive coloured quarks. We exploit the correspondence with a theory in which a massive coloured scalar circulates in the loop to avoid a proliferation in the number of terms in the result. In addition, we use momentum twistors and high precision numerical evaluations to simplify the expressions. The analytic results in this paper, in terms of spinor products, allow construction of an efficient numerical program to calculate the amplitude.


Introduction
At the Large Hadron Collider (LHC) the primary mechanism for producing and detecting Higgs bosons is the process gg → h. This process is mediated, in the Standard Model, by a loop of massive coloured fermions. Since the Yukawa coupling is proportional to the fermion mass the predominant contribution is the result of the coupling of the top quark to the Higgs boson. In the limit in which only a very heavy top quark contributes, the corresponding amplitude is independent of the top quark mass; this gives rise to an effective field theory (EFT) in which the loop of heavy top quarks is replaced by an effective Lagrangian, where g s is the strong coupling constant, v is the vacuum expectation value of the Higgs field, G µν is QCD field strength, and h is the Higgs boson field. The EFT in eq. (1.1)

JHEP05(2020)079
has been used to compute higher-order corrections to the inclusive cross-section -most recently up to next-to-next-to-next-to-leading order [1,2] -as well as rates for the production of Higgs bosons in association with up to three additional jets up to next-to-leading order [3,4]. The effective field theory description is expected to break down when, for example, the transverse momentum of produced gluons is of order of the top quark mass. This breakdown has most recently been investigated at NLO in ref. [5]. This kinematic regime is beginning to be explored at the LHC [6] and can give important information about the mediators in the loop that couple to the Higgs. For such configurations it is therefore important to make use of a superior calculation in which the full dependence on the top quark mass is retained. Such a calculation also allows a direct quantification of the breakdown of the EFT approach.
Analytic results for the Higgs+3 parton amplitude in the full theory have been known for a long time [7,8]. Corresponding results for Higgs+4 parton amplitudes have been obtained in refs. [9,10], although in both cases expressions for at least some of the amplitudes were too long to report. In addition there are several automatic procedures than can provide numerical results for one-loop amplitudes [11][12][13][14]. The aim of this paper is to present compact amplitudes for all contributing processes, 0 →qqq q h , (1.4) retaining all mass effects. Compact analytic results for the 0 → ggggh case when all the gluons have positive helicity have been published in ref. [15]. Although our result is therefore not new per se, it is the first time that a compact publishable analytic result has been obtained for all gluon helicities. A calculation with compact analytic formulae allows examination of the structure of the amplitude for all values of the fermion mass. It also has the potential to lead to faster and more stable numerical evaluation of the amplitude. This would be a boon to calculations requiring this amplitude in all regions of phase space, such as recent NLO predictions for Higgs boson plus 1-jet production in the full theory [5] and at large transverse momentum [16,17]. Although the results are quite compact, given the number of integral coefficients, this paper is not easy to read. However we believe that it is detailed enough that readers wishing to implement this amplitude in a numerical program will find enough information to do so in our paper.
2 Structure of the calculation

Definition of colour amplitudes
The amplitude for the production of a Higgs boson and n gluons can be expressed in colour-ordered sub-amplitudes as follows:

JHEP05(2020)079
where the sum with the prime, {1,2,...,n} , is over all (n − 1)! non-cyclic permutations of 1, 2, . . . , n and the t matrices are the SU(3) matrices in the fundamental representation normalized such that, tr(t a t b ) = δ ab . (2.2) m is the mass of the quark circulating in the loop. Because of Bose symmetry it is sufficient to calculate one permutation, and the other colour sub-amplitudes can be obtained by exchange. For the particular case at hand with four gluons eq. (2.1) becomes, Squaring the amplitude eq. (2.3) for a fixed helicity configuration and summing over colours we find colours |H gggg where N is the dimensionality of the SU(N ) colour group, i.e. N = 3, and the labels for the helicity configuration (as explicitly shown in eq. (2.3)) have been suppressed. The amplitude for the production of a Higgs boson, an antiquark, quark and two gluons is similarly decomposed into colour-ordered amplitudes as follows, The colour structure δ c 3 c 4 δ j 2 j 1 /N is also present in individual diagrams but makes no net contribution to the one-loop amplitude. In this paper we will give results for the colour-ordered amplitude H 34 4 . It is straightforward to obtain H 43 4 from this through the parity operation (complex conjugation) and permutation of momentum labels. Squaring and summing over colours yields, 6) where the labelling of the helicity configuration shown in eq. (2.5) has again been suppressed.
The four-quark amplitude takes the form,

JHEP05(2020)079
where the helicities of the quarks are fixed by those of the antiquarks.Performing the sum over colours we then have, for the case in which the quark lines have different flavours. For the case of identical quarks we first introduce, The sum over the colours for the identical case is then, where, as indicated, the term on the second line only contributes when the quarks have the same helicity.

Decomposition into scalar integrals
The colour-ordered sub-amplitudes can be expressed in terms of scalar integrals. For the 0 → ggggh sub-amplitude we have, The scalar bubble (B 0 ), triangle (C 0 ), box (D 0 ) and pentagon (E 0 ) integrals, and the constant r Γ , are defined in appendix A.μ is an arbitrary mass scale, and r are the rational terms. The rank of a Feynman integral is defined to be the number of powers of the loop momentum in the numerator. A scalar Feynman integral has no powers of the loop momentum in the numerator, and is hence of rank zero. All scalar integrals are well known and readily evaluated using existing libraries [18][19][20]. The sums in the above equation scan over groupings of external gluons. Thus, for example, the sum for the scalar triangle integrals will contain a term c 1×234 which multiplies the scalar triangle integral C 0 (p 1 , p 234 ; m) where p 234 = p 2 + p 3 + p 4 . The reduction in eq. (2.11) is written in n dimensions, although at the end the amplitude is finite. The individual bubble integrals contain ultra-violet singularities that are regulated using dimensional regularization.

JHEP05(2020)079
In four dimensions the pentagon integral can be reduced to a sum of the five box integrals obtained by removing one propagator [21][22][23], Explicit forms for the pentagon reduction coefficients, C  The factor |S 1×2×3×4 | is the determinant of the matrix, (2.14) As a result of eq. (2.12), in four dimensions the integral basis given in eq. (2.11) is overcomplete. The full amplitude can be described by the box, triangle and bubble integrals alone (+ rational terms). This is the specific choice made in this paper, but the other choice to keep the redundant basis of eq. (2.11) is also perfectly viable. In this paper we will work in a basis without pentagon integrals, but the box coefficients will in part display vestiges of their pentagon origin, through effective pentagon coefficients and the presence of the pentagon-to-box reduction coefficients, C 1×2×3×4 . This will be explained in detail in section 5. Our decomposition of the sub-amplitudes is thus, (2.15) which also applies for the 0 →qqggh sub-amplitude H 34 4 since it contains no pentagon diagrams in the first place.

Unitarity methods
The modern treatment of one-loop amplitudes containing massive particles was pioneered 25 years ago in ref. [24]. Since that early paper a whole set of tools and methods have been invented to deal with one-loop amplitudes (for an introduction and comprehensive review see refs. [25] and [26] respectively). We shall apply many of them in order to arrive at the simplest form for the Higgs + 4 parton amplitude. This paper will present compact expressions for the coefficients in the four dimensional version of eq. (2.11) where the scalar pentagon integral has been expressed as a sum of box integrals. The coefficients will be expressed in terms of spinor products. Our notation for spinor products is reported in appendix B.
As we shall see below, there is an intimate connection between the full one-loop calculation with a massive fermion and a suitably normalized calculation performed with the Higgs boson coupling to four partons via a loop of colour-triplet, massive scalar particles. The latter calculation with scalar intermediaries has two advantages. First, the scalar calculation is completely free of Dirac algebra, allowing more compact expressions to be maintained throughout the calculation. This is useful if one can show the identity of the coefficients of scalar integrals between the scalar and the fermionic theories. Second, the scalar calculation, unlike the fermionic calculation, can be performed in the m → 0 limit. If it can be shown that: 1. the result for a particular coefficient in the scalar theory is identical to the result in the fermionic theory, 2. that particular coefficient is also independent of the mass, the value of the m → 0 limit is established. In order to perform the reduction to scalar integrals indicated in eq. (2.11) we use unitarity techniques to isolate the contribution of boxes [27], triangles [28] and bubbles [29][30][31]. Since bubble integrals satisfy both criteria enumerated above, their coefficients are most easily calculated using a massless internal scalar loop.
Integrals that do not give rise to rational terms are said to be cut-constructible. In general, x-point integrals are cut constructible in four dimensions if the rank r satisfies r < max[(x − 1, 2)]. Thus rank-3 pentagons, rank-2 boxes, rank-1 triangles and bubbles are cut constructible. Integrals that are not cut-constructible give rise to the rational terms (r) in eq. (2.11). In our case the rational terms can instead be obtained by using already-computed results for the mass-dependent coefficients of triangle integrals [32].

Simplification techniques
We now briefly describe two further techniques that are useful to help simplify the results obtained using unitarity methods. Both methods exploit the fact that the kinematics of our process can be expressed as massless 6-point kinematics by decomposing the momentum of the Higgs boson into two light-like momenta, which for definiteness we call p 5 , p 6 .
Reduction of the analytic forms to simpler expressions is aided by the use of momentum twistors [33][34][35][36]. In this formalism each particle is described by a 4-component momentum JHEP05(2020)079 twistor Z(λ, µ), where λ is the usual two-component holomorphic Weyl spinor (with i j = λ α λ α ) and µ is a two-component object related to dual momentum coordinates [33]. Antiholomorphic spinors (λ i , with [i j] =λαλα) are obtained from these via the identity, To describe an n-particle scattering amplitude there are thus 4n momentum twistor components, of which only 3n − 10 are independent due to a U(1) symmetry for each particle and overall Poincaré symmetry. We thus need 8 momentum-twistor variables (x 1 . . . x 8 ) to describe our 6-point kinematics, which we choose to parametrize as, The spinors involving our four massless partons are then given by, where we see explicitly that the variables x 3 and x 4 involve momenta p 5 and p 6 and effectively decouple in our case [36]. In order to use this parametrization we first need to remove the overall phase of the coefficient corresponding to the helicities of the external gluons, for example by multiplying by 1 2 2 3 4 2 for the all-plus amplitude. The advantage of this approach is that the amplitude can now be simplified using straightforward algebra, without needing to account for momentum conservation and Schouten identities to manipulate spinor strings. In this way overall factors can easily be identified and the true denominator structure of the coefficients established. The second method we adopt is to use high precision floating-point arithmetic to simplify our analytic expressions [37]. The study of singular and doubly singular limits in JHEP05(2020)079 complex phase space allows us to explore the singularity structure of the integral coefficients. The integral coefficients can then be reconstructed by solving linear systems for the rational numerical coefficients of generic spinor trial functions. The helicities of the gluons impose constraints on the structure of the trial functions. This method is particularly useful when unitarity techniques result in lengthy expressions that are hard to treat using twistor variables, such as in the case of some triangle and bubble coefficients. It is also useful to bypass the algebra involved in removing artefacts of loop-momentum parametrizations, such as square roots and massless projections of non-lightlike external momenta [28,38]. This paper is the first instance where this method has been applied in the presence of massive particles.

Higgs boson production mediated by a coloured scalar
Consider a complex scalar field φ in the triplet representation of colour SU(3) coupled to a gluon field and to the Higgs boson, h. The part of the QCD Lagrangian involving the field φ is The partial correspondence with the fermion theory emerges when setting λ = −4m 2 /v. We will calculate colour-ordered sub-amplitudes for the production of a Higgs boson and four gluons mediated by a scalar loop. For the Higgs + 4 gluon amplitude this is, and the colour amplitudes have the decomposition, Thus the tilde indicates that we are referring to an amplitude mediated by a scalar field. For the case of the triangle coefficients we divide the coefficient into two pieces, to separate JHEP05(2020)079 the mass dependence. Thus for both the fermion-and scalar-mediated loops we have, In eq. (3.2) we have chosen the normalization so that in some cases the coefficients c i×j = c i×j and in addition in all cases c (2) i×j =c (2) i×j , b i =b i and r =r. Thus we can perform certain parts of the calculations in the scalar theory. We further have thatb i are independent of the mass m, so that they can be calculated in the massless scalar theory. In addition, for the case at hand the rational terms are fully fixed byc (2) i×j , where the sum runs over all non-zeroc (2) i×j for a particular helicity.

Relationship of the fermion theory to the scalar theory
In order to elucidate the relationship between the fermion and scalar theories [39] it is instructive to review the steps previously used to demonstrate their similarity using the second order formalism [40]. Following ref. [40] we define the quantity A to represent the combination of the numerator part of a fermion propagator with momentum and a gluon-quark-antiquark vertex at which momentum p 1 flows out along the gluon line (3.7) The first term on the right hand side of eq. (3.7) already resembles the vertex for a gluon (on-shell or off-shell) interacting with a scalar field. Now consider the integrand of a triangle diagram for a Higgs boson coupled to two gluons as shown in figure (1a), where our notation for the loop momenta and propagator factors D( ) is given in appendix A. The minus sign is included because of the fermion loop and an overall factor of JHEP05(2020)079 1/m has been included for convenience. By using the decomposition in eq. (3.7) twice, the expression in eq. (3.8) can be brought into the form, where B is a four-by-four matrix-valued function. It contains a convection term (as expected for a scalar field) and a spin term, Note that the spin term is independent of the loop momentum. The result for the companion triangle as shown in figure (1b) is Adding both diagrams, dropping vanishing terms and exploiting the cyclicity of the trace, the integrand appearing in the full amplitude is, .
If we drop the terms involving the commutators of gamma matrices, the fermionic loop (after removing an overall factor of m) can be written as the effect of a (suitably normalized) scalar triangle, with 3-point and 4-point (seagull) vertices. The full fermionic theory requires the inclusion of the additional spin flip terms, given by commutators. These additional terms do not involve the loop momentum and are thus of lower rank. Note also that there is no explicit mass dependence in eq. (3.12). Dependence on the mass m will be generated by the reduction to scalar integrals. In contrast to the full fermionic theory, we may also consider the scalar theory in the massless case. Iterating this argument for a larger number of gluons, it can be shown that separation of the full fermionic theory into a suitably normalized scalar theory, plus spin terms involving gamma matrix commutators of rank lower by two powers of , continues to hold. The scalar theory is obtained by dropping all of the spin terms, cf. eq. (3.10). Thus the full amplitude can be written as the sum of the scalar theory and a correction of lower rank, ∆F , Fermion theory = Scalar theory + ∆F (3.13) JHEP05(2020)079 The difference between the scalar theory and the full fermionic theory ∆F is of lower rank. Although the scalar theory contains rank-4 pentagons, rank-3 boxes, and rank-2 triangles, ∆F contains only rank-2 pentagons, rank-1 boxes and rank-0 triangles. This has several important consequences: 1. ∆F is cut constructible.
2. ∆F gives no contribution to bubble integral coefficients. Bubble integral coefficients can thus be calculated in the scalar theory.
3. ∆F gives no contribution to the m 2 contributions to triangle coefficients, c i×j .
4. ∆F gives no contribution to certain m 0 triangle coefficients, c i×j .
Exploiting these facts provides the following identities for the ggggh case. (3.14) i.e. all triangles that do not have the Higgs boson as an external leg can be calculated fully in the scalar theory. In appendix C we reproduce several results for tree graphs involving massive scalars, Higgs bosons, gluons and a quark-antiquark pair. These results are useful in applying unitarity to calculate loop diagrams.
4 Coefficients for H 1234 4 (g + , g + , g + , g + ; h) The reduction of a scalar pentagon integral into a sum of box integrals shown in eq. (2.12) applies in four dimensions. Therefore an extraction of the scalar pentagon integral coefficient must be performed by making use of unitarity methods in d dimensions. To this end the loop momentum is expressed most generally as, where represents the excursion beyond four dimensions, with 2 = −µ 2 . Putting the propagators on shell determines α, β, γ, δ and µ 2 . Parametrizing the amplitude at hand with the decomposition in eq. (4.1) we thus find the pentagon coefficient, where the trace functions are defined by, and γ R/L = (1 ± γ 5 )/2. The identities on the far right of eq. (4.3) holds only for lightlike p i . The value of µ 2 is fixed by the constraint 2 = m 2 (implying m 2 + µ 2 = −γδ s 12 ) but we have left the coefficient in eq. (4.2) in this form in order to emphasise the d-dimensional nature of the result. We choose to write our amplitude in terms of an effective pentagon coefficientê that corresponds to the four-dimensional limit µ 2 → 0, This leads to a very compact form for the complete amplitude [15], (4.5) Although this expression, that includes the scalar pentagon integral, is very simple, we do not follow this approach for the other helicity choices in the following sections. Instead, since in the end our aim is to produce a numerical code to calculate this amplitude, we feel that the structure is more straight-forward working only in terms of boxes, triangles and bubble coefficients. Adopting this approach also for this amplitude means that the minimal set of coefficients that must be specified corresponds to the ones shown in the first and third columns of table 1. Note that coefficients of integrals that could in principle appear but that are not specified in this table (and in subsequent tables in later sections) vanish. There is some ambiguity in the naming convention for the integral coefficients. For a given colour-ordered amplitude the external legs will appear in cyclic or anti-cyclic order, cf. eq. (2.3). For example, d 1×2×34 in table 1 could equally well be written as d 43×2×1 . Our convention is that the vertex containing more than one gluon, should appear last in the name of the coefficient. Where the compound vertex is in the centre, we have chosen the cyclic ordering. Table 1. Minimal set of integral coefficients for 1 +  In presenting these results we have adopted the notation defined in eq. (3.4) to separately quote the mass-independent and mass-dependent parts.

Rational terms
1×234 (4 + , 1 + , 2 + , 3 + ) = 4 s 1234 1 2 2 3 3 4 4 1 (4.11) 5 Coefficients for H 1234 4 (g + , g + , g + , g − ; h) Following the procedure outlined in section 4 to obtain the pentagon coefficients yields, for this particular helicity combination, Note that the last identity, eq. (5.2), only applies for the case of lightlike momenta. At this point we could follow the same strategy as in the previous section and take the limit µ 2 → 0 to obtain an effective pentagon coefficientê. However the appearance of the factor tr 5 {1 2 3 4} in the denominator of eq. (5.1) is unpalatable since it is an unphysical singularity. In the amplitude its presence is compensated by corresponding factors in box coefficients, and separating contributions in this way can lead to a loss of numerical precision.
As an alternative solution, we also choose to modify the coefficient itself in such a way that these factors are eliminated. We do so by noting that the denominator can be expressed as, where G has already been introduced in eq. (2.14). In fact, rearranging that equation and making use of spinor notation we have, We may use this equation to eliminate the factor of tr 5 {1 2 3 4} 2 in the pentagon coefficient indicated in eq. (5.1). The additional term that is introduced contains a factor of |S 1×2×3×4 |; this neatly cancels the denominator factor involved when reducing the pentagon integral to boxes (cf. eq. (2.13)) such that factors of 1/tr 5 {1 2 3 4} 2 can be explicitly absorbed into the box coefficients and cancelled.
In this way we arrive at effective pentagon coefficients, The correspondence betweenê {1 + ×2 + ×3 + ×4 − } defined here, and e {1 + ×2 + ×3 + ×4 − } given in eq. (5.1) is clear from the first term on the right-hand side of eq. (5.4). We will see that the box coefficients take a particularly simple form when written this way. The minimal set of remaining integral coefficients to determine this amplitude is shown in the first and third columns of  Table 2. Minimal set of integral coefficients for 1 + determined by using symmetry properties of the amplitude and relabelling momenta. These are mostly straightforward except for the coefficient b 341 where we have found the relation,

Triangles
For the case of the triangle coefficients we divide the coefficient into two pieces, to separate the mass dependence, see eq. (3.4). For many coefficients the c (2) i×j term is equal to zero. For these cases, the full result is given by the m 0 term and we shall omit the superscript.

b 1234
Since the full amplitude is finite in four dimensions, one of the coefficients is uniquely determined in terms of the remainder. We thus have, where we have suppressed momentum and helicity labels on the right-hand side for brevity.

Rational terms
Coefficients for H 1234 4 (g + , g − , g + , g − ; h) For this helicity combination the coefficients of the scalar pentagon integrals contain a factor of 1/tr 5 {1 2 3 4} 4 and we must modify the pentagon integral coefficients in a similar fashion as described for the + + +− configuration in section 5. All coefficients can then be written in terms of a single function,

JHEP05(2020)079
Coefficient Related coefficients Coefficient Related coefficients Table 3. Minimal set of integral coefficients for 1 + The minimal set of integral coefficients that must be calculated for the colour ordering H 1234 4 is shown in the first and third columns of table 3, for example the bubble coefficients are given by: (6.5) The calculation of the coefficients of other colour orderings requires the use of + + − − functions which are given in the next section.
6.1 Boxes  Note that here the symmetrization only applies to the terms in square brackets. JHEP05(2020)079 where ∆ 3 is given by eq. (B.2).

c
where ∆ 3 is given by eq. (B.2). Similarly to eq. (5.32) we have, suppressing momentum and helicity labels on the right-hand side for brevity.

JHEP05(2020)079
Coefficient Related coefficients Coefficient Related coefficients Table 4. Minimal set of integral coefficients for 1 + 7 Coefficients for H 1234 In this case, as for the all-plus helicity amplitude, there are no factors of 1/tr 5 {1 2 3 4} 2 in the pentagon integral coefficients. Therefore the effective pentagon integral coefficients simply correspond to the µ 2 → 0 limit, as in the + + ++ case. We thus have, The minimal set of integral coefficients that must be computed in this case is shown in the first and third columns of table 4. For the colour ordering H 1234 4 the complete set of related coefficients is given in table 4, for example bubble coefficients are given by: (7.4) The calculation of the coefficients of other colour orderings requires the use of + − + − functions which are given in the previous section.   The full coefficient for this scalar integral is defined in terms of the coefficient with a scalar running in the loop,c

Boxes
In turn the coefficient for a scalar loop is given by, 12×34 is given by eq. (6.13) with the appropriate permutation of arguments.

Rational terms
together with the related coefficients that can be obtained from the base set.
8 Coefficients for H 34 4 (q + , q − , g + , g + ; h) The coefficients that must be computed for this amplitude are shown in table 5.

b 1234
The final bubble coefficient is given by the relation, where momentum and helicity labels have been suppressed on the right-hand side.

Rational terms
Coefficients for H 34 4 (q + , q − , g − , g + ; h) The coefficients that must be computed for this amplitude are shown in table 5.  The final bubble coefficient is given by the relation,

Boxes
where momentum and helicity labels have been suppressed on the right-hand side. 9.13) 10 Coefficients for H 34 4 (q + , q − , g + , g − ; h)

Rational terms
Most of the coefficients for this amplitude can be easily obtained from those for H 34 4 (q + , q − , g − , g + ) by performing the following operation: However, for some coefficients, this procedure effectively changes the colour ordering of the gluons in the sub-amplitude. For this reason it is necessary to specify the four coefficients shown in table 5. Results for the base set of coefficients are given explicitly here.  10.5) 11 Amplitude for 0 →qqqqh

Triangle
The amplitude for the ggh process with off-shell gluons with momenta k 1 , k 2 has been given in ref. [9]. Thus our result is exactly given by ref. [9] and is only included here for completeness.
Our results for the amplitudes were checked using an in-house implementation of the D-dimensional unitarity method [44] and also against a previous unitarity-based calculation [10]. Complete agreement was found at the amplitude level. Our results for the squared matrix elements are also in full agreement with those obtained using the code OpenLoops 2 [14]. A comparison of the evaluation time of squared matrix elements against both the previous code implemented in MCFM [45][46][47] and OpenLoops 2 indicates a speed-up by at least an order of magnitude over previously-available results. Our results have been implemented in version 9.1 of the code MCFM, 2 that includes a calculation of the full matrix elements for the three partonic processes under consideration. Results are given in the subdirectory MCFM-9.1/src/ggHgg mass of that release and the result for a particular equation can be found by searching for the text "Implementation"; every file that gives the result for an integral coefficient contains a comment of the form: ! Implementation of eq.~(x.xx) from arXiv:2002.04018 v2 The results of this paper will be useful for improving calculations of the h+jet process at NLO in the full theory. The analytic forms presented here will also be useful in their own right, for understanding the structure of gauge-theory amplitudes. In particular the simplification of our results due to the choice of effective pentagon integral coefficients may have deeper origins that remain to be explored.

B Spinor algebra
All results are presented using the standard notation for the kinematic invariants of the process, s ij = (p i + p j ) 2 , s ijk = (p i + p j + p k ) 2 , s ijkl = (p i + p j + p k + p l ) 2 , (B.1) and the Gram determinant, We express the amplitudes in terms of spinor products defined as,

C.3 Three gluons
The spurious poles in the original BCFW form of the amplitudes [48] can be eliminated and the amplitudes rewritten in the following form [49],