Two-loop leading-colour QCD helicity amplitudes for Higgs boson production in association with a bottom-quark pair at the LHC

We compute the two-loop QCD helicity amplitudes for the production of a Higgs boson in association with a bottom-quark pair at a hadron collider. We take the approximations of leading colour and work in the five flavour scheme, where the bottom quarks are massless while the Yukawa coupling is non-zero. We extract analytic expressions from multiple numerical evaluations over finite fields and present the results in terms of an independent set of special functions that can be reliably evaluated over the full phase space.


Introduction
Precise theoretical predictions are an essential ingredient in the search for subtle deviations from the Standard Model (SM) at current and future collider experiments. The enormous amount of data gathered at the Large Hadron Collider (LHC) is already challenging the theoretical precision for many processes and the pressure will increase as data continues to pour in. It has been clear for some time [1][2][3] that, for a large class of two-and threeparticle final states, at least the Next-to-Next-to-Leading Order corrections in Quantum Chromodynamics (NNLO QCD) will be required for fully differential cross sections. The high multiplicity final states, particularly those with many kinematic scales, pose the biggest technical challenge owing to the enormous analytic and algebraic complexity. Nevertheless, thanks to a huge effort across the theoretical community, new tools and methods have been developed that have produced the first NNLO QCD predictions for 2 → 3 processes [4][5][6], most notably for 3-jet production [7]. This remarkable progress has been driven both by the advancements in the scattering amplitude computation and by efficient subtraction schemes. By now, all the two-loop master integrals -which are one of the important ingredients in the amplitude computation -are available for the fully massless five-particle processes [8][9][10][11][12][13], allowing for several two-loop QCD amplitudes to be derived analytically [14][15][16][17][18][19][20][21][22][23][24][25], improving on previous results that were obtained numerically [26][27][28][29]. These new results have been achieved thanks to technological breakthroughs in the method of differential equations [30][31][32][33][34], integral reduction algorithms [35][36][37][38][39] and the use of finite-field arithmetic to tame the algebraic complexity of multi-leg and multi-scale problems [40][41][42][43][44][45].
A different class of 2 → 3 processes that is of great interest involves an external massive leg. The two-loop planar helicity amplitudes for W + 4-parton scattering (contributing to the prediction for W + 2-jet production at NNLO QCD) have been previously studied numerically [46]. All two-loop integrals needed for this type of process are available for the planar contributions [8,[47][48][49], and the first analytic result was derived for W bb production in the leading colour, massless b-quark and on-shell W -boson approximations [50]. Very recently, complete analytic results for one of the non-planar integral topologies have also become available [51,52].
In this article we consider the two-loop amplitudes relevant for the production of a Higgs boson in association with a bottom-quark pair, i.e. pp → bbH, in the leading colour approximation. bbH production at the LHC has been a subject of great phenomenological interest due to its potential in directly measuring the bottom-quark Yukawa coupling. In the Standard Model (SM), the coupling strengths of the Higgs boson to the fermions and vector bosons are proportional to their mass, causing the rate of the bbH production to be suppressed with respect to, for example, Higgs production in gluon fusion (gg → H) or vector boson fusion (pp → Hjj), associated production with a vector boson (pp → V H), and associated production with a top-quark pair (pp → ttH). In addition, the presence of btagging further suppresses the bbH production rate. In some new physics scenarios, such as the Two Higgs Doublet Models (2HDM's) and the Minimal Supersymmetric Standard Model (MSSM), the bottom-quark Yukawa coupling can be dramatically enhanced, resulting in a considerable increase of the bbH production rate [53,54]. Thus, the study of the bbH production will allow to constrain supersymmetric models and other extensions of the SM that modify the bottom-quark Yukawa coupling. Recent studies on the interplay between bbH signal and backgrounds can be found in Refs. [55][56][57].
The theoretical approach to obtaining predictions for the pp → bbH process has been subject of much discussion in the community. This is due to the fact that, in the presence of bottom quarks, a theoretical prediction can be computed in either the four-flavour scheme (4FS) or the five-flavour scheme (5FS). In the 4FS computation, bottom quarks are treated as massive and they do not contribute to the parton distributions functions (PDFs), hence only appearing in the final state. Large logarithms of the form log(m b /Q) with Q ∝ m H arise when the integration over the bottom-quark phase space is performed, and such contributions may spoil the convergence of the perturbative series. These large logarithms can be resummed to all orders by introducing the bottom-quark parton distribution functions. The 5FS approach stems from this prescription, where the bottom-quarks are included in the parton distribution functions, allowed to appear in the initial state, and treated as massless. We refer the reader to Ref. [58] for further discussion on the 5FS and 4FS approaches. In 5FS the inclusive bbH production (where the tree level process is bb → H) has been computed up to N3LO QCD [59][60][61][62][63][64][65][66][67], while for the case where a single bottom quark is observed NLO QCD [68], weak [69] and SUSY QCD [70] corrections are available. In 4FS the bbH production has been calculated up to NLO QCD [71][72][73][74][75][76], and the supersymmetric QCD corrections [77,78] are also known. There have also been efforts in matching the 5FS and 4FS calculations to obtain accurate predictions across the entire kinematic region [79][80][81][82][83]. A first step towards a massive version of the five-flavour scheme (5FMS) has been devised to naturally connect the 4FS and 5FS approaches [84,85].
In this work we compute the two-loop QCD corrections to the gg → bbH, qq/qq → bbH, bb/bb → bbH, bb → bbH andbb →bbH reactions in the 5FS. These two-loop amplitudes enter the computation of pp(bb) → H at N4LO, pp → b(b)H at N3LO when one b-jet is tagged, and pp → bbH at NNLO when two b-jets are required in the final state. We note that beyond NLO, for the computation with massless bottom quarks, a flavoured jet algorithm [86] would have to be employed when identifying the b-jets, since the use of conventional k T or anti-k T jet algorithms would render the fixed order computation infrared unsafe. We further remark that the two-loop amplitudes for bbH production derived in this work can also be used in the computation of Higgs decaying into a bottom-quark pair in 5FS, by crossing initial partons to the final state. Specifically, they will contribute in the N4LO H → bb, N3LO H → bbj and NNLO H → bbjj computations. In addition, by crossing the bb pair to the initial state and the gg/qq pair to the final state we obtain the contribution of the bottom-quark initiated channel to H + 2j production (bb → Hjj).
We present analytic results for the finite remainders after ultraviolet (UV) and infrared (IR) poles have been subtracted. This is possible using a basis of independent special functions recently identified in the context of W bb production [50]. We obtain numerical results valid across the full phase space by applying the generalised series expansion approach [47,87,88] to the differential equations satisfied by the special functions appearing in the finite remainders.
This paper is organised as follows. We begin by describing the structure of the bbH amplitudes at leading colour in Section 2, followed by a discussion of the methodology used in deriving the analytic expression of the amplitudes in Section 3. A description of the basis of special functions is presented in Section 4 and a number of validations that we performed on the results derived in this work are discussed in Section 5. We present benchmark numerical evaluations together with evaluations on a physical phase-space slice in Section 6. Finally we draw our conclusions in Section 7.

Structure of the pp → bbH Amplitudes at Leading Colour
We compute the two-loop QCD corrections in the leading colour approximation for the following subprocesses where all momenta are taken as outgoing, We work in the five-flavour scheme, where the bottom quark is taken as massless while its Yukawa coupling to the Higgs boson is kept finite, so that where m H is the mass of the Higgs boson. The kinematics is described by six independent scalar products, which we choose as The latter is connected to the Gram determinant of the four linearly independent momenta, The colour decomposition of the L-loop amplitudes in the leading colour approximation is given by where n = m α s /(4π), α s = g 2 s /(4π), m = i(4π) e − γ E , T a are the fundamental generators of SU (N c ) normalised such that tr(T a T b ) = δ ab , while g s and y b are the strong coupling constant and bottom-quark Yukawa coupling, respectively. We further decompose the partial amplitudes at one and two loops based on their closed fermion loop contributions, as  The pole structure of the unrenormalised amplitudes in the 't Hooft-Veltman (tHV) scheme at one and two loops is given by

8)
whereÂ (1) is the unrenormalised one-loop amplitude normalised to the tree-level amplitude. s 1 and s 2 are the bottom-quark Yukawa renormalisation constants, and their expressions can be found in Appendix A. We used a mixed renormalisation scheme where the strong coupling α s and the bottom-Yukawa coupling y b are renormalised in the MS scheme, while the bottom-quark mass and wave function are renormalised in the on-shell (OS) scheme. This allows to keep y b finite while taking the bottom-quark mass smoothly to zero (m OS b → 0) [67]. Such a mixed renormalisation scheme can be used so long as pure QCD corrections are considered. In fact, using the MS scheme to renormalise y b allows us to better control the convergence of the perturbative corrections by resumming the large logarithms that appear in the OS scheme by running y b to a scale close to the Higgs mass. In the presence of electroweak (EW) corrections, however, the relationship between y b and m b must be imposed to guarantee the cancellation of UV singularities [55].
The I 2 ( ) operator is defined by while the I 1 ( ) operators for bbH production in both the gg and the qq channels are given at leading colour by 14) where The β function coefficients and anomalous dimensions are given in Appendix A. The finite remainder of the L-loop partial amplitude is then obtained by subtracting the poles P (L) (which include both the ultraviolet and infrared singularities) from the unrenormalised partial amplitude A (L) and setting to 0, The partial finite remainders F (L) inherit from the partial amplitudes the decomposition in powers of n f , The full finite remainders F (L) are obtained from the partial ones F (L) through a colour decomposition analogous to that given in Eq. (2.7) for the bare amplitudes,

Tree-Level Amplitudes
The tree-level amplitudes can be obtained using the BCFW recursion relations [89,90] within the spinor helicity formalism. In thebbggH case, we choose to shift the momenta of gluons 3 and 4, while in thebbqqH case we choose particles 1 and 4 to avoid shifting the momenta of adjacent quarks of the same flavour. Moreover, we ensure that the shifted brackets [î , |ĵ] do not belong to particles of helicities i − , j + . These choices are necessary for the validity of the recursion relations as they prevent the shifted amplitude from having poles at infinity. For thebbggH channel we obtain the following non-vanishing tree-level partial amplitudes,

Amplitude Reduction and Reconstruction
We derive the analytic form of the two-loop QCD helicity amplitudes using the framework discussed in Refs. [25,91], which combines the Feynman diagram approach with numerical sampling over finite fields and functional reconstruction techniques. We generate a set of Feynman diagrams contributing to both thebbggH andbbqqH processes using QGRAF [92], and perform diagram filtering, topology identification and colour decomposition using a collection of in-house Mathematica and Form [93,94] scripts. In the leading colour approximation there are 749 (264) Feynman diagrams contributing to the two-loopbbggH (bbqqH) QCD amplitudes, including all closed fermion loop contributions. The numerators of the loop amplitudes are then computed for each independent helicity configuration, and the spinor-helicity manipulations are carried out using the library Spinney [95]. The helicity dependent loop numerators are written in terms of spinor products ( ij , [ij]), scalar products (p i · p j , k i · k j , k i · p j ) and spinor strings ( i|k i |j], i|p 5 |j], where p i 's and k i 's are the external and loop momenta, respectively. The angle |i and square |i] brackets denote the usual holomorphic and anti-holomorphic bi-spinors built from the external massless momenta. The momenta appearing between them in the spinor strings are intended as the matrices obtained by contracting the corresponding four-momenta with the vector of the Pauli matrices, e.g. i|p 5 |j] = i|σ µ |j] p µ 5 . The loop numerators are therefore expressed as linear combinations of monomials of loopmomentum dependent objects that multiply coefficients which depend only on the external kinematics. These coefficients contain spinor products and loop momentum independent spinor strings, and such objects not only are not independent (due to momentum conservation and Schouten identity), but they are also incompatible with the finite field arithmetic framework.
In order to allow for the use of finite field sampling in the computation of the helicity amplitudes, we make use of an explicit rational parametrisation of the external kinematics. This is constructed with help of the momentum twistor [96] and spinor-helicity formalisms. To obtain a configuration for the off-shell five-particle configuration we begin with a massless configuration for six particles. While the exact parametrisation is not important, the form presented in Ref. [97] is a useful starting point. We can think of the massless process as the result of the off-shell leg decaying into a pair of massless particles. There are 8 independent variables in the six-particle process, but we can reduce the degrees of freedom by choosing one of the decay direction to be collinear with one of the other massless legs in the fivepoint kinematics. Since the momentum twistors are associated with complex momenta we consider the angle and square bracket spinors products to be independent. We can write this explicitly as follows: the six massless momenta q µ i are related to the off-shell five particle momenta p µ i by We impose the constraints q 2 q 6 = 0 and [q 2 q 6 ] = 0 to reduce the independent variables to 6. We then apply a transformation onto the following choice of variables: where tr ± (ij · · · kl) = 1 2 tr((1 ± γ 5 ) / p i / p j · · · / p k / p l ). Some features of this parametrisation are particularly convenient. For example, the only dimensionful quantity is x 1 , and tr 5 is rational. However, such a choice does break some symmetries in the problem and so different helicity configurations will have different complexities. For the processes considered here the polynomial complexity was manageable using this form. Applications to other off-shell five-particle processes may require further thought.
Having set up a rational parametrisation of the external kinematics, the helicity dependent loop numerators are constructed analytically and ready to be further processed. This is the starting point of our numerical algorithm in the finite field setup. In order to write the loop amplitude in terms of Feynman scalar integrals, we first define the integral families for the maximal topologies, which are the topologies with the maximum number of loop propagators (8 in this case). There are 15 maximal topologies for both thebbggH and bbqqH processes, and they are shown in Fig. 3 of Ref. [50]. Each diagram topology appearing in the amplitude can be mapped onto at least one of the maximal topologies. After that is done, the objects which depend on the loop momenta in the helicity dependent loop numerators are written in terms of the 11 propagators associated with the chosen maximal topology. These mapping procedures are performed numerically within the FiniteFlow framework [43]. External kinematics that involve in the mapping are already expressed in terms of momentum twistor variables. At this stage, the loop amplitude is expressed as a linear combination of scalar integrals that will further be reduced to a set of master integrals. The coefficients associated with these scalar integrals are functions only of the external kinematics in the form of momentum twistor variables.
The integrand is now ready to be reduced to a set of master integrals using the Integration-by-Parts (IBP) reduction method [98] within the finite field setup. The IBP relations are generated using LiteRed [99] and solved numerically over finite fields using the Laporta algorithm [100]. The IBP reduction is performed directly to a set of master integrals with uniform transcendental (UT) weight identified in Ref. [47]. The UT master integrals are further decomposed into a basis of special functions f that is built out of the UT master integral components as proposed in Ref. [50]. We refer to this basis of special functions as the master integral function basis. Subtracting the UV and IR poles from the bare helicity amplitudes and performing Laurent expansion in the dimensional regulator = (4 − d)/2, we arrive at the following expression for the finite remainders, where m i (f ) are monomials of the special functions and r i are rational functions of the momentum twistor variables x. It is important to note that the definition of the UT master integrals involve three square roots. One is related to the pseudo-scalar invariant tr 5 , which captures the parity-odd part of the spinor expressions and is therefore present already in the rational coefficients. This square root is rationalised by the momentum-twistor inspired parameterisation (3.2). The other two square roots are associated with the three-mass triangle Gram determinants, where p ij = p i + p j . They are not rationalised by our parameterisation and may therefore be problematic for the finite field setup. We overcome this issue by absorbing the three square roots in the definition of the UT master integrals, which is possible because they appear only as overall normalisations of the latter. As a result, they are contained in the monomials m i (f ) in Eq. (3.3) and do not appear in the amplitude reconstruction. We note that in our computation, while the tr 5 originating from the UT master integral definition is absorbed in m i (f ), the tr 5 already present in the rational coefficients is rationalised and present in the amplitude reconstruction. At this stage we have an algorithm which computes the coefficients of the special function monomials, r i (x) in Eq. (3.4), numerically over finite fields for each of the independent helicity configurations of the two processesbbggH andbbqqH. We emphasise that from the start of our numerical algorithm until the evaluation of r i (x), the computation is done within one system of FiniteFlow graphs. The last step which remains to be done is the functional reconstruction of the rational coefficients. This task is made challenging by the high polynomial degrees, which are shown in the third column of Tables 1 and 2. We tackle the complexity of the reconstruction through the strategy already used in Refs. [25,50]. We refer to Ref. [25] for a thorough discussion, and give here only a brief outline.
The first step of the strategy consists in fitting the linear relations among the rational coefficients. These linear relations are then used to express the rational coefficients in terms of a set of linearly independent coefficients, which are chosen to have the lowest possible degrees. The degrees of the independent rational coefficients are given in the fourth column of Tables 1 and 2.
The second step of our reconstruction strategy consists in determining the factors in the denominators of the rational coefficients. The analytic structure of the special functions is determined by a set of algebraic functions of the kinematics called letters. In other words, the singularities and branch cuts of the special functions are located on the hypersurfaces where any of the letters vanishes (or goes to infinity). It is therefore natural to expect that the rational coefficients which multiply the special functions should feature poles which are similarly linked to the letters, given for this problem in Ref. [47]. Inspired from the letters we can therefore guess the factors in the denominators of the rational coefficients. Our guess is given by where the indices i, j in the spinor expressions vary from 1 to 4, while the indices i, j, k, l in the Mandelstam invariants vary from 1 to 5. The various free indices in each of the factors are understood to be different from each other. We then match an ansatz made of the factors in Eq. (3.5) against the rational coefficients reconstructed on a random univariate slice of the phase space. This allows to determine entirely the denominators, as well as some factors in the numerators.
Having determined the denominators of the rational coefficients we can proceed to the third and last step of our reconstruction strategy: a univariate partial fraction decomposition on the fly. We find that, with the parameterisation of the external kinematics given b bggH helicity configurations  by Eq. (3.2), partial fractioning with respect to x 5 is the most convenient choice. The partial fraction decomposition of multivariate functions has recently drawn a lot interest as a powerful tool to simplify the analytic expressions of scattering amplitudes. A number of new approaches have been proposed, which make use of algebraic geometry techniques to handle the multivariate case efficiently [17,[101][102][103][104]. These algorithms are however applied to known analytic expressions, namely after the rational reconstruction has been performed. Conversely, the univariate partial fraction decomposition can be performed within the finite field setup to simplify the reconstruction of the rational coefficients. The algorithm makes use of the knowledge of the denominators to construct an ansatz for the decomposition into partial fractions with respect to the chosen variable. The coefficients of this ansatz can then be solved for and reconstructed within the finite field workflow. In order to reconstruct the coefficients of the partial fraction decomposition we perform a further matching of the factors from the ansatz (3.5). The determined factors are then removed, and the remaining functions are reconstructed with FiniteFlow's reconstruction algorithm. The degrees of the rational coefficients which remain to be reconstructed are shown in the fifth column of Tables 1 and 2. Note that after partial fractioning the denominators are not entirely determined from the ansatz (3.5). The univariate partial fraction decomposition in fact introduces spurious factors in the denominators. The latter could be determined as well, but we refrain from doing so as it does not reduce the complexity of the reconstruction.
Following the strategy outlined above we reconstructed the partial finite remainders for b bqqH helicity configurations  the independent mostly-plus helicity configurations of the subprocessesbbggH andbbqqH.
Reconstruction data are shown in Tables 1 and 2. The mostly-minus helicity finite remainders can be obtained by parity conjugation. Moreover, the finite remainders for the helicity configuration + + −+ can be obtained by swapping 1 ↔ 2, 3 ↔ 4 in the + + +− finite remainders, as discussed in Section 2.1 for the tree-level amplitudes. This symmetry follows from the colour structure and thus holds at any loop order. We could therefore evaluate the + + −+ finite remainders by permuting numerically the + + +− ones, as we do in order to get all the other helicity configurations from the independent ones. By permuting numerically we mean that we obtain the values of the permuted finite remainders by evaluating the un-permuted ones at permuted points. This is possible because our approach to the evaluation of the special functions, which we discuss in Section 4, handles the analytic continuation to any region automatically. Each numerical permutation therefore amounts to one more evaluation for each point. The permutation that takes from + + +− to + + −+ is however peculiar, as it is covered by the basis of special functions defined in Ref. [50]. We can thus reconstruct the finite remainders for both helicity configurations directly in terms of the same basis of special functions. This is much more convenient, as it reduces the number of permutations which need to be carried out numerically, this way decreasing the global evaluation time of the finite remainders. For this reason we reconstructed the analytic expression of the + + −+ finite remainders as well, and include the + + −+ configuration in the independent helicity set in the following sections. The relation with the + + +− configuration constitutes a non-trivial check of our results, which we discuss in Section 5.

A Custom Basis of Special Functions for the Finite Remainders
The one and two-loop finite remainders are expressed as combinations of rational coefficients -functions of the momentum twistor variables (3.2) -and monomials of square roots and elements of the master integral function basis {f (w) i }. The latter were classified in Ref. [50] so as to span the cyclic permutations of the planar five-particle integrals with one massive off-shell leg up to two loops. The function space of the finite remainders is however simpler than that of the integrals and of the amplitudes. This becomes particularly clear when we express the special functions in terms of Chen's iterated integrals [105] (see the notes [106] for a thorough discussion of their properties). Given a set of logarithmic integration kernels {d log W i }, where the arguments W i are algebraic functions of the external kinematics, Chen's iterated integrals can be defined iteratively as where we denote cumulatively by s the dependence on the external kinematics, s 0 is an arbitrary reference point, and γ is an arbitrary contour in the phase space going from γ(0) = s 0 to γ(1) = s. The iteration starts with [] s 0 (s) := 1, and the number of iterated integrations -n in Eq. (4.1) -is called transcendental weight. At two loops up to order 0 iterated integrals with transcendental weight up to four are required. In our case, the integration kernels are given by the letters {W i } of the alphabet defined in Ref. [47].
Although they may look rather abstract, Chen's iterated integrals offer several important advantages. Most importantly, they implement automatically all the functional relations which would otherwise be hidden in a representation in terms of other types of functions, such as Goncharov polylogarithms. This property was exploited in Ref. [50] to construct the master integral function basis {f The expression in terms of Chen's iterated integrals of the master integral basis functions can be obtained through their definition in terms of master integrals components given in Ref. [50], and the canonical differential equations for the master integrals given in Ref. [47] (see Ref. [50] for a thorough discussion). Another important benefit of using iterated integrals is that their analytic structure is beautifully manifest: an iterated integral may have singularities or branch points only where one of its letters vanishes or diverges. Expressing the special functions in the finite remainders in terms of Chen's iterated integrals therefore highlights their analytic structure. Indeed this unveils important simplifications: certain letters, which are present in the expressions of the master integrals, are absent in the finite remainders. In other words, certain branch cuts of the integrals drop out of the finite remainders. We observe the same pattern noticed for the W bb amplitudes in Ref. [50]. One letter, W 49 = tr 5 , is present in the bare amplitudes but absent in the finite remainders. This pattern of the pseudo-scalar invariant tr 5 has already been observed explicitly in many massless two-loop five-particle finite remainders [15][16][17][18][19][20][107][108][109][110][111], and is linked to an underlying cluster algebra structure [112]. In addition, six letters, W i with i = 16,17,27,28,29,30, are present in the two-loop integrals but drop out already at the level of the bare amplitudes (truncated at order 0 at two loops). It would be of great interest to uncover the physical principle underlying these simplifications.
The simplifications of the analytic structure mentioned above require the interplay among various elements of the master integral function basis {f (w) i } of Ref. [50]. In other words, the separate terms of the finite remainders may have spurious branch cuts which cancel out in the sum. It is therefore convenient to construct a new, ad hoc basis of special functions where the properties of the finite remainders are manifest. In addition to being more elegant from the theoretical point of view, such a basis is also much more convenient from the practical point of view. Evaluating the master integral function basis {f (w) i } in fact requires handling integration kernels -letters -which eventually cancel out in the objects we are interested in evaluating. It is desirable that these cancellations take place analytically rather than numerically, so that the spurious kernels are avoided altogether. For this reason we need a new basis of special functions where these properties are manifest, and an approach to evaluate it which by-passes the master integral function basis {f (w) i } and its unnecessary complexity.
We thus construct a new basis of special functions, which we label by {h  i } of Ref. [50] by identifying relations among the higher weight functions and products of lower weight ones which were originally missed. We achieved this following the approach of Ref. [50], but evaluating the boundary values with higher accuracy.
In order to evaluate numerically the finite remainder function basis {h (w) i } we apply the method of the generalised power series expansion [87]. This approach has already found several successful applications to the evaluation of master integrals from the differential equations they satisfy [47,52,87,[113][114][115][116]. Here, following Ref. [50], we apply it directly to the basis of special functions. The main advantage is that having a basis of special functions allows to subtract from the bare amplitudes the IR and UV poles analytically, which makes it possible to reconstruct the finite remainders directly over finite fields.
The method of the generalised power series expansion can be applied to the finite remainder and master integral function basis because they too, like the master integral bases they stem from, satisfy systems of differential equations in the canonical form [34]. This follows from the fact that the functions in the bases {h for some rational constant coefficients c (i) I 1 . Second, the derivative of a pure weight-w func-tion is a uniform transcendental function with weight w − 1. This follows straightforwardly from the differential of Chen's iterated integrals, where we note that the derivatives of log W i are algebraic functions and hence have transcendental weight 0. The transcendental purity of the functions in the basis therefore implies that the vector of all weight-w functions in the finite remainder basis, 2 , . . . T , satisfies a differential equation of the form for w > 1, where b j are constant rational matrices. The matrices b (w) j can be determined easily from the iterated integral expression of the functions through Eq. (4.3). In general the derivatives of the weight-w functions may involve weight-(w − 1) functions which are not needed to express the finite remainders. These must be included in the basis as well, in order to be able to write down the differential equations. So, strictly speaking, the finite remainder function basis {h where is an auxiliary parameter with transcendental weight −1, so that the vector h is pure with transcendental weight 0. The full basis h then satisfies a system of differential equations in the canonical form [34], where the constant matrices a i are strictly upper block triangular, with the blocks given by the matrices b (w) i in the differential equations (4.4) for the various weights. This system of differential equations is however much simpler than that for the master integrals. First, it contains only the letters which do appear in the finite remainders, i.e. a i = 0 for i = 16, 17, 27, 28, 29, 30, 49. Second, while the differential equations for the master integrals contain information about all the orders in , the system (4.6) encodes only those orders which are relevant for the finite remainders, i.e. up to 4 . The constant matrices a i are in fact nilpotent with degree 5, i.e. (a i ) 5 = 0 for any i, which follows from their strictly upper triangular form.
We neglect them here to simplify the presentation.

We evaluate the finite remainder function basis {h (w)
i } by solving the system of differential equations (4.6) numerically with the method of the generalised power series expansions [87]. For this purpose we make use of the public Mathematica package Dif-fExp [88]. In order to compute the boundary values, i.e. the values of the functions at some base point, we relate the finite remainder functions to master integral components using the definition of the master integral function basis [50], and evaluate the latter through their analytic expression in terms of Goncharov polylogarithms [117][118][119] provided in Refs. [8,48,49]. We evaluate the Goncharov polylogarithms numerically with the C++ library GiNaC [120]. We provide in ancillary files the differential equations (4.6) and the values of the finite remainder functions at six points, one for each 2 → 3 scattering region with massless incoming particles, with (at least) 200-digit accuracy. Using this information, the generalised power series expansion method allows to evaluate the finite remainder function basis reliably anywhere in the phase space. This technique in fact handles the analytic continuation automatically, so that also the permutations required to evaluate the complete finite remainders starting from the partial ones -as shown in Eqs. (2.7) for the amplitudesas well as the other helicity configurations can be performed numerically straightforwardly. We assessed the reliability of this evaluation approach by checking the convergence of the finite remainders close to the spurious poles of the rational coefficients. We discuss this and other checks in the next section.

Further Validation
The finite remainders are defined by subtracting the expected UV and IR poles from the bare amplitudes. Therefore, the fact that all the poles in cancel out, so that our expressions for the finite remainders are indeed finite, is already a strong check of the validity of our results. We have also checked that our amplitudes are correctly normalised by comparing the spin and colour averaged squared tree-level amplitudes with full colour dependence against MadGraph [121]. On top of that, we performed a number of additional checks, which we discuss in the next subsections.

Direct computation of the squared finite remainders
We checked the helicity amplitudes derived in this paper against a squared matrix element computation, carried out independently following the approach used in the previous work on ud → W bb [50]. In the squared matrix element computation the two-loop amplitude is interfered with the tree-level amplitude analytically. After manipulating the Dirac traces, the loop numerators contain only scalar products (k i · k j , k i · p j , p i · p j ), that can be mapped easily onto propagator denominators. This way we derive an analytic form of the twoloop squared matrix element written as a linear combination of scalar Feynman integrals. The squared matrix element is then processed further through IBP reduction to obtain the master integral representation, followed by mapping of the master integrals onto the basis of special functions, subtraction of UV and IR singularities, and finally Laurent expansion in . All these steps are carried out numerically over finite fields within the FiniteFlow framework. The squared matrix element computation uses the Conventional Dimensional Regularisation (CDR) scheme, where internal and external live in d = 4 − 2 dimensions. We find complete agreement with this approach and the helicity amplitudes in the tHV scheme at the level of the squared finite remainders.

Renormalisation scale dependence
The analytic expressions of the one-and two-loop finite remainders have been derived with the renormalisation scale µ set to µ = 1. The dependence of the finite remainders on the renormalisation scale can be determined by restoring the µ dependence of the strong and Yukawa couplings (α s → α s µ 2 , y b → y b µ ), and by replacing (−s ij ) − → (−µ 2 /s ij ) in the I 1 operators defined in Eqs. (2.13) -(2.14), which enter the pole terms in Eqs. (2.10) - (2.11). In order to capture the scale dependence of the finite remainders we define the difference where the dependence on the kinematic variables is understood. The difference δF (L),i (µ 2 ) is entirely determined by the finite remainders evaluated at µ 2 = 1 -which we reconstructed analytically -and logarithms of µ 2 . For thebbggH finite remainders it is given by while for thebbqqH finite remainders it is given by To check that our results for the finite remainders have the correct scale dependence, we evaluate them at two kinematic points related by a rescaling by some positive factor a, s = (s 12 , s 23 , s 34 , s 45 , s 15 , p 2 5 ) , s = a s = (as 12 , as 23 , as 34 , as 45 , as 15 , ap 2 5 ) .

(5.12)
We then verify numerically that the finite remainders satisfy the following relation, where we extended the notation of the finite remainders and tree-level amplitudes to take into account their dependence on the phase space point s.

Relation between + + +− and + + −+
The partial finite remainders for the single-minus helicity configurations, + + +− and + + −+, are related by a permutation of the external particles, (5.14) Permuting the special functions is however non-trivial, and may in general require analytic continuation. The permuted special functions are defined by where σ denotes the permutation (12345) → (21435) of the external particles, and {p i } denotes the dependence on the external momenta. In order to check the relations (5.14) analytically, we need to express the permuted functions (σ • h) in terms of the ones in the standard orientation, h. To perform this operation in an algorithmic and robust way we resort to the system of differential equations (4.6) satisfied by the finite remainder function basis. The permuted functions in fact satisfy the permuted differential equations Permuting the letters W i is however trivial, as they only involve algebraic functions. Since the alphabet is closed by construction under this permutation, we obtain an explicit system of differential equations for the permuted special functions, The latter can straightforwardly be solved in terms of Chen's iterated integrals (see Ref. [50] for a thorough discussion). In order to be comparable with the solution of the system (4.6), which defines the finite remainder function basis h, we must make sure that the same boundary point is used when solving both systems of differential equations in terms of iterated integrals. The boundary values can be obtained numerically with arbitrary precision using the expressions in terms of Goncharov polylogarithms of Refs. [8,48,49], as discussed at the end of Section 4. Using the differential equations (5.17) we can then express the permuted finite remainder special functions in terms of the ones in the original orientation through Chen's iterated integrals. Once that is done, the right-hand sides of Eqs. (5.14) are written in terms of the same special function basis as the left-hand sides. Since the rational coefficients can be permuted analytically trivially, it is immediate to verify that our results for the one-and two-loop finite remainders satisfy the relations given by Eqs. (5.14).

Convergence near spurious poles
The rational coefficients in the finite remainders contain spurious poles, namely poles which are not related to the physical singularities of the amplitudes. When evaluating at a phase space point infinitesimally close to a spurious pole -but at a finite distance from all physical poles -the rational coefficients diverge, whereas the finite remainders must stay finite. This can only occur through large numerical cancellations involving both the rational coefficients and the special functions. Verifying this behaviour constitutes both a check on our analytic expressions for the finite remainders and a stress test of our evaluation approach for the special functions. We do it as follows. First, we use the factor matching strategy discussed in Section 3 with the ansatz given by Eq. (3.5). This allows us to determine that the spurious poles in the rational coefficients are associated with the following factors in their denominators, 18) where i, j ∈ {1, 2, 3, 4} with i = j. Next, for each factor in Eq. (5.18), we construct a one-dimensional slice of the phase space, parametrised by a parameter δ, such that P i = δ. As we probe the small-δ region all the other factors in Eq. (5.18) and the factors associated with the physical singularities -(p i +p j ) 2 and p i ·p j -must stay finite, i.e. they must neither diverge nor vanish. This ensures that we target a specific spurious pole, rather than multiple at once, and that we stay away from the physical singularities. Finally, we evaluate the finite remainders on these slices for increasingly small values of δ. We evaluate the special functions with 64-digit accuracy. Our analytic expressions for the finite remainders exhibit the expected behaviour: as δ approaches zero the rational coefficients diverge, while the finite remainders converge to finite values.

Results
We provide the analytic expressions of the independent one-and two-loop finite remainders in the ancillary files. We present them as combinations of linearly independent rational coefficients and monomials of the square roots and of the finite remainder basis functions.
The rational coefficients are expressed in terms of the momentum twistor variables defined in Eq. (3.2). We evaluate the finite remainder function basis numerically by integrating the differential equations (4.6) with the method of the generalised series expansions. We also provide Mathematica scripts which illustrate how to evaluate the finite remainders interfered with the tree-level amplitudes for all the partonic channels contributing to the process pp → bbH, which we label as The interference between the finite remainders and the tree-level amplitudes summed over colour and helicity is given at leading colour by colour helicity where α = 3 for gg and α = 2 for all the other channels, and the reduced squared finite remainders, H (L) , are given explicitly for each channel by We evaluate the permutations of the finite remainders numerically. The analytic continuation is performed by adding a small positive imaginary part to each s ij and to p 2 5 , which is done automatically by DiffExp.
To facilitate future checks, we provide here the benchmark evaluation at a physical point corresponding to the pp → bbH scattering process, The values of the bare two-loop amplitudes normalised by the tree-level amplitudes,Â (L) = A (L) /A (0) , for the independent mostly-plus helicity configurations of the subprocesses 0 → bbggH and 0 →bbqqH are given in Tables 3 and 4. Table 5 shows the values of the two-loop reduced squared finite remainders H (2) . Analogous tables for the one-loop amplitudes are given in Appendix B.
To demonstrate the suitability of our results for phenomenological applications we present the evaluation of the finite remainders interfered with the tree-level amplitudes on a univariate slice of the phase space. For this purpose we parametrise the momenta for the   Table 3: Numerical values of the barebbggH partial amplitudes at two loops (normalised to the tree-level amplitude) at the kinematic point in Eq. (6.11) for the four independent helicity configurations and the various closed fermion loops contributions.
scattering processes (2.1) in terms of angles and energy fractions of the final state as   Table 5: Numerical values of the two-loop reduced squared finite remainders H (2) defined in Eqs. (6.3)-(6.9) at the kinematic point in Eq. (6.11) for the various closed fermion loops contributions and the scattering channels specified in Eq. (6.1).
while p 5 is fixed by momentum conservation, p 5 = −p 1 − p 2 − p 3 − p 4 . We have chosen the particle with momentum p 1 to be produced at an angle of π/2 with respect to the beam axis. Requiring that the Higgs boson is on-shell, p 2 5 = m 2 H , constrains the angle θ as To restrict the kinematics to a one-dimensional slice we choose 14) The reality of the angle θ then restricts the free parameter y 2 to the interval y 2 ∈ [ 39 100 , 39 40 ]. In order to evaluate the finite remainders for all the processes shown in Eq. (6.1) we need to compute the finite remainder special functions at 16 permutations of each phase space point, as can be seen explicitly in Eqs.(6.3) -(6.9). We do this by integrating the system of canonical differential equations (4.6) with DiffExp using rationalised values of the invariants. For each permutation of each point we compute the "distance" (technically speaking, the number of segments into which the integration contour is divided by DiffExp) from the six reference points provided in the ancillary files, and choose the closest one as initial point for the integration. In Figure 3 we plot the values of the one-and two-loop reduced squared finite remainders for all the processes defined in Eq. (6.1), as functions of the parameter y 2 .
We stress that the purpose of the plots in Figure 3 is merely to demonstrate that our results for the finite remainders can be evaluated reliably in the physical scattering region. Nothing can be inferred about the convergence of the perturbative expansion at the cross section level. One interesting feature which can be appreciated from Figure 3 is the appearance of a loop-induced peak in the finite remainders for the channelqq. The peak is absent at tree level for the same channel and up to two loops for qq. The latter channel is related toqq by the exchange 3 ↔ 4 of the external particles. We observe that the peak stems from the values of the finite remainder function basis, while the rational coefficients are not enhanced in that region. In order to pinpoint more precisely the origin of this phenomenon, we construct the explicit analytic expressions for some of the special functions which exhibit the peak, starting from their iterated integral expression obtained by solving the system of canonical differential equations (4.6). For instance,     qq for L = 1, 2, is peaked at y 2 ≈ 0.5566, where s 24 reaches its minimum absolute value on the slice. The tree-level amplitudes for the channel gg instead do have poles at s 23 = 0, which can be seen explicitly in Eqs. (2.21). Since H (L) gg with L = 1, 2 receive contributions from the partial finite remainders in both the standard orientation and with the swap 3 ↔ 4, as shown in Eq. (6.3), their plot in Figure 3 (a) exhibits this peak already at tree level. The same holds for the bb andbb channels, as can be seen in Figures 3 (d) and (f). Also in Figure 3, we observe divergences at y 2 = 39/100 for the processes bb andbb. This divergence is associated with the propagator 1/s 12 , which can only appear in processes with two pairs of bottom quarks. In Eqs. (6.6) and (6.7) we can see the qq fermion pairs can appear with momenta p 1 and p 2 , which is not the case for other processes. All the other features of the plots in Figure 3 can be similarly understood in terms of tree-level propagators going on shell.

Conclusions
In this article we have presented an analytic form for the two-loop QCD corrections to the process pp → bbH. To the best of our knowledge this is the first complete set of helicity amplitudes provided for a 2 → 3 scattering process with an off-shell leg. A special function basis for the finite remainder was identified obeying a canonical form differential equation. The method for constructing such a basis was presented in the recent pp → W bb computation [50] by three of the present authors. The finite remainders were then extracted from multiple evaluations over finite fields using a rational phase-space and IBP reduction. We obtained relatively compact results after determining the linear relations between the rational coefficients and performing a univariate partial fraction decomposition on the fly. The final expressions were evaluated using the method of generalised series expansions as implemented in the DiffExp code [88].
The expressions have been validated in a number of ways and we observe that they exhibit a smooth behaviour in all scattering regions. Evaluation times appear to be suitable for phenomenological applications, especially since we have not tried to optimise the route through the phase-space evaluations as has been done in other applications of the method [47,52,[114][115][116].
The techniques presented here show promise for applications to other important scattering processes such as pp → V + 2j and pp → H + 2j, although the non-planar sector is not entirely known and remains a high priority.
The analytic expressions together with Mathematica scripts to evaluate them numerically are provided in the ancillary files accompanying this article.

A Renormalisation Constants
The one-and two-loop bottom-quark Yukawa renormalisation constants entering the UV counterterm in Eqs. (2.10) and (2.11) are while the β function coefficients and anomalous dimensions are

B One-Loop Results
We present in Table 6 the numerical values of the one-loop bare helicity amplitudes at the kinematic point given in Eq. (6.11), evaluated through O( 2 ) for the different closed fermion loop contributions defined in Eq. (2.8). These one-loop amplitudes are required in the computation of the two-loop pole terms in Eqs. (2.10) and (2.11). Furthermore, Table 7 shows the tree-level reduced squared amplitudes and one-loop reduced squared   Table 7: Numerical values of the tree-level reduced squared amplitudes H (0) and one-loop reduced squared finite remainders H (1) defined in Eqs. (6.3)-(6.9) at the kinematic point in Eq. (6.11) for the closed fermion loops contributions and the scattering channels specified in Eq. (6.1).