One-loop weak corrections to Higgs production

We compute mixed QCD-weak corrections to inclusive Higgs production at the LHC from the partonic process gg→Hqq¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ gg\to Hq\overline{q} $$\end{document}. We start from the UV- and IR-finite oneloop weak amplitude and consider its interference with the corresponding one-loop QCD amplitude. This contribution is a Oαsα\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}_s\alpha \right) $$\end{document} correction to the leading-order gluon-fusion cross section, and was not numerically assessed in previous works. We also compute the cross section from the square of this weak amplitude, suppressed by Oα2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}^2\right) $$\end{document}. Finally, we consider contributions from the partonic process gq → Hq, which are one order lower in αs, as a reference for the size of terms which are not enhanced by the large gluon luminosity. We find that, given the magnitude of the uncertainties on current state-of-the-art predictions for Higgs production, all contributions computed in this work can be safely ignored, both fully inclusively and in the boosted Higgs regime. This result supports the approximate factorisation of QCD and weak corrections to that process.


Introduction
In the quest towards an ever more accurate prediction for the inclusive Higgs production cross section at hadron colliders, one of the major tasks is the computation of fixed-order corrections in the perturbative expansion in powers of the Standard Model (SM) couplings. Our understanding of pure QCD corrections, which are known to be very important for this process, has reached an unprecedented level of accuracy in recent times. A milestone in this programme was achieved with the computation of the third correction term in the expansion in the strong coupling α s of the cross section for Higgs production via gluon fusion in the infinite top mass limit [1,2]. In a typical setup for the LHC running at a centre-of-mass energy of 13 TeV, this contribution shifts the prediction for the total cross section upwards by roughly 3% [3].
On the other hand, weak corrections to the leading-order (LO) inclusive Higgs cross section also need to be considered. In the same setup mentioned before, the first weak term turns out to increase the total gluon fusion cross section by a significant 5% [4][5][6]. Since next-to-leading-order (NLO) QCD corrections can be as large as the leading contribution, the motivation to investigate mixed first-order QCD and first-order weak corrections is very strong. Although the exact size of this term is at present unknown, various approximations have been considered in the literature. The first estimate to appear was based on the JHEP05(2019)002 argument that mixed QCD-weak effects on the inclusive Higgs production cross section are well approximated by combining the purely weak term and the full QCD series in a multiplicative fashion [7]. Following this factorisation approach, the authors of ref. [3] reported the mixed QCD-weak corrections to be approximately 3% of the full result, and conservatively estimated the uncertainty stemming from non-factorisable contributions to be 1% of the total. The estimates of [3,7] are obtained by considering the unphysical limit m H m W , m Z . The gluon induced interference contributions discussed in our work are suppressed in this limit by two powers of the weak boson masses with respect to the leading order O(α 2 S α) cross section, which we verified by explicit calculation. The theoretical uncertainty associated to each of the other main error sources (determination of parton distribution functions, truncation of the QCD perturbative series, and missing quark-mass effects) is currently of the same order. It is therefore highly desirable to remove the ambiguity due to the factorisation approximation.
Important steps have recently been made in this direction. Thanks to the calculation of the three-loop mixed QCD-weak correction to Higgs boson gluon fusion for arbitrary masses of the W , Z, and Higgs bosons [8], an estimate of the cross section in the softvirtual approximation was obtained [9]. An independent work considered three-loop matrix elements in the limit of massless vector bosons instead, and combined them with a different class of two-loop real-emission contributions [10]. The estimates obtained using these approximations support the validity of the factorisation approach, since they include some non-factorisable effects and find that these are numerically small.
In order for the full mixed QCD-weak term to become available, however, two pieces of the puzzle are still missing. On the one hand there is the formidable challenge of computing two-loop matrix elements with an extra real emission for arbitrary W , Z, and Higgs masses. On the other hand, there are UV-and IR-finite one-loop weak contributions to the production of the Higgs in association with two partons, which feature more complicated kinematics but whose one-loop integrals are well understood. Although in general corrections with fewer or soft real emissions are expected to dominate within the inclusive cross section [9], the contributions with two extra hard partons are formally of the same order and may disrupt the approximate factorisation of weak and QCD corrections because of their final-state kinematic structure.
In the present paper, we address this issue by carrying out the exact inclusive computation of the contribution to mixed QCD-weak corrections from the one-loop partonic subprocess gg → Hqq. We stress that this contribution features one-loop pentagon topologies which appear only in matrix elements with (at least) two real emissions, that do not fit in a factorised picture and that have not been assessed before.
The paper proceeds as follows. In section 2, we discuss the different contributions that enter our computation, we categorise them and identify potential competing mechanisms which are formally of the same order or slightly higher. Although the computation of the required matrix elements is straightforward using standard public codes for one-loop calculations, the computation of the pieces of cross sections we are interested in requires the renormalisation of parton distributions and the subtraction of initial-state collinear singularities. Given the very special features of the process examined, these steps require JHEP05(2019)002 some care and are thus described in section 3. Finally, we report and discuss numerical results.

Classification of contributions
In order to classify contributions to the Higgs inclusive cross section, it is useful to write its mixed QCD and weak expansion as where the prefactor α 2 s α 1 is chosen so as to match the couplings factorised by the leadingorder loop-induced gluon-fusion contribution to inclusive Higgs production. Notice that we group all squared couplings that are not strong, including the Yukawa of the top quark, under the label α, in view of their comparable strength and of the electroweak gauge relations often rendering their separate factorisation ambiguous. The corrections often labelled "QCD N m LO" and "(electro)weak N n LO" are then denoted by σ (m,0) pp→H+X and σ (0,n) pp→H+X , as they become impractical when addressing the mixed cases σ (m,n) pp→H+X . With such a notation in mind, the expected naive parametric suppression from the couplings, which counts α s ∼ 10 −1 and α ∼ 10 −2 , simply reads σ (m,n) pp→H+X ∼ 10 −m−2n . In order to discuss interference terms, we also find it useful to introduce a similar notation for amplitudes: where we denote by g all couplings that are not g s . As mentioned above, weak and QCD corrections are expected to factorise to a certain degree, such that This approximation is valid under the assumption that the main contributions to the mixed QCD-weak cross section are to be attributed either to soft gluons or Sudakov weak logarithms. If one is to assess violations of this factorisation, the expansion term σ (1,1) pp→H+X must be computed exactly. We now set out to discuss the many contributions this term receives.
In this work, we only consider weak corrections involving the W and Z bosons, as these dominate over the genuine electroweak corrections (i.e. unresolved photon exchange or emission) to contributions where the Higgs is produced from massive quark loop lines that are not the top-quark.
Also, the gluon initiated processes are expected to be the dominant contributions at the LHC, where quark parton distribution functions (PDFs) are small in comparison to the gluon one for the typical values of the Bjorken x's probed by the kinematics involved. We therefore neglect all contributions to σ (1,1) pp→H+X that factorise parton luminosities with The results reported in this work are highlighted with a green background, while those addressed in refs. [7,9] are denoted in blue. Together, these form the complete σ at least one quark. To get a reference for the size of these terms that we do not compute, we report numerical results also for σ (0,1) gq→Hq . 1 Weak corrections stemming from the interference with leading QCD production modes are often subject to kinematic suppressions that renders them smaller than what is naively expected from their factorised couplings. For this reason, we also report the pieces of the cross sections σ gq→Hq . These form a gauge-invariant subset of higher-order contributions.
Our work reports on the contribution σ (1,1) gg→Hqq for the first time and, together with the results from refs. [7,9], it completes the computation of σ (1,1) gg→H+X . We now proceed to list in table 1 all amplitudes building σ (1,1) gp→H+X . We now turn to discussing the Feynman diagrams building the amplitudes A gq→Hq that contribute to the cross sections presented in this work. The amplitude A (2,0) gg→Hqq is built from the diagrams depicted in figure 1 where the Higgs is produced via weak vector boson fusion and interfered with the leading QCD gluon-fusion diagrams shown in figure 3.
Diagrams of the class 1d and 1e, where the Higgs is produced via gluon-fusion, feature a Z-boson propagator 2 which however does not yield any Breit-Wigner resonance as they JHEP05(2019)002 gg→Hqq . In diagrams 1a, 1b and 1c, the Z boson can be interchanged with a W boson. Diagrams 1d and 1e are contributions to the production of a Higgs in association with a Z boson and are only included in the computation of σ (1,1) pp→H+X , and not that of σ (1,2) pp→H+X . Diagrams of the class 1f are only present for the process gg → Hbb.
are interfered against the QCD diagrams of figure 3. We must nonetheless regulate the Z-boson propagator pole, which motivates our use in this computation of the complexmass scheme [11,12] with finite widths for the internal top quark and unstable weak gauge bosons. These diagrams 1d and 1e are however ignored when considering their squared contribution to σ (1,1) pp→H+X , since in this case they are best accounted for in the narrowwidth approximation as the LO prediction for associated Higgs production, i.e. σ (1,1) gg→HZ (also reported in this work).
Finally, diagrams of the class 1f are specific to the third-generation quarks where the Higgs can also be emitted from the top-quark running in the loop. This contribution is analogous to that of the heavy quarks in the two-loop electroweak corrections to Higgs production investigated in ref. [13] and, for this reason, we found it interesting to report our results separately for the processes gg → Hqq, with q ≡ u, d, c, s, and gg → bbH.

Initial-state collinear singularities
All of the one-loop amplitudes considered in this paper are free of explicit ultraviolet and infrared divergences that can arise from the integration over the loop momenta. In other words, working in dimensional regularisation with D ≡ 4 − 2 , their analytic expressions do not contain explicit poles in the dimensional regulator . However, matrix elements may feature non-integrable infrared divergences in regions of the phase space which correspond to unresolved configurations. In order to discuss this issue, we concentrate on the amplitude A (0,2) gg→Hqq as it constitutes the main focus of the present work.
In principle, the process gg → Hqq presents infrared divergences when the quarkantiquark pair in the final state is collectively soft, and/or when one or both of the quarks JHEP05(2019)002 gq→Hq respectively. Diagrams 2a and 2b also appear in the reduced matrix elements factorised by the collinear subtraction local counterterms of eqs. (3.2) and (3.1). Notice that diagrams belonging to the class 2d are specific to the process gb → Hb. In all cases, the full top-quark mass dependence is retained. qq→H (that is, of order g 3 ) which is identically zero. Indeed, the triangle one-loop diagrams for qq → H require a mass insertion for the chirality flip and therefore vanishes for massless onshell quarks. This explains why the interference involving the amplitude A (0,2) gg→Hqq only requires the subtraction of single-unresolved infrared limits, while the interference built upon the amplitude A (−1,2) gq→Hq does not require IR subtraction at all.
The same observations can be made, perhaps more intuitively, by inspecting the representative Feynman diagrams depicted in figure 1. It is straightforward to see that propagators of massless partons which do not belong to closed loops can go on-shell only in the graphs of type 1b and 1c. In the case of the diagram 1b, this happens when antiquark d 5 becomes collinear to gluon g 2 such that the hard scattering subgraph corresponds to diagram 2a. By contrast, in the kinematic limit where quark d 4 is collinear to gluon g 1 and quarkd 5 is collinear to gluon g 2 , both non-loop propagators of graph 1c are singular. The subgraph that describes the hard scattering process, however, evaluates to zero for massless quarks as explained before, thus avoiding the singularity. In the limit where only one of the quarks is collinear to an incoming gluon, the hard part of diagram 1c matches that of graph 2b.

JHEP05(2019)002
From the observations drawn so far, we conclude that for the local subtraction of implicit singularities it is sufficient to consider standard NLO initial-collinear counterterms. These subtraction terms are to be added back, analytically integrated over the unresolved degrees of freedom yielding explicit poles in the dimensional regulator . These poles cancel against those part of the parton distribution function (PDF) renormalisation counterterms, as guaranteed by collinear beam factorisation, thus rendering the complete computation finite.
The formal expression which describes this subtraction procedure and the combination with PDF renormalisation counterterms reads: where the dependences on the factorisation and renormalisation scales µ F and µ R as well as on the kinematic inputs for the matrix elements have been suppressed for brevity. The sums run over the four permutations π that are obtained exchanging the quark and the antiquark in the final state and/or the two initial-state gluons among themselves. The symbol C ij denotes the local counterterm for particles i and j going collinear and C ij its counterpart analytically integrated over the unresolved degrees of freedom. The observable functions are indicated with J, and ∆ ik is the PDF renormalisation kernel for parton with flavour i to change into species k before entering the hard process. The notationφ Hq indicates reduced kinematics of lower multiplicity which are obtained by mapping a pair of collinear partons to a massless parent. The concrete expressions of all subtraction ingredients closely follow ref. [14] and are presented more explicitly in appendix A, where we also explicitly show that our subtraction counterterms correctly regulate the relevant collinear singularities.

Setup of the computation and numerical results
The amplitudes A gq→Hq that factorise a Higgs coupling to weak bosons were first computed analytically (for massless quarks only) in ref. [15], in the different context of NLO QCD corrections to weak vector-boson fusion. In the present case and as indicated in table 1, in order to obtain contributions to σ (1,1) gg→Hqq and σ (0,1) gq→Hq , these amplitudes must be interfered against their corresponding QCD analog.
Nowadays such one-loop amplitudes are readily available from many automated oneloop matrix-element generators. However, a high degree of flexibility is necessary in order to be able to select the relevant diagrams and interferences, and to construct the appropriate subtraction terms. This motivates our choice of generating the relevant one-loop squared amplitudes using MadLoop [16], part of MadGraph5 aMC@NLO [17] (henceforth abbreviated MG5aMC), as it can efficiently generate and interfere [18] [19,20] and OneLOop [21], or alternatively COLLIER [22], for performing one-loop reductions and for the evaluation of the scalar one-loop master integrals. We present in appendix C some details about the generation procedure as well as benchmark numbers in order to facilitate the reproduction of our results. Moreover, we have cross-checked MadLoop's numerical implementation of the amplitudes A gg→Hqq against a completely independent and analytical computation described in appendix B.
As already mentioned, we choose to renormalise all unstable particles in the complexmass scheme [11,12] and consider the SM input parameters given in table 2, for the LHC at a collision energy of 13 TeV.
The numerical Monte-Carlo integration as well as the necessary IR subtraction procedure, presented in eqs. (3.1) and (3.2) as well as in appendix A, have been implemented in a private extension of MG5aMC currently under development. The poles in the dimensional regulator have been checked to cancel as expected. 3 Moreover, we have validated our code by comparing NLO QCD cross sections against results from MG5aMC for the processes pp → Z and pp → H, the latter in the Higgs Effective Theory.
Our results are presented in table 3. Along with the different contributions to the inclusive cross section for Higgs production, we also report the semi-inclusive cross sections for the production of a Higgs boson with transverse momentum larger than 400 GeV. The motivation to consider this boosted Higgs regime is twofold. On one side, it mimics typical experimental selection cuts used to reduce backgrounds and study new physics effect prominent in that regime. On the other side, it selects a region of phase space where real emissions are typically hard and the relative importance of the corrections computed in this work may in principle be enhanced.
We find that the squared contributions of order O α 2 s α 3 can be suppressed compared to their O α 3 s α 2 counterpart by less than what is expected by their parametric ratio α/α s .    Table 3. Fully and semi inclusive cross sections at LHC13 obtained with SM input parameters given in table 2 for the processes gg → Hqq and gg → Hbb, as (partial) contributions to the corrections of order O α 3 s α 2 and O α 2 s α 3 to inclusive Higgs production. For the contributions of order O α 2 s α 3 labelled "no-HZ", the diagrams of the class 1d and 1e are ignored, as they are best accounted for in the narrow-width approximation as the LO contribution to the process gg → HZ, which is also shown. We also report the O α 2 s α 2 and O α s α 3 contributions from the quarkinitiated processes qg → Hq and bg → Hb. Finally, we consider the boosted regime, in which the Higgs transverse momentum is required to be at least 400 GeV. In each bracket separated by a dashed line, the upper number corresponds to the scale choice µ R = µ F = m H /2 while the lower one corresponds to µ R = µ F = m H .

JHEP05(2019)002
This is for example the case for the processes involving b quarks, and it can be explained by the kinematic suppressions interfering contributions are typically subject to.
Also, contributions of order O(α m s α n ) with m + n = 4 are numerically more relevant than those with m+n = 5 in spite of their suppression by one quark luminosity. The quarkinitiated weak corrections are however still small in comparison with the whole σ pp→H+X , and can thus be safely neglected as already observed in ref. [23]. These two observations reinforce the conclusion that the contributions to σ gg→Hqq that are computed in this work and which have been neglected up to this point are of similar (ir)relevance to that of other neglected terms of weak origin.
The cross section σ gg→Hbb from only final-state b quarks reveals that contributions featuring Higgs production from the internal top quark line (see figure 1f) are comparable and of opposite sign to that of emissions from internal weak bosons. This fact contrasts with the study of ref. [13] of the two-loop amplitude A (0,2) gg→H where it was instead found that Higgs emissions from internal top quarks only contribute to less than 2% of the complete amplitude at this order, and could thus be safely ignored in the computation of the threeloop amplitude A (2,2) gg→H of refs. [8][9][10]. Indeed, the higher partonic collision energy probed by A (1,2) gg→Hbb enhances contributions from internal top-quark Higgs emissions, even more so in the boosted regime. Similarly, the same mechanism enables the bottom-quark initiated contribution σ bg→Hb at the same level as that of the channels initiated by each other valence quark flavour.
The squared contribution σ (α 2 s α 3 ,no-HZ) gg→Hbb omits the diagrams 1d and 1e featuring a Z boson decay since it is best accounted for in the narrow-width approximation. It is however clear that the extent to which one should consider Higgs production in association with an on-shell Z boson depends on the particular observable considered. We chose to report here the quantity σ (α 2 s α 2 ,Γ Z =0) gg→HZ only to serve as an upper bound to this contribution.
Notice that the contribution σ (α 2 s α 3 ,no-HZ) gg→Hqq is negative, despite involving squared amplitudes. This originates from the finite logarithms in the PDF renormalisation term ∆ gq and integrated counterterm C gq , stemming from dimensional regularisation. Our results also highlight that considering a subset of higher-order corrections and factorising only a particular combination of initial-state flavours typically yields a large dependence on the factorisation scale. This is especially true for squared amplitude contributions and the boosted regime, for which the chosen fixed scales proportional to the Higgs mass (as it is tailored to the prediction of the inclusive Higgs production cross section) are not well suited in light of the significantly larger collision energies probed. We choose to report here absolute factorisation scale dependency, given that some contributions can be accidentally close to zero 4 for one of the two scale choices. A more detailed analysis of the sensitivity of these contributions to the factorisation scale is beyond the scope of this work.
The overall magnitude of all contributions computed here is such that they can be safely neglected in light of the total size of the gluon fusion cross section 48.58 pb +2.22 pb

JHEP05(2019)002
(theory) [3] with contributions of weak origins that are estimated to be of the order of 2 pb, with a theoretical uncertainty in the range of 200 fb [9,10]. The aggregated sum σ total of all contributions computed in this work is only meant to serve as qualitative highlight of that fact. Our results then further support the factorisation approximation when accounting for mixed weak and QCD corrections to inclusive Higgs production.
The hierarchy of the various terms is altered when considering the boosted Higgs regime, where the kinematic suppression of gluon-initiated interference contributions is strong enough to make them of the same order or smaller than their squared counterpart. All interference contributions also become negative in this case, while the square term σ (α 2 s α 3 ,no-ZH) gg→Hqq is now positive as hard real emissions become dominant. Overall, none of the contributions computed plays a significant role in that scenario either, given that the pure QCD contribution is estimated in ref. [24] to be 25 fb with a large theoretical uncertainty exceeding 20%.

Conclusion
The large QCD corrections to inclusive Higgs production at LHC13 calls for accounting for mixed weak and QCD corrections in a multiplicative scheme, that is assuming their complete factorisation. In light of the accuracy sought-after for this process, it is important to assess the validity of this factorisation assumption by explicitly computing σ (1,1) pp→H+X , namely the mixed QCD and weak correction of order O(α 3 s α 2 ) to the Higgs inclusive cross section.
To this end, two groups [8][9][10] computed σ (1,1) gg→H+X and found that it supports the hypothesis that weak corrections factorise. These works however neglected the quarkinitiated components as well as the "double-real" channel gg → Hqq and we confirm here that these terms can be safely neglected, amounting to about 5% of the total mixed weak and QCD corrections. We verified that our conclusions also apply when imposing that the Higgs transverse momentum lies above 400 GeV. The interference nature of the contributions σ (1,1) gg→Hqq and σ (1,1) gq→Hq renders them prone to kinematic suppressions, and we indeed found that the square of the one-loop weak amplitudes involved can be larger than naively expected from their parametric suppression of α/α s . The selective nature of the contributions computed in this work is such that they feature a large factorisation scale dependency, further stressing that their inclusion would require to also consider all other partonic channels.
Besides further establishing the validity of the hypothesis assumed when accounting for weak corrections to inclusive Higgs production, our work also showcases the novel flexibility brought by recent developments in the realm of automated one-loop matrix element generation and Monte-Carlo integration for higher-order computations.

JHEP05(2019)002
from the European Research Council (ERC) under grant agreements No. 772099 (JetDynamics) and No. 694712 (PertQCD). We thank the Galileo Galilei Institute for Theoretical Physics for the hospitality and the INFN for partial support during the completion of this work.

A Initial-collinear counterterms
In this appendix, we detail all ingredients that are necessary for the subtraction of implicit singularities outlined in eqs. (3.1)-3.2, and we demonstrate that the matrix element for g 1 g 2 → q 3q4 H 5 is correctly regulated in the regions of phase space close to unresolved configurations. As already announced in section 3, the construction and notation follow closely ref. [14].
Let us begin with the expression of the local initial-collinear counterterms in eq. (3.1). In general, in order for these counterterms to approximate the matrix element point by point in the phase space, spin correlations need to be taken into account. Suppressing the coupling orders, we define where s and s respectively specify the spin of the quark which enters the reduced amplitude and the corresponding conjugate. The factor ω(q)/ω(g) accounts for the different averaging on the initial state spins and colours in the matrix elements and equals N c /(N 2 c − 1) in four spacetime dimensions. The symbolP ss gq denotes the final-final qg splitting function and the variable z in eq. (A.1) is computed using For the process at hand, since the parton which enters the hard process after the splitting is always a quark, spin correlations are absent as indicated by δ ss in eq. (A.2). The momentum mapping that we use to determine the reduced phase-space point φ qH is the one used for two initial-state partons in Catani-Seymour dipole subtraction (see section 5.5 of [25]). 5 Finally, the Heaviside θ function at the end of eq. (A.1) controls the region of phase-space where the counterterm is active through the parameter y 0 , which determines the range for the variable y ≡ 2p 1 · p 3 /Q 2 . At this point, all the elements needed to check that eq. (3.1) only features integrable singularities have been presented. In order to validate our subtraction and assess the JHEP05(2019)002 numerical stability of the integrand which is built from the interference of two one-loop amplitudes, we start from a random resolved kinematic configuration and examine the behaviour as different collinear limits are approached. We control the distance from any given unresolved limit using a scaling variable λ, which is engineered to approach the singular configuration at a pace such that the phase-space volume between λ and λ + dλ is proportional to λ itself. For a kinematic configuration with a centre-of-mass energy of 1 TeV and in the case of a collinear pair, the typical invariant mass is then O(1 GeV) for λ = 10 −6 and O(1 MeV) for λ = 10 −12 . Under the same conditions and in the case of two collinear pairs, the typical invariant mass of each of them is O(1 GeV) for λ = 10 −12 and O(1 MeV) for λ = 10 −24 instead. For the sake of concreteness, we consider the partonic subprocess g 1 g 2 → b 3b4 H 5 , with the understanding that all qualitative features are identical in the case of light quark flavours in the final state.
In figure 4, we display the behaviour of the matrix element interference and its four initial-collinear counterterms as a function of λ for a given starting kinematic configuration. The left panels simply show the ratio of counterterms to the matrix element. In the right panels, we plot their sum weighted by λ, which is representative of the contribution to the total integral coming from a neighbourhood of λ. We therefore expect the integral to be convergent if this quantity tends to zero when λ → 0. In figure 4a we consider the limit C (1, 3), where the matrix element is approximated by the counterterm C(1, 3) and all other terms in the sum over π of eq. (3.1) are regular. The cases of C(2, 4), C(1, 4) and C(2, 3) are fully analogous. In figure 4b we study the limit of two collinear pairs C(1, 3)C(2, 4) which, as discussed in section 3, does not require any additional treatment since the matrix element for qq → H at order O g s g 2 is zero. Finally, in figure 4c we consider the limit C (3,4) to confirm that the matrix element of order O α 2 s α 2 for gg → bbH does not feature a non-integrable divergence when the two quarks in the final state are collinear. We note that the figures in this section can be sensitive to the numerical stability parameters of MadLoop which, among other things, control when to switch to a slower quadruple precision evaluation. Further discussion of this technical aspect is however beyond the scope of this work and we limit ourselves to reporting here that all Monte-Carlo integrations performed in this work could be successfully carried out using MadLoop's default parameters. 6 Incidentally, we observe that in order to obtain results at the level of precision needed for this work it is not necessary to introduce a technical cutoff.
The integral of the collinear counterterm over the unresolved phase space has been computed in [14] and reads 6 We note however that it proved to be necessary to employ an estimate of MadLoop's accuracy based on the comparison of two separate numerical evaluations that differ by a Lorentz transformation of the kinematic inputs (by setting MadLoop's parameter NRotations DP to 1).

JHEP05(2019)002
where Q = p 1 + p 2 for the reduced process g 1 q 2 → q 3 H 4 and we have defined The explicit pole in the dimensional regulator featured by this integrated counterterm is cancelled by the contribution from PDF renormalisation, which is given by where the relevant Altarelli-Parisi splitting kernel reads 7 We have confirmed that in our implementation of this subtraction scheme the sum of (3.1) and (3.2) does not depend on y 0 , which provides a non-trivial cross-check of O 0 terms.

B Analytical computation of the amplitudes A
where in the following we will denote the form factor by F µ 1 µ 2 , suppressing the spinor indices. We separate couplings of the quarks to the weak gauge bosons according to and refer to g V as the vector and to g A as the axial coupling constant. In the following, we restrict ourselves to the case of d quarks in the final state for concreteness. We will first discuss the computation of A (0,2),VV gg→ddH ∝ |g V | 2 and then argue that this piece is sufficient to determine the complete amplitude.
In order to compute the form factors, we first generate all contributing diagrams with QGraf [26] and perform the color-, Dirac-and Lorentz algebra in Mathematica, whereas 7 Note that in our notation the first subscript indicates the parton extracted from the hadron according to its PDF, and the second one denotes the parton that enters the hard process.

JHEP05(2019)002
the γ traces are performed using FORM [27]. For more compact expressions, we choose an axial gauge for the external gluons g 1 and g 2 , such that with the physical polarization sum As a direct consequence of the gauge choice, terms in the form factors proportional to p µ 1 1 , p µ 1 2 , p µ 2 1 or p µ 2 2 can be set to zero, since they will not contribute to the amplitude. Internal gauge bosons and quarks are treated in the Feynman gauge.

B.1 Tensor reduction to scalar integrals
The form factor F µ 1 µ 2 can be written as where the S µ 1 µ 2 i denote tensor integrals. We can reduce the tensor integrals S µ 1 µ 2 i to scalar integrals S k : To achieve the above decomposition, we use a particular flavour of Passarino-Veltman tensor reduction [28]. The reduction of the tensor integrals is discussed here by writing only the integrand numerators N µ 1 µ 2 (k), keeping in mind that the identities exclusively hold at the integral level.
As a first step we strip off the external Lorentz structures factorising the loop momentum N µ 1 µ 2 (k) and write where the tensor coefficients c i only involve γ-matrices, external momenta p i and the metric tensor g. The tensor reduction is performed with the fully symmetric tensor numerators N (α 1 α 2 ...αn) (k) = k α 1 k α 2 · · · k αn . Performing the loop integration of the tensor integral will result in Lorentz tensors t is chosen to match the convention of OneLOop [21], which is also used to evaluate the remaining non-trivial scalar master integrals.
In the evaluation of the expressions above, some care is needed in order to evaluate multi-valued functions on their physical Riemann sheet. The convention for numerical implementations of such functions is that the value assigned on the cut is the one coming around the finite endpoint of the cut in a counter-clockwise direction [33]. The Feynman prescription, however, dictates to replace s with s + iη and take the limit η ↓ 0, which gives if the right-hand side respects the convention. The same holds for s > 0 in the expansion of the massless bubble. It is easy to see that for the dilogarithm in (B.14), instead, the physical sheet coincides with the conventional one for all x = 1.

B.3 Relations between the axial and vector parts of the amplitude
In the previous section we discussed the computation of the vector part A (0,2),VV gg→ddH . In what follows, we restrict the discussion to a single quark family with a diagonal CKM matrix, The generalisation to all families of light quarks is straightforward and purely combinatorial. Since there are no closed fermion loops, we do not have to worry about ambiguous traces of γ 5 in 4 − 2 dimensions and we may take the D-dimensional γ 5 to be anticommuting. 9 JHEP05(2019)002 Figure 5. The relevant γ matrix structures for A (0,2),AA gg→ddH . Figure 5a corresponds to triangle diagrams, figure 5b corresponds to box diagrams and figure 5c corresponds to pentagon diagrams. Within the purely axial amplitude A (0,2),AA gg→ddH both weak couplings are ∝ g A γ µ γ 5 . The γ chains that appear in the amplitude are shown in figure 5. One always needs to do an even number of anticommutations to arrive at γ 5 γ 5 = 1 from which immediately follows that The axial-vector piece A (0,2),AV gg→ddH features the γ chains shown in figure 6. These chains represent the cases where only the vertex closest to the outgoing d quark (of momentum p 3 ) contributes with an axial coupling. 10 It is easy to see that an uneven number of anticommutations is needed to bring γ 5 to the beginning of every spinor chain appearing in the process. The form factor for the AV part of the amplitude is therefore given by We thus conclude that the complete weak amplitude can be determined from its purely vector piece.

B.4 Supplementary material
The notation employed for the supplementary material is the following: we write every form factor in the supplementary material as the scalar product

JHEP05(2019)002
The interference between two amplitudes A andÃ in this notation then reads where B is the structure matrix obtained summing over all colours, spins and polarisations: The supplementary material contains the vector T µν,ab s 1 s 2 ,lm and the vector S for the QCD background, the VV and the AV part of the weak amplitude. We furthermore provide the structure matrices B for A A * (2,0) gg→ddH , which are sufficient to reproduce analytically the one-loop mixed QCD-weak matrix element for light quarks (excluding Higgs-strahlung contributions).

C Validation material
In order to facilitate the reproduction of our results, we provide below the numerical result for the matrix element M gg→Hbb summed (averaged) over final (initial) state helicity and colour configurations for the following two kinematic points and α s = 0.118 (other SM parameters set to the values indicated in table 2, unless otherwise stated).
The matrix elements computed are free of any explicit IR or UV divergence, so that the specific -dependent normalisation factor considered in MadLoop's conventions is irrelevant in this case. For the two kinematic points shown in table 4, we find: The first two matrix element evaluations given in table 5 are exactly those used for obtaining the results of table 3. The next six correspond to simplified setups that are only meant to ease comparisons against independent computations. More specifically, the matrix element denoted M (α 3 s α 2 ,Γ t,W ± ,Z =0,W ± @[1a, 1b, 1c],V V ) gg→Hdd corresponds to the case where: • all widths are set to zero (then using on-shell renormalisation conditions) • only the diagrams from the classes 1a, 1b and 1c with a W ± in the loop are kept • only the vector part of the two W ± interactions is considered.
The definition of the last five matrix elements of the table is fully analogous, with 'AV+VA' indicating that the amplitude includes exactly one vector-like and one axial coupling of the electroweak boson to the quarks.
For each matrix element we checked numerical evaluations of the analytic result for 100 phase-space points and compare them against MadLoop evaluations. We found perfect agreement at the level of the 10 th digit on average.
The above matrix elements can readily be generated by MadLoop (from within MG5aMC v2.6+) using commands similar 11 to the following which generates the matrix element M Evaluation of analytic result 5.365688093207525e-12 7.085433478511895e-11 Table 5. Benchmark evaluations of various matrix elements comparing numerical results from MadLoop against an independent analytical derivation of the amplitude, presented in appendix B (see text for details).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.