Two-loop helicity amplitudes for H + jet production to higher orders in the dimensional regulator

: In view of the forthcoming High-Luminosity phase of the LHC, next-to-next-to-next-to-leading (N 3 LO) calculations for the most phenomenologically relevant processes become necessary. In this work, we take the ﬁrst step towards this goal for H + jet production by computing the one-and two-loop helicity amplitudes for the two contributing processes, H → ggg , H → q ¯ qg , in an eﬀective theory with inﬁnite top quark mass, to higher orders in the dimensional regulator. We decompose the amplitude in scalar form factors related to the helicity amplitudes and in a new basis of tensorial structures. The form factors receive contributions from Feynman integrals which were reduced to a novel canonical basis of master integrals. We derive and solve a set of diﬀerential equations for these integrals in terms of Multiple Polylogarithms (MPLs) of two variables up to transcendental weight six.


Introduction
A little over a decade ago, the Higgs boson was discovered after analysing the data collected during Run I at the Large Hadron Collider (LHC) at CERN [1,2].The discovery was achieved while the collider was running at reduced centre-of-mass energies of 7 and 8 TeV and with only a small fraction of the total dataset which will be accumulated during its entire runtime.Indeed, it is expected that the forthcoming High-Luminosity phase of the LHC (HL-LHC) will yield a dataset corresponding to 3 ab −1 of integrated luminosity for pp collisions at 14 TeV [3].
Since its discovery, the Higgs boson has been at the centre of the experimental effort at the LHC [4].Studying its properties improves our understanding of electroweak symmetry breaking (EWSB), the mechanism which is believed to be responsible for the generation of the masses of fermions and weak gauge bosons.The dominant production channel for the Higgs boson at the LHC is gluon fusion.In the Standard Model, the Higgs coupling to two gluons is mainly mediated through a loop of top quarks, making it a loop-induced process already at Leading Order (LO).For this reason, computing higher order perturbative corrections to Higgs production in the full theory quickly becomes prohibitive.
It was realized long ago that radiative corrections can increase the LO Higgs crosssection in gluon fusion by as much as O(100%) [5,6] and describing its production in hadron collisions thus requires higher-order calculations.Indeed, quite recently the computation of the fully inclusive Higgs cross-section with full dependence on the top quark mass has been pushed to next-to-next-to-leading-order (NNLO) [7] using numerical techniques to handle the required two-and three-loop scattering amplitudes.
A useful alternative to performing calculations with full dependence on the top mass is to work in the heavy top quark mass limit M t → ∞.Under the assumption that the top quark is the largest scale involved in the calculation, one can integrate out the top mass and formulate an effective Lagrangian for the Hgg coupling [8][9][10].In this description, the top-quark loop mediating the Hgg interaction shrinks to a point and calculations start at tree level, involving only massless partons.In this limit, inclusive [11][12][13] as well as fully differential [14] predictions for Higgs boson production via gluon fusion are known up to next-to-next-to-next-to-leading-order (N 3 LO).
One of the most promising observables to study the EWSB mechanism is the Higgs transverse momentum, see for example [15].To remain differential in external radiation, one must study the production of a Higgs boson in association with (at least) one resolved jet.For this process, results in the heavy top quark limit are currently known up to NNLO [16][17][18][19][20], and the residual theoretical uncertainty can be estimated at around O(5%).The heavy top quark limit approximation is valid for transverse momentum of the Higgs which is lower than two times the top quark mass, p T < 2 m t .For higher p T , the heavy quark loop is resolved and finite mass corrections are needed.Results with finite top mass at NLO were first obtained for the very high transverse momentum kinematic region, p T ≫ 2m t , [21], and numerically for general kinematics [22,23].More recently, following the calculation of the relevant two-loop master integrals [24,25], the full NLO analytical calculation has been completed [26].
Taking into account the wealth of currently available data, as well as the expected high-luminosity phase, N 3 LO calculations will become essential [27] in order to perform phenomenological studies at the 1% level at the LHC.A key ingredient to extend the current calculations to N 3 LO are the four-point amplitudes for the production of a Higgs boson and a parton in parton-parton collisions.Working in the effective theory described above, one needs to compute tree level, one-loop, two-loop and three-loop amplitudes for the scattering of three massless partons and one massive scalar.On top of their phenomenological importance, the structure of these amplitudes to higher loops is also of formal interest and has been the subject of thorough investigation.In this context it is worth noticing that contrary to the naive expectation, up to two loops, the finite remainder amplitudes for the decay of a Higgs boson to three gluons have been shown to be expressible in terms of just classical polylogarithms [28].
Starting at one loop, amplitudes exhibit singularities which can be regulated in the framework of dimensional regularization.In D = 4−2ϵ dimensions, the loop amplitudes are computed as a Laurent expansion in ϵ.An N 3 LO computation requires one-loop amplitudes up to O(ϵ 4 ), two-loop amplitudes up to O(ϵ 2 ) and three-loop amplitudes up to O(ϵ 0 ).Our goal in this paper is to provide the first ingredient for such a calculation, namely the twoloop amplitudes up to order ϵ 2 .These amplitudes were previously obtained in [29] to order ϵ 0 , providing a key ingredient to the calculation of NNLO QCD corrections to Higgs+jet production [16][17][18][19][20] and the Higgs boson transverse momentum distribution [19,30].
While our overall approach is relatively standard, we include multiple new elements which help us organize the calculation more efficiently in view of a subsequent extension to three loops.First of all, we cast the relevant amplitudes into a compact tensorial basis and construct new helicity projectors to extract the corresponding helicity amplitudes directly from the Feynman diagrams [31,32].We then employ standard integration-by-parts identities [33,34] to express these amplitudes in terms of so-called master integrals, and evaluate them using the differential equation method [35][36][37][38].At two loops, the master integrals were computed up to transcendental weight four more than two decades ago [39,40].Since then, substantial advances were made in the understanding of mathematical structures underlying Feynman integrals and their associated differential equations.
Indeed, about a decade ago it was realized that a special class of Feynman integrals, dubbed local integrals [41,42], plays a crucial role in representing scattering amplitudes, in particular in the case of (planar) N = 4 Super Yang-Mills theory (SYM).These integrals only feature singularities of the logarithmic type and exhibit uniform maximum transcendentality [42,43], which has for a long time been conjectured to characterize scattering amplitudes in N = 4 SYM [44].While these properties do not translate in an obvious way to non-supersymmetric theories like QCD, it has been shown that such integrals can still substantially simplify the calculations of scattering amplitudes within the Standard Model.Namely, integrals of this type fulfil particularly simple systems of differential equations, in the so-called canonical form [45]. Canonical sets of equations are more easily solved and provide a direct handle on the analytic properties of the corresponding integrals in various singular regions.Importantly, the solution of a canonical system with rational coefficients is straightforwardly expressed in terms of a well-understood class of functions, the Multiple Polylogarithms (MPLs) [46][47][48].
Among the initial applications of this formalism were planar ladder-type integrals in H+jet production up to three loops [49], which include the planar two-loop integrals up to O(ϵ 2 ) as a subset.Here we consider the full set of two-loop integrals: planar and nonplanar.We construct a pure basis of uniform transcendental weight and demonstrate that it can be solved in terms of MPLs to any order in ϵ.At variance with [39,40], we show that using a generalized set of regularity conditions on the canonical master integrals, all boundary conditions required to fix the solution of the differential equations can be inferred in terms of a small number of one-scale two-and three-point functions which are known in closed form in the literature.In view of applications at N 3 LO, we limit ourselves to perform the calculation explicitly to order ϵ 2 , which corresponds to transcendental weight six.
The structure of the paper is as follows.The effective coupling of the Higgs boson to light partons and the definition of kinematics is given in Section 2. In Section 3, we formulate the general structure of the amplitude and use projectors to obtain tensor coefficients and construct the helicity amplitudes.In Section 4, we describe in detail the construction of pure bases for the two-loop integral families and solve their canonical differential equations analytically.The ultraviolet renormalization and the subtraction of infrared singularities of our amplitudes are discussed in Section 5. Finally, the crossing of the helicity amplitudes to all appropriate kinematic configurations and the necessary analytic continuation of the relevant multiple polylogarithms is discussed in Section 6.We conclude with a brief summary in Section 7.

The effective Lagrangian
The Higgs boson interacts with Standard Model particles with a coupling strength proportional to their mass, and therefore cannot couple directly to gluons or massless quarks.Nevertheless, starting at one loop, the Higgs can interact with gluons through virtual loops of massive quarks.In the limit of a very heavy quark mass, m q → ∞, one can show that this coupling becomes independent of m q and an effective theory can be formulated by integrating out the corresponding quark from the full theory [8][9][10].While most quarks have relatively small masses compared to the typical energy scales of scattering processes at the LHC and can often be considered as massless, the same is not true for the top quark.Since it is the heaviest of all known Standard Model particles, this effective field theory (EFT) works extremely well for the top quark, at least as long as all scales involved are smaller than twice its mass [50,51].
In the following, we will work in the EFT with a matter content of N f = 5 massless quarks and one very heavy quark, the top quark, integrated out.In this case, the effective Lagrangian becomes where G µν a is the field strength tensor of the gluons and H is the Higgs field.From dimensional analysis, the effective coupling λ has inverse mass dimension.It was shown long ago [52,53] how to perform the matching of this effective theory to the Standard Model Lagrangian [6,54,55].

Kinematics
We are ultimately interested in computing the amplitude for the production of a Higgs boson and a hadronic jet in parton-parton annihilation at the LHC.For simplicity, we start considering the problem in the crossed kinematics which corresponds to the decay of a Higgs boson into three partons.There are two relevant partonic channels, namely the decay into three gluons and into a quark, anti-quark and a gluon H(p 4 ) → q(p 1 ) + q(p 2 ) + g(p 3 ). (2. 3) The amplitudes for the production processes can then be obtained through an analytic continuation of the decay kinematics [56], see Section 6 for more details.
The Mandelstam invariants are defined as and satisfy the conservation equation where M H is the mass of the Higgs particle.It is more convenient to work with dimensionless ratios such that (2.5) implies the relation, In the decay kinematic region, all these invariants are non-negative.This, together with (2.6), defines the corresponding kinematic region (2.8)

Tensor Decomposition
Following earlier work on this process [29], we decompose the amplitudes M ggg and M qqg as where we used ϵ i to denote the polarization vectors of external gluons.The above tensors can be expanded perturbatively in the QCD coupling constant α s as T µ (q, q, g) = λ √ 4πα s T a ij T (0) µ (q, q, g) + α s 2π T (1)  µ (q, q, g) where the coefficients S Given the external states and their possible helicity and spin configurations, the amplitudes S µνρ and T µ can only depend on a limited number of tensor structures.These structures can be further constrained by exploiting symmetries and choosing a gauge or, following [29], by enforcing gauge invariance through the Ward identities.In this Section, we aim to find such a tensor basis in order to be able to work with the scalar coefficients of the amplitudes with respect to this basis, known as the form factors.The form factors, in turn, are obtained by applying projector operators on the full amplitude expanded in Feynman diagrams.We derive a basis of 4 tensor structures for H → ggg and a basis of 2 tensor structures for H → q qg, and the corresponding projectors in Subsections 3.1 and 3.2, respectively.
Since we are ultimately interested in fixing the helicities of external states and computing the associated helicity amplitudes, we introduce the spinor-helicity formalism in Subsection 3.3 and derive a set of genuinely independent helicity amplitudes.These are written as a unique spinor factor times a scalar coefficient which is a linear combination of the form factors.The same linear combination of form factor projectors also defines a helicity amplitude projector.

Tensor decomposition for H → ggg
Let us start considering the decay of a scalar boson into three massless spin-one particles.The most general tensor structure one can build using the four-vectors associated with the external particles is There are four helicity configurations for this amplitude, so we expect four independent tensor structures in D = 4. Indeed, the transversality conditions p i • ϵ i = 0, i = 1, 2, 3 and the cyclic gauge choice ϵ 1 • p 2 = 0, ϵ 2 • p 3 = 0, ϵ 3 • p 1 = 0 restrict the tensor structures considerably: where we relabelled the coefficients {A i , ..., D i } as the form factors G i and defined the basis The form factors can be obtained from a Feynman diagram decomposition of the amplitude S (i) µνρ at any loop order by applying suitable projectors P i defined as where polarization vectors satisfy the gauge-fixed polarization sum, with reference vectors defined above.The projectors can in turn be decomposed in terms of the dual of the tensor basis in (3.5): To work out the projectors explicitly, we insert the decompositions (3.8) and (3.5) into the definition (3.7) and obtain the requirement which is satisfied by the coefficients c In particular for the amplitude H → ggg with external states in D dimensions [32], we get (3.10)

Tensor decomposition for H → q qg
Similarly, it is easy to see that the most general tensor decomposition for two external spinors and a four-vector is where we sum over odd products of gamma matrices.In the first sum, all indices but one are contracted amongst each other or with external momenta while in the second, there are no indices left, and the polarization vector contracts with an external momentum vector.
In general, one might expect the length of the spinor chains to be bound by the loop order.
By simple enumeration, one can see that at tree level there can only be one gamma matrix, at one loop up to three gamma matrices and at two loops up to five.Nonetheless, it is easy to see that with the momenta at hand, one cannot build any spinor chain with more than one Dirac γ matrix.Enforcing the transversality condition p 3 • ϵ 3 = 0 and gauge choice ϵ 3 • p 1 = 0, one is left with the decomposition Hence we define the two tensor structures The form factors are extracted from the amplitude using the projectors P i with where polarization vector of the gluon satisfies the gauge-fixed polarization sum, with reference vector defined above.Following the strategy outlined in Subsection 3.1, we get for the two projectors

Helicity amplitudes
From the tensors found above, we can easily obtain compact expressions for the relevant helicity amplitudes.We work in the 't Hooft-Veltman scheme and consider external states as four-dimensional and assume fixed helicity states.In massless QCD, both gluons and quarks have two helicity configurations, λ i = ±, and the amplitudes can be written as The two helicity states of a four-component massless spinor are projected out through and we fix the spinor-helicity bracket representation of incoming fermions as and of incoming anti-fermions as For massless vector bosons, incoming states of positive and negative helicity are given by where the reference momentum r is an arbitrary light-like vector such that r • p ̸ = 0.
Outgoing states have the same representation with switched helicities.Let us start by considering the decay H → ggg.There are two independent helicity configurations, which we choose to be M +++ ggg and M ++− ggg , while the other helicity amplitudes are obtained through parity conjugation and by relabelling of the gluon momenta.Starting from the decomposition (3.5) with the basis (3.6) and applying the definitions (3.18)-(3.20),we can cast the independent helicity amplitudes in terms of spinor products, where the coefficients α and β are simple linear combinations of the original form factors: The helicity projectors for α and β are built by replacing the G i with the corresponding form factor projectors from (3.10): Let us now consider the other partonic process H → q qg.In this case, there is only one independent helicity amplitude, which we choose to be M LR+ q qg .In this helicity configuration, the first tensor in (3.12) vanishes, while the second tensor yields where γ = s 12 F 2 . (3.28) Again, the helicity projector is easily obtained by replacing F 2 with the associated projector from (3.15), Just like full amplitudes in (3.2) and (3.3), the helicity amplitude coefficients α, β, γ also have a perturbative expansion, for Ω = α, β, γ.The overall colour factors for the two processes are T α = T β = f a 1 a 2 a 3 and T γ = T a 3 i,j .The tree-level helicity amplitudes are well-known and the coefficients evaluate to

Master integrals
The form factors F i and G i , which relate to the helicity amplitudes via (3.24), (3.28), receive contributions from all relevant Feynman diagrams at a given perturbative order.All treelevel, one-and two-loop diagrams with H and ggg or q qg external states are generated with standard QCD vertices and the effective Higgs interactions using QGRAF [57].We shift each diagram to a kinematic crossing of one of the auxiliary topologies presented in Table 1 using Reduze2 [58,59].After inserting Feynman rules and evaluating the Dirac and Lorentz algebra in FORM [60], the contribution of each diagram to the form factors in (3.5) and (3.12) can be written as a combination of scalar integrals of the form with k l the loop momenta, D i ∈ {P i } ∪ {N i } the internal propagators, and γ E = 0.577 . . . the Euler-Mascheroni constant.The topology of each scalar integral is uniquely identified by the propagators which appear in the denominator (a i ≥ 1).This information can be used to compute the diagram's sector ID within one of the auxiliary topologies in Table 1 as a binary nine-bit number.The scalar integrals defined in (4.1) satisfy integration-by-parts (IBP) identities [33,34], allowing them to be expressed in terms of a minimal set of so-called master integrals.Note that the physical measure in the bare amplitude is (e) NPL 487 loop compared to the integration measure in (4.1), which then requires a simple conversion factor when inserting the solutions of the master integrals into the amplitude.
To evaluate the master integrals, we use the method of differential equations [35][36][37][38], augmented with a canonical basis [45] for the master integrals.In particular, we reduce all scalar integrals appearing in the amplitudes using Reduze2 [59] and Kira [61,62] directly to a canonical basis which involves 89 master integrals.Our basis is different from the one considered in [39,40], but relations amongst the two sets of integrals can easily be established.There are 4 planar and 2 non-planar sectors with the maximum number of propagators, t = 7, depicted in Figure 1.85 of the 89 IBP master integrals are subsectors of the non-trivial top sectors (d)-(f).They are kinematic crossings of the 16 planar and 8 non-planar topologies depicted in Figures 2 and 3, which will be explained in more detail after introducing the canonical basis.

Canonical basis
Multiple public packages exist for the derivation of a candidate canonical basis like CANON-ICA [63], Fuchsia [64], or DlogBasis [65].We use DLogBasis to provide a list of candidate UV-finite integrals with unit leading singularities for the planar family.For the non-planar family, we constructed canonical candidates by studying the leading singularities of the relevant master integrals, following a loop-by-loop approach.Note that for all planar and some non-planar integrals, canonical candidates can be obtained just by rescaling a single integral in the sector by its maximal cut, i.e. a solution to the homogeneous part of its differential equation [66].This can be confirmed by studying the maximal cut of the specific canonical integral in the Baikov representation [67].The canonical bases sufficient for the computation of all integrals in the families PL and NPL are presented in the supplementary material.The amplitude, however, contains diagrams which produce integrals in kinematic crossings of the families in Table 1.To obtain a basis that is sufficient to represent the physical amplitude, we proceed from the lowest sectors in applying crossings to the canonical integrals in PL and NPL, and appending to our basis those canonical integrals whose reduction contains new masters.This yields a minimal set of crossed and uncrossed canonical master integrals, sufficient for physical applications, whose generic topologies are depicted in Figures 2 and 3. Note that we obtain results first in the Euclidean region, where all invariants are negative s ij < 0. For simplicity, we also set M 2 H = −1 in explicit formulas below.
Up to now, we only considered the non-trivial top sectors (d)-(f) in Figure 1, but the amplitude requires the reduction of integrals from all six topologies.While integrals in (a) and (b) can be reduced to integrals considered earlier, (c) contains four new master integrals, I 86 , . . ., I 89 which need to be added to the basis.They are the integrals (the numbering refers to the canonical basis in the supplementary material): and The full expressions up to O(ϵ 6 ) for I 86 as well as I 87 and its two crossings are given in the supplementary material.In this way, one can complete the set of 89 master integrals sufficient to reduce any integral in this process.The full canonical basis is given in the supplementary material accompanying this paper.

Solution of differential equations
For the purpose of computing the master integrals, we consider first the 16 integrals in the tree of top sector (d) in PL and the 36 integrals in the tree of the non-planar top sectors  (e), (f) in NPL.Separately for the two auxiliary topologies, we compute the derivatives of the candidate canonical combinations, and insert the IBP reduction to obtain differential equations in the following form where A i , B i are sparse matrices of rational numbers.It is obvious from this form that the solutions for the canonical combinations can be expressed in terms of MPLs [46] with the alphabet {y, z, y − 1, z − 1, y + z, 1 − y − z}, which are usually written out in terms of a fibration in either y or z.We recall here that MPLs are defined as iterated integrals over rational functions with G(x) = 1.In this context, n is referred to as the transcendental weight of the polylogarithm.By construction, the method of differential equations cannot be used to compute purely one-scale integrals directly.These must instead be obtained from alternative methods and added to the system of differential equations.At two loops, there are five such two-and three-point functions: and their kinematic crossings.Analytical expressions for all these one-scale integrals are known in closed form in the dimensional regulator ϵ [38,68].
We approach the solution of the differential equations for the remaining integrals as follows.At each order n in ϵ, we consider the vector of master integrals ⃗ I (n) (y, z) and by choice integrate first the equations in the variable y.If the set of differential equations in the two variables are consistent, we obtain a partial solution which differs from the full solution ⃗ I (n) (y, z) only by a function of the other variable, ⃗ f (z), The intermediate result is then substituted into the set of equations in z and solved for ⃗ f (z), which fixes the final solution up to a numerical constant As stated above, (4.15) cannot depend on y and the spurious dependence must cancel from the right-hand side of the equation.
Figure 3: The non-planar topologies in the canonical basis.

Fixing boundary conditions
The remaining numerical constants ⃗ c can in principle be fixed by evaluating the integrals at special kinematical points.We opted to obtain all boundary conditions without any additional computations simply by imposing a set of regularity conditions on the general solution (4.14).As suggested in [39,40], a large set of boundary conditions can be obtained exploiting regularity of the master integrals at various pseudo-thresholds.This can be achieved in practice by multiplying the differential equations (4.4) and (4.5) by those letters which correspond to pseudo-thresholds of the integrals, and taking the limit as follows: left-hand side : lim right-hand side : lim Requiring the right-hand side to vanish yields non-trivial relations between the integrals in the limit x → l i .Namely, if the rational factor in question appears in the homogeneous term of the differential equation for a given master integral, its value at a regular kinematic point can be typically related to other integrals in the same sector and its subtopologies.This approach is sufficient for the planar topology.
In order to impose these regularity conditions, we used PolyLogTools [69] to manipulate multiple polylogarithms up to weight 5, evaluate the required limits, and perform changes in the fibration basis.Beyond weight 5, we had to carry out the required fibrations ourselves, building on the implementation of differentiation and integration of MPLs in PolyLogTools.In particular, we differentiated the integrals with respect to the variable we intend to fibrate into and obtained linear combinations of MPLs of weight 5, which can be treated with automated routines.The result can subsequently be integrated back and expressed in the required form, up to an integration constant.All constants can be fixed by comparing the original and fibrated function at a kinematic point, and reconstructing their difference as a number of the appropriate weight using the PSLQ algorithm [70].In practice, we find that our definition of the variables y and z, (2.6), allows us to consider just the 3 limits y → 1, y → 0, and y → −z.
In the NPL topology, integrals possess branch points when y → 0, z → 0 and 1−y−z → 0. For this reason, the strategy outlined above can only be applied to a small number of letters and one can show that these conditions are not enough to fix all remaining constants.Taking inspiration from [65,71], we also consider the singular limits (i.e.genuine thresholds of the master integrals).In terms of the Mandelstam variables, these correspond to the additional limits s → 0, t → 0 and u → 0. Crucially, our canonical basis consists of UVfinite integrals.Hence, we only need to regulate IR divergences and can assume that ϵ < 0. Within this condition and keeping ϵ fixed, if a linear combination of integrals develops a singular behaviour for one of the singular limits, this must correspond to a spurious UV-type divergence.Since no new UV divergences can appear when the kinematical invariants take special values, we may impose that such spurious divergences do not occur at the kinematic points where one of the letters vanishes.In this way, we obtain yet more equations between the boundary constants.Explicitly, the solution to the DE near the point y → l i is

.18)
The matrix exponential contains elements with terms of the type (y − l i ) aϵ .If a > 0, such an expression diverges in the limit y → l i .Since all our integrals must be finite at this kinematic point, the constants ⃗ I y=l i ought to take specific values such that these terms cancel.
We also studied the asymptotic behaviour of our master integrals with asy [72], which is implemented in FIESTA [73], and verified the ϵ-dependence in the limit y → l i .Note that in [40], it was required to integrate one order in ϵ higher than necessary to enforce the nonappearance of spurious singularities at the previous order, which is no longer required.We checked our solutions for all top sector master integrals numerically against pySecDec [74] for several Euclidean points up to weight six and found perfect agreement.
In Appendix B, we present the computation of the one-loop master integrals to order O(ϵ 4 ), which also serves as a simple example of some of the techniques discussed in this section.

UV renormalisation and IR regularisation
The bare helicity amplitudes (3.30) contain ultraviolet (UV) as well as infrared (IR) divergences that manifest as poles in the Laurent expansion in the dimensional regulator ϵ.The former are treated in the MS scheme by expressing the amplitudes in terms of the renormalized couplings, α s ≡ α s (µ 2 ) and λ ≡ λ(µ 2 ), evaluated at the renormalization scale µ 2 .The resulting amplitudes still contain IR singularities, which will be cancelled analytically by those occurring in radiative processes of the same order [75,76].Their structure is universal and it was originally determined up to two loops by Catani [77,78].These results were later systematised and extended to general processes and up to three loops in [79][80][81][82][83][84][85][86][87].
In this Section, we present the necessary steps and formulae to perform the UV renormalization and subtraction of the IR poles.This allows us to obtain the one-loop and two-loop finite remainders, which we decompose according to their colour structure.In contrast to [29], where the IR subtraction was performed in the Catani scheme, we followed a subtraction scheme based on Soft-Collinear Effective Theory [82,83], which can be more naturally extended to higher loops.In subsection 5.5 we provide conversion formulas between the two different schemes.

Ultraviolet renormalization
We start by denoting all unrenormalized quantities with a superscript U and then replace the bare coupling α U with the renormalized strong coupling α s ≡ α s (µ 2 ), evaluated at the renormalization scale µ 2 , where S ϵ = (4π) ϵ e −ϵγ E and µ 2 0 is the mass parameter in dimensional regularization introduced to maintain a dimensionless coupling in the bare QCD Lagrangian density.The explicit form of the first two β-function coefficients β 0 , β 1 reads with the QCD colour factors, The effective coupling λ is renormalized as follows, The renormalized coefficients of the UV-finite but IR-divergent amplitudes can be written in terms of the i-loop contribution to the unrenormalized coefficients Ω (i), U as ) ) For the remainder of the paper, we will set µ 2 = µ 2 0 = M 2 H for simplicity.

Infrared factorization
Since the IR poles of l-loop amplitudes in gauge theories factorize in colour-space in terms of lower loop amplitudes, the IR poles can be subtracted multiplicatively as Contrary to multiplicative renormalization of UV divergences, the bold notation in (5.9) indicates that, in general, Z and Ω are operators and vectors in colour space, respectively.The all-order nature of (5.9) makes this approach particularly advantageous for generalizations to higher orders in perturbation theory.Solving a renormalization group equation for Z, one finds where P is the path-ordering symbol, meaning that the colour operators are ordered from left to right in decreasing values of µ ′ .As originally proposed by Catani [78], the anomalousdimension matrix Γ for amplitudes with n QCD partons up to two loops is entirely governed by the dipole colour correlations operator where γ cusp is the cusp anomalous dimension and γ i is the anomalous dimension of the i-th external particle, the latter depending on the nature of the particle.The cusp anomalous dimension carries the information about overlapping soft and collinear divergences, while γ i only involves collinear divergences associated to the i-th parton.The perturbative expansions up to two loops for the anomalous dimensions are listed in Appendix A. The coupling constant is evaluated at the renormalization scale α s = α s (µ).The colour operators T a i are related to the SU (N ) generators and their action on the i-th coloured parton is defined following the convention in [78] as: where i labels the particle on which the operator is acting.Importantly, colour conservation can be rephrased as i It follows from these definitions that the repeated action of one operator evaluates to a Casimir, where We also define the expansions, where one can drop the bold notation in the derivative because the resulting operator is always diagonal in colour space: (5.17) In our specific case, expanding (5.9) up to two loops, IR divergences in the renormalized two loop amplitudes can be expressed in terms of the renormalized tree and one-loop amplitudes multiplied by appropriate operators ) Ω Ω (0) , (5. 19) Ω Ω (0) − I Ω Ω (1) , ( With the definitions given in (5.10)-(5.16),we can express the subtraction operators in (5.19)-(5.20)as where the process dependence is implicit and It is now manifest that the IR operators in SCET scheme have only pole terms, without any O(ϵ 0 ) contributions.The cancellation of IR poles according to (5.20) is a strong analytic check on multiloop amplitudes.We shall now follow the SCET approach to derive explicit expressions for the subtraction operators.
Accordingly, we can drop the bold notation and we find expressions relevant to the three helicity coefficients α, β, γ: with and the anomalous dimension coefficients defined as in Appendix A.
We stress here that there is an imaginary part arising from the logarithms above whenever the corresponding invariant is positive, which depends on the kinematical region considered.

Results for the helicity amplitudes
We computed the renormalized amplitudes for the decay process up to order ϵ 4 at one loop and up to order ϵ 2 at two loops.Working in the SCET subtraction scheme, we derived the finite remainder for all the helicity amplitudes.Decomposed according to their colour structure, they read Ω + N 2 F F (2) Ω . (5.30) The same structure holds for the renormalized amplitudes, though the respective coefficients will still contain poles in ϵ.
In the supplementary material, we provide the coefficients of the renormalized amplitudes for the decay processes and of the finite remainder for the Higgs decay kinematics as well as for all crossings [56] relevant to H+jet production processes.

Conversion to the Catani scheme
The two-loop helicity amplitudes for the decay of a Higgs boson into three partons were first computed in [29].The authors obtained the finite remainder by subtracting the IR singularities according to the original Catani prescription [78].Here we give the conversion rules between the two subtraction schemes, which also served as a cross-check of our results.
Following closely the notation of [29], we write the subtraction operators I Ω of (5.19) in Catani scheme as with, (5.32) The second-order operator can be built starting from the one-loop operator as 1) Ω (ϵ) I C,( 1) 1) 1) Ω (ϵ) , (5.33) where we introduced the constant (5.34) The remaining term in (5.33) involves the operator H Ω (ϵ) and produces only a single pole in ϵ.Its explicit form is Ω . (5.35) The constant H Ω is renormalization scheme and process dependent and in our case it reads where in the MS scheme the constants H q , H g are (5.39) When describing the SCET subtraction scheme, we pointed out that the subtraction operators have no finite O(ϵ 0 ) contribution.
In contrast, the Catani operators contain coefficients at higher order in the dimensional regulator, generated by the ϵ expansion of the resummed coefficient defined in Eq. (5.32).The tree-level amplitude is finite.One-loop amplitudes have poles starting from order O(1/ϵ 2 ) and two-loop amplitudes have poles starting from order O(1/ϵ 4 ).We indicate with Ω (1) (1) (1) 0 . (5.42) In deriving the above rules, we made use of the fact that I C, (1) and I SCET, (1) have the same pole structure.Consequently, terms multiplying (I C,(1) ) or (I C,(1) ) in (5.42) vanish.This can be easily understood by inspecting the origin of the poles in the one-loop cancellation in (5.19).
For both the decay and the production kinematics, we verified that the result given in [29] is correctly reproduced converting our finite remainder with the above rules.

Analytic Continuation
So far, we have described our calculation making explicit reference to decay processes, for which all the kinematic invariants are positive.We considered the decay of a Higgs boson into three gluons, H → ggg, and into a quark-antiquark pair and gluon, H → q qg.In view of applications to LHC physics, we are interested in the production regions, in which a Higgs is produced together with a parton: gg → Hg, qg → Hq, qg → H q and q q → Hg.
The general strategy for performing the analytic continuation for the MPLs appearing in 2 → 2 scattering involving 4-point functions with one external off-shell leg and massless propagators was outlined in detail in [56].Our aim is to describe how this strategy can be applied to the process at hand.Referring to Fig. 1 of reference [56] and using the same labels for the various kinematic regions, the goal is to find a procedure to analytically continue MPLs from the decay region (1a) to the three production regions, (2a), (3a), (4a).Whenever a particle is crossed from the initial state to the final state (or vice versa), two of the invariants become negative and one remains positive, representing the centre-of-mass energy of the incoming partons.We recall here that results in the decay region (1a) are expressed as MPLs of the variables y and z defined in (2.6), which fulfil the constraints 0 < z < 1 , 0 < y < 1 − z.Within these bounds, our set of MPLs are real and no branching point is crossed.
To describe the analytic continuation to the scattering kinematics, let us consider the case of region (3a).In (3a), the kinematic constraints are z < 0 , 1 − z < y < +∞ and MPLs must be evaluated across a branch cut, developing an imaginary part.Physically, the particle with momentum p 2 is crossed and the process is p 1 + p 3 → p 2 + p 4 .Crucially, there exists a change of variables which maps region (3a) linearly back into the decay region (1a) and can be implemented analytically on our MPLs.In fact, by defining the new variables u and v satisfy again 0 < u < 1 , 0 < v < 1 − u.By re-expressing the MPLs in terms of u and v, the imaginary part can be made explicit in terms of multiple zeta values and real-valued MPLs.The same manipulations can be performed in regions (2a) and (4a) with different definitions of the (u, v) variables.In general, the v variable is the reciprocal of the centre-of-mass energy and so its definition depends on which particle is crossed from the final to the initial state.Instead of performing the analytic continuation in all the three regions (2a), (3a) and (4a), we found it simpler to first consider suitable crossings of the amplitude in the decay Process Parity Related In the production region, there are eight different helicity configurations.Thanks to parity symmetry, we can limit ourselves to consider just four of them (see the left column of Table 2) and relate them to the other four (right column).Note that in continuing to the region (3a), the momenta p 1 and p 3 are always in the initial state.We computed the first three amplitudes with a combination of crossings and analytic continuation as follows (an extra helicity flip due to time reversal is always understood after continuation to (3a)): The fourth amplitude, M +−,+ gg,g , can be derived from M −+,+ gg,g by the crossing p 1 ↔ p 3 , which implies v → v and u → 1 − u − v. Since no branch cut is crossed under this transformation, no new analytic continuation is needed and we can conclude that M +−,+ gg,g (v, u) = M −+,+ gg,g (v, u)| u→1−u−v .(6.7) Let us consider now the decay into a quark-antiquark pair and a gluon: M LR+ : H → q L (p 1 ) + qR (p 2 ) + g + (p 3 ) .(6.8) In the production region there are 12 non-zero independent helicity configurations, but only 3 of them need to be computed, while the others are related by parity and charge conjugation.Our choice of the independent configurations is according to the left column of Table 3.
verified that our results for the finite remainder of the amplitude match the literature [29], and analytically continued the helicity amplitudes to all kinematic regions.The infrared structure of the result was inspected in the frameworks of SCET and the Catani infrared factorization formula.This updated result marks the first step towards computing the N 3 LO corrections to H+jet production.

B One-loop master integrals
In order to obtain a result for the two-loop amplitudes up to O(ϵ 2 ), we ought to calculate the one-loop master integrals and amplitudes up to O(ϵ 4 ).Denoting Clearly requiring the non-appearance of terms y ϵ in the solution replicates the condition (B.9).The strength of the second approach is that it applies to regular as well as singular limits.Expanded solutions for all four integrals can be found in the supplementary material.
the i-loop contributions to the amplitude.The SU (3) group generators are normalised as Tr(T a T b ) = δ ab /2.

Figure 1 :
Figure 1: The 4 planar and 2 non-planar sectors with t = 7, labelled with the propagators P i and N i , respectively, as listed in Table 1.Their names indicate the integral family and sector ID.The sectors (d)-(f) contain new master integrals and are the focus of this Section.Sectors (a) and (b) are entirely reducible to master integrals found in the trees of the former three top sectors.Sector (c) can be reduced to masters from (d)-(f) and four additional integrals which are easily computed from one-loop results.

Figure 2 :
Figure 2: Planar topologies which appear in the canonical basis.The dashed line is the massive leg, dots represent squared propagators.A propagator in brackets denotes an integral with this propagator in the numerator.
order ϵ n of the renormalized amplitude, the SCET operator and the Catani operator, respectively.From the subtraction formulae (5.19)-(5.20), it is easy to obtain the following conversion rules for the finite remainders Ω

Table 1 :
The two auxiliary topologies of momenta which label the propagators of every diagram appearing in the amplitudes.PL labels exclusively planar diagrams, whereas NPL accommodates also all non-planar sectors.