One-loop weak corrections to Higgs production
Valentin Hirschi Simone Lionetti Armin Schweitzer


EW

[
    Endre Süli
    December 30, 2023
=====================

fancy



We compute mixed QCD-weak corrections to inclusive Higgs production
at the LHC from the partonic process g g → H q q̅.
We start from the UV- and IR-finite one-loop weak amplitude
and consider its interference
with the corresponding one-loop QCD amplitude.
This contribution is a
𝒪(α_sα)
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 𝒪(α^2).
Finally, we consider contributions
from the partonic process g q → H q,
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 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 <cit.>.
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% <cit.>.

On the other hand, EW corrections to the LO inclusive Higgs cross section also need to be considered.
In the same setup mentioned before, the first EW term turns out to increase the total gluon fusion cross section by a significant 5% <cit.>.
Since NLO QCD corrections can be as large as the leading contribution,
the motivation to investigate mixed first-order QCD and first-order EW 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 argument
that mixed QCD-EW effects on the inclusive Higgs production cross section
are well approximated by combining the purely EW term
and the full QCD series in a multiplicative fashion <cit.>.
Following this factorisation approach, the authors of ref. <cit.> reported the mixed QCD-EW 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 <cit.> 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 𝒪(α_S^2α)  cross section, which we verified by explicit calculation.
The theoretical uncertainty associated to each of the other main error sources
(determination of PDF, 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-EW correction
to Higgs boson gluon fusion for arbitrary masses
of the W, Z, and Higgs bosons <cit.>,
an estimate of the cross section in the soft-virtual approximation
was obtained <cit.>.
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 <cit.>.
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-EW 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 EW 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 <cit.>,
the contributions with two extra hard partons are formally of the same order
and may disrupt the approximate factorisation of EW 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-EW corrections from the one-loop partonic subprocess g g → H q q̅.
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 <ref>,
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 some care
and are thus described in Section <ref>.
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 EW expansion as

    σ_pp→ H+X = ∑_m,nσ^(m,n)_pp→ H+X,
    	    where    σ^(m,n)_pp→ H+X∝α_s^m+2α^n+1,

where the prefactor α_s^2α^1 is chosen so as to match the couplings factorised by the leading-order 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^mLO” and “(electro)weak N^nLO” 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:

    A^(i,j)_a b→ H+X∝ g_s^i+2 g^j+1,

where we denote by g all couplings that are not g_s.

As mentioned above,
EW and QCD corrections are expected to factorise to a certain degree,
such that

    σ^(m,n)_pp→ H+X∼σ^(m,0)_pp→ H+X·σ^(0,n)_pp→ H+X.

This approximation is valid under the assumption that the main contributions
to the mixed QCD-EW cross section are to be attributed
either to soft gluons or Sudakov EW 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 PDF 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 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.
[Note that our initial-state notation g q encompasses in this context both permutations g q and q g.]

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 σ^(0,2)_gg→ H q q̅ and σ^(-1,2)_gq→ H q
built from the square of the amplitudes A^(0,2)_gg→ H q q̅ and A^(-1,2)_gq→ H q.
These form a gauge-invariant subset of higher-order contributions.

Our work reports on the contribution σ^(1,1)_gg→ Hq q̅ for the first time and, together with the results from refs. <cit.>, it completes the computation of σ^(1,1)_g g → H + X.
We now proceed to list in Table <ref> all amplitudes building σ^(1,1)_gp→ H+X.




We now turn to discussing the Feynman diagrams building the amplitudes A^(2,0)_gg→ H q q̅, A^(0,2)_gg→ H q q̅, A^(1,0)_gq→ H q and A^(-1,2)_gq→ H q that contribute to the cross sections presented in this work.

The amplitude A^(2,0)_gg→ H q q̅ is built from the diagrams depicted in Fig. <ref> where the Higgs is produced via weak vector boson fusion and interfered with the leading QCD gluon-fusion diagrams shown in Fig. <ref>.

Diagrams of the class <ref> and <ref>, where the Higgs is produced via gluon-fusion, feature a Z-boson propagator[The diagram analogous to <ref> with a photon instead of the Z-boson is exactly zero in virtue of Furry's theorem.] which however does not yield any Breit-Wigner resonance as they are interfered against the QCD diagrams of Fig. <ref>.
We must nonetheless regulate the Z-boson propagator pole, which motivates our use in this computation of the complex-mass scheme <cit.> with finite widths for the internal top quark and unstable weak gauge bosons.
These diagrams <ref> and <ref> 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 narrow-width 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 <ref> 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. <cit.> and, for this reason, we found it interesting to report our results separately for the processes g g → H q q̅, with q ≡ u,d,c,s, and g g → b b̅ H.










§ 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→ H q q̅ as it constitutes the main focus of the present work.

In principle, the process g g → H q q̅ presents infrared divergences
when the quark-antiquark pair in the final state is collectively soft,
and/or when one or both of the quarks are collinear
to the direction of an incoming gluon.
However, thanks to the factorisation properties of QCD,
in double-unresolved configurations the amplitude A^(0,2)_gg→ H q q̅
can be approximated by universal factors times the reduced amplitude  A^(-2,2)_qq̅→ H (that is, of order g^3) which is identically zero.
Indeed, the triangle one-loop diagrams for q q̅→ 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 Fig. <ref>.
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 <ref> and <ref>.
In the case of the diagram <ref>, this happens
when antiquark d̅_5 becomes collinear to gluon g_2
such that the hard scattering subgraph corresponds to diagram <ref>.
By contrast, in the kinematic limit where
quark d_4 is collinear to gluon g_1
and quark d̅_5 is collinear to gluon g_2,
both non-loop propagators of graph <ref> 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 <ref> matches that of graph <ref>.

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 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:

    σ^(m,n)_g g → H q q̅ =    ∫_0^1 x_1∫_0^1 x_2 f_g(x_1) f_g(x_2) {
           ∫Φ_H q q̅[
    ℳ^(m,n)_g g → H q q̅𝐉(ϕ_H qq̅)
    - ∑_π𝐂_g q⊗ℳ^(m-1,n)_g q → H q𝐉(ϕ̃_H q )
    ]
    
           
    + ∑_π∫Φ_H q∫_0^1 ξ[
    	𝐂_gq(ξ)
    	+ Δ_g q(ξ)
    ]
    ℳ^(m-1,n)_g q → H q𝐉(ϕ_H q)
    },

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 𝐂_ij denotes the local counterterm
for particles i and j going collinear
and 𝐂_ij its counterpart analytically integrated over the unresolved degrees of freedom.
The observable functions are indicated with 𝐉,
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. <cit.>
and are presented more explicitly in appendix <ref>,
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^(0,2)_gg→ H q q̅ and A^(-1,2)_gq→ H q
that factorise a Higgs coupling to weak bosons 
were first computed analytically (for massless quarks only) in ref. <cit.>,
in the different context of NLO QCD corrections to weak vector-boson fusion.
In the present case and as indicated in Table <ref>, in order to obtain contributions
to σ^(1,1)_gg→ H q q̅ and σ^(0,1)_gq→ H q,
these amplitudes must be interfered against their corresponding QCD analog.

Nowadays such one-loop amplitudes are readily available 
from many automated one-loop 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  <cit.>, 
part of  <cit.> (henceforth abbreviated ), as it can efficiently generate and
interfere <cit.> arbitrary one-loop amplitudes in the SM and beyond.
 uses  <cit.> and  <cit.>, or alternatively  <cit.>, 
for performing one-loop reductions and for the evaluation of the scalar one-loop master integrals.
We present in appendix <ref> 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 's numerical implementation of the amplitudes A^(2,0)_gg→ H q q̅ and A^(0,2)_gg→ H q q̅
against a completely independent and analytical computation described in appendix <ref>.

As already mentioned, we choose to renormalise all unstable particles in the complex-mass scheme <cit.> and consider the SM input parameters given in Table <ref>.



The numerical Monte-Carlo integration as well as the necessary IR subtraction procedure, presented in Eqs. <ref> and <ref> as well as in appendix <ref>, have been implemented in a private extension of  currently under development.
The poles in the dimensional regulator ϵ have been checked to cancel as expected.
[This check of course only considers the convoluted term of Eq. <ref> as our computation involves no virtual contribution. Also, for the pole cancellation to occur, it is important to restrict the initial state contributions to gluons only, as poles from the beam factorisation terms Δ_q g and Δ_q q remain uncanceled given that we ignore the corresponding real-emission subprocesses.]
Moreover, we have validated our code by comparing NLO QCD cross sections
against results from  for the processes p p → Z and p p → H,
the latter in the Higgs Effective Theory.




Our results are presented in Table <ref>.
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 α_s^2α^3 can be suppressed compared to their α_s^3α^2 counterpart by less than what is expected by their parametric ratio α/α_s.
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 α_s^mα^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 quark-initiated EW corrections are however still small in comparison with the whole σ^(α_s^3α^2)_pp→ H+X, and can thus be safely neglected as already observed in  ref. <cit.>.
These two observations reinforce the conclusion that the contributions to σ^(α_s^3α^2)_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 σ^(α_s^3α^2)_gg→ Hbb̅ from only final-state b quarks reveals that contributions featuring Higgs production from the internal top quark line (see Fig. <ref>) are comparable and of opposite sign to that of emissions from internal weak bosons.
This fact contrasts with the study of ref. <cit.> 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 three-loop amplitude A^(2,2)_gg→ H of refs. <cit.>.
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 σ^(α_s^2α^2)_bg→ Hb at the same level as that of the channels initiated by each other valence quark flavour.

The squared contribution σ^(α_s^2α^3,no-HZ)_gg→ Hbb̅ omits the diagrams <ref> and <ref> 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 σ^(α_s^2α^2,Γ̅_Z=0)_gg→ HZ only to serve as an upper bound to this contribution.

Notice that the contribution σ^(α_s^2α^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 𝐂_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[This is for example the case in σ^(α_s^2α^3,no-ZH)_gg→ Hbb̅|_μ_F=m_H where Higgs emissions from weak bosons are close to equal and opposite in sign to emissions from internal top-quarks, and in σ^(α_s^2α^3,no-ZH)_gg→ Hqq̅|_μ_F=m_H/2 where the cancellation occurs between the hard reals and the logarithms in ξ part of the integrated counterterms 𝐂_gq. ]
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 mixed QCD-EW corrections, which is estimated to be of the order of 2 pb with an associated uncertainty in the range of 200 fb <cit.>.
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 EW 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 
σ^(α_s^2α^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. <cit.> 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 EW 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 EW correction of order 𝒪(α_s^3α^2) to the Higgs inclusive cross section.

To this end, two groups <cit.> computed σ^(1,1)_gg→ H+X and found that it supports the hypothesis that EW corrections factorise. These works however neglected the quark-initiated components as well as the “double-real” channel g g → H q q̅ and we confirm here that these terms can be safely neglected, amounting to about 5% of the total mixed EW 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→ H qq̅ and σ^(1,1)_gq→ H q renders them prone to kinematic suppressions, and we indeed found that the square of the one-loop EW 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.



§ ACKNOWLEDGMENTS


We are greatly indebted to
C. Anastasiou,
V. del Duca,
N. Deutschmann,
E. Furlan,
B. Mistlberger
and C. Specchia
for fruitful discussions.
This project has received funding
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.


equationsection




§ INITIAL-COLLINEAR COUNTERTERMS



In this appendix, we detail all ingredients that are necessary
for the subtraction of implicit singularities
outlined in Eqs. <ref>–<ref>,
and we demonstrate that the matrix element for g_1 g_2 → q_3 q̅_4 H_5
is correctly regulated in the regions of phase space
close to unresolved configurations.
As already announced in Section <ref>,
the construction and notation follow closely ref. <cit.>.

Let us begin with the expression of the local initial-collinear counterterms
in Eq. <ref>.
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

    [𝐂_g_1q_3⊗ℳ_q̅_13 g_2 →q̅_4 H_5]
    	(ϕ_q_3 q̅_4 H_5)
    	≡ (8πα_sμ^2ϵ)
    	1/s_g_1 q_3P̂_g_1q_3^ss'(1/z)
    	ω(q̅_13)/ω(g_1)
    ×ℳ_q̅_13 g_2 →q̅_4 H_5^ss'
    	(ϕ̃_q̅_4 H_5)
    	θ(y_0-y),

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_c^2-1) in four spacetime dimensions.
The symbol P̂_gq^ss' denotes the final-final qg splitting function

    P̂_gq^ss'(z)
    	≡δ^ss' C_F [ 1+(1-z)^2/z - ϵ z ],

and the variable z in Eq. <ref> is computed using

    z ≡Q· (p_1-p_3)/Q· p_1,

with Q = p_1+p_2.
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. <ref>.
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 <cit.>).
[
Note that this mapping involves recoiling against all final state particles,
and it would not be efficient for studying differential Higgs observables.
]
Finally, the Heaviside θ function at the end of Eq. <ref>
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. (<ref>)
only features integrable singularities have been presented.
In order to validate our subtraction
and assess the 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 λ+λ
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 1 GeV for λ=10^-6
and 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 1 GeV for λ=10^-12
and 1 MeV for λ=10^-24 instead.
For the sake of concreteness,
we consider the partonic subprocess g_1 g_2 → b_3 b̅_4 H_5,
with the understanding that all qualitative features are identical
in the case of light quark flavours in the final state.



In Fig. <ref>, 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 Fig. <ref> 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. (<ref>) are regular.
The cases of C(2,4), C(1,4) and C(2,3) are fully analogous.
In Fig. <ref>
we study the limit of two collinear pairs C(1,3)C(2,4) which,
as discussed in Section <ref>,
does not require any additional treatment
since the matrix element for q q̅→ H
at order g_s g^2 is zero.
Finally, in Fig. <ref> we consider the limit C(3,4)
to confirm that the matrix element of order α_s^2α^2
for g g → b b̅ H 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  
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 's default parameters[We note however that it proved to be necessary to employ an estimate of 's accuracy based on the comparison of two separate numerical evaluations that differ by a Lorentz transformation of the kinematic inputs (by setting 's parameter NRotations_DP to 1).].
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 <cit.> and reads

    𝐂_gq(ξ) =
    	α_s/2π S_ϵ[μ_R^2/Q^2]^ϵT_R/C_F{
    		[ξ^2+(1-ξ)^2] [ -1/ϵ
    		+ ln(1-ξ)  (1+θ[ξ-(1-y_0)])
    		
    
    		+ ln(y_0)θ[(1-y_0)-ξ]
    		]
    		+2ξ(1-ξ)
    	}
    	+ϵ,

where Q = p_1+p_2 for the reduced process g_1 q_2 → q_3 H_4
and we have defined

    S_ϵ≡(4π)^ϵ/Γ(1-ϵ).

The explicit pole in the dimensional regulator ϵ 
featured by this integrated counterterm
is cancelled by the contribution from PDF renormalisation,
which is given by

    Δ_gq(ξ)
    	=
    	α_s/2π S_ϵ1/ϵ[μ_R^2/μ_F^2]^ϵ
    	P_gq(ξ),

where the relevant Altarelli–Parisi splitting kernel reads
[
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.
]

    P_gq(ξ) = T_R [ξ^2+(1-ξ)^2].

We have confirmed that in our implementation of this subtraction scheme
the sum of (<ref>) and (<ref>) does not depend on y_0,
which provides a non-trivial cross-check of ϵ^0 terms.




§ ANALYTICAL COMPUTATION OF THE AMPLITUDES A^(2,0)_GG→ H Q Q̅ AND A^(0,2)_GG→ H Q Q̅



The analytic validation is performed by computing the form factors for the QCD background depicted in the diagrams Fig. <ref> and the EW contributions with sample diagrams shown in Figs. <ref> to <ref>.
While we retain the full quark mass dependence for the QCD background,
we assume massless quarks for the EW contributions.
The computation is performed in D=4-2ϵ dimensions.
However, due to the special reduction of the scalar pentagon integrals, the final result is only valid for the provided order ϵ^0 in the dimensional regulator (see Sec. <ref>).

The amplitude for the EW process g_1 g_2→ q_3 q̅_4 H_5 may be written in the general form

    A^(0,2)_g g → q q̅ H   =A^(0,2),_g g → q q̅ H+A^(0,2),_g g → q q̅ H+A^(0,2),_g g → q q̅ H+A^(0,2),_g g → q q̅ H
       =ε_μ_1(p_1)ε_μ_2(p_2)u̅^s_3(p_3)v^s_4(p_4)ℱ_s_3 s_4^μ_1 μ_2,

where in the following we will denote the form factor by ℱ^μ_1 μ_2, suppressing the spinor indices.
We separate couplings of the quarks to the EW gauge bosons according to

    u d W^+∝ g_γ^μ + g_γ^μγ^5,     
    	u d W^-∝ g_^* γ^μ + g_^* γ^μγ^5   and
        q q Z∝ g_,Zγ^μ+g_,Zγ^μγ^5,

and refer to g_ as the vector and to g_ 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),_g g → d d̅ H∝ |g_|^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  <cit.> and perform the color-, Dirac- and Lorentz algebra in Mathematica, whereas the γ traces are performed using  <cit.>.
 For more compact expressions, we choose an axial gauge for the external gluons g_1 and g_2, such that

    p_1·ε_1(p_1)=0        p_2·ε_1(p_1)=0        p_2·ε_2(p_2)=0        p_1 ·ε_1(p_2)=0,

with the physical polarization sum

    ∑_polarizationε^*_μ(p_i)ε^_ν(p_i)
    	=-g_μν+p_1μ p_2ν+p_2μ p_1ν/p_1 · p_2    for  i=1,2.

As a direct consequence of the gauge choice, terms in the form factors proportional to p_1^μ_1, p_2^μ_1, p_1^μ_2 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. 




 §.§ Tensor reduction to scalar integrals


The form factor ℱ^μ_1μ_2 can be written as

    ℱ^μ_1μ_2
    	=∑_i α_i 𝒮^μ_1μ_2_i
    	=∑_i α_i ∫N^μ_1 μ_2(k)/D_1… D_m_i[D]k,

where the 𝒮^μ_1μ_2_i denote tensor integrals.
We can reduce the tensor integrals 𝒮^μ_1μ_2_i to scalar integrals S_k:

    𝒮^μ_1μ_2_i=∑_j T_j^μ_1μ_2(p_1,p_2,p_3,p_4)S_j.

To achieve the above decomposition, we use a particular flavour of Passarino–Veltman tensor reduction <cit.>.
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
 
    N^μ_1μ_2(k)=
     	c_0^μ_1μ_2
     	+c_1^μ_1 μ_2; α_1 k_α_1
     	+c_2^μ_1 μ_2; α_1 α_2 k_α_1k_α_2
     	+c_3^μ_1 μ_2; α_1 α_2 α_3 k_α_1k_α_2k_α_3+ …,

 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 Ñ^(α_1α_2…α_n)(k)=k^α_1k^α_2⋯ k^α_n. Performing the loop integration of the tensor integral

    𝒮^(α_1α_2…α_n)_i
    =∫Ñ^(α_1α_2…α_n)(k)/D_1… D_m_i[D]k
    =∑_j t_j^(α_1α_2…α_n)(g,p_1,p_2,p_3,p_4)c_j,

will result in Lorentz tensors t_j^(α_1α_2…α_n) which are also completely symmetric in the internal Lorentz indices α_i. The symmetric tensor t^(α_1…α_n) is given by

    t^(α_1…α_n)=1/n!∑_σ∈Σ_nt^α_σ(1)…α_σ(n),

where Σ_n is the symmetric group of order n, e.g.

    t_p_1,p_2^(α_1α_2)=p_1^(α_1 p_2^α_2) 
    	=1/2(p_1^α_1p_2^α_2+p_1^α_2p_2^α_1).

Reduction with respect to a fully symmetric tensor basis reduces the number of tensor structures to be considered in the Ansatz significantly, while still remaining completely algorithmic.

A further simplification arises from the fact that the only underlying scalar topology depending on the full external kinematics p_1,p_2,p_3,p_4 is the pentagon displayed in Fig. <ref>.
Every other diagram will yield scalar integrals with reducible external kinematics yielding results depending on a reduced set of Mandelstam variables only.
In general we group all diagrams into families characterized by their minimal set of external momenta and perform the tensor reduction for the tensor numerators Ñ^(α_1α_2…α_n)(k) in each family separately with respect to the reduced external kinematics.
This approach keeps the intermediate expressions obtained from tensor- and IBP reduction very compact.
The biggest matrix we have to invert is a 24×24 matrix for the rank 3 Lorentz tensor integrals of the pentagon diagram.
To perform the analytic matrix inversion we employ the computer algebra system  <cit.>, which takes below a minute on one core of a modern computer.
In order to obtain the form factor in terms of scalar integrals we then insert the solutions back into (<ref>).

With the tensor decomposition described, we are able to express the form factors for the QCD background and the vector-vector part of the EW contributions in terms of the following 20 tensor structures:

    T_1^μ_1μ_2   =p_1^_γ^μ_1γ^μ_2,   
    	T_2^μ_1μ_2   =p_2^_γ^μ_1γ^μ_2,    
    	T_3^μ_1μ_2   =g^μ_1 μ_2p_1^_,    
    	T_4^μ_1μ_2   =g^μ_1 μ_2p_2^_,
     
    	T_5^μ_1μ_2   =γ^μ_1 p_3^μ_2,    
    	T_6^μ_1μ_2   =γ^μ_1 p_4^μ_2,    
    	T_7^μ_1μ_2   =p_1^_p_2^_γ^μ_1 p_3^μ_2,    
    	T_8^μ_1μ_2   =p_1^_p_2^_γ^μ_1 p_4^μ_2,
     
    	T_9^μ_1μ_2   =γ^μ_2 p_3^μ_1,    
    	T_10^μ_1μ_2   =p_1^_p_2^_γ^μ_2 p_3^μ_1,    
    	T_11^μ_1μ_2   =γ^μ_2 p_4^μ_1,    
    	T_12^μ_1μ_2   =p_1^_p_2^_γ^μ_2 p_4^μ_1,
     
    	T_13^μ_1μ_2   =p_1^_ p_3^μ_1 p_3^μ_2,    
    	T_14^μ_1μ_2   =p_2^_ p_3^μ_1 p_3^μ_2,    
    	T_15^μ_1μ_2   =p_1^_ p_3^μ_1 p_4^μ_2,    
    	T_16^μ_1μ_2   =p_2^_ p_3^μ_1 p_4^μ_2,
     
    	T_17^μ_1μ_2   =p_1^_ p_3^μ_2 p_4^μ_1,    
    	T_18^μ_1μ_2   =p_2^_ p_3^μ_2 p_4^μ_1,    
    	T_19^μ_1μ_2   =p_1^_ p_4^μ_1 p_4^μ_2,   
    	T_20^μ_1μ_2   =p_2^_ p_4^μ_1 p_4^μ_2.





 §.§ Evaluation of scalar integrals



The IBP reduction of the remaining scalar integrals is performed using the program  <cit.>.
We decompose the scalar pentagon integrals appearing as master integrals following ref. <cit.>.
This decomposition relates the pentagon in 4-2ϵ dimensions to a linear combination of all boxes obtainable by pinching one of the propagators, and the pentagon in 6-2ϵ dimensions multiplied by a prefactor of order ϵ.
Since the pentagon in six dimensions is finite, the additional term involving the six-dimensional pentagon is of order ϵ and can be omitted for the computation at hand.

We find that all form factors are finite diagram-by-diagram, but order ϵ coefficients of the bubbles appear explicitly in the final amplitude.
[This is particular to our approach and originates from the IBP reduction of the scalar integrals.]
The relevant coefficients are given by

    b_0(s;m^2,0)=h(ϵ)∫1/(k^2-m^2)(k-p)^2[D]k
        =
    	(μ^2/m^2)^ϵ[ 1/ϵ
    		+2-(x-1) ln(1-x)/x.
        .+(4+π^2/6+(1-x)/2 x(2 _2(-x/1-x)-ln^2(1-x)+4 ln (1-x)))ϵ+ϵ^2],
    	
       b_0(s;0,0)=h(ϵ)∫1/k^2(k-p)^2[D]k
    	=(-s/μ ^2)^-ϵ[1/ϵ+2+4 ϵ +ϵ^2], 
       a_0(m)=h(ϵ)∫1/k^2-m^2[D]k
    		=	m^2 (μ ^2/m^2)^ϵ[ 1/ϵ+1+(1+π ^2/6) ϵ + ϵ^2],

where x=s/m^2.
The normalization

    h(ϵ)
    	≡μ^2ϵ/i π^2-ϵΓ(1-2ϵ)/Γ(1-ϵ)^2 Γ(ϵ+1),

is chosen to match the convention of  <cit.>,
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
<cit.>. 
The Feynman prescription, however, dictates to replace s with s+iη
and take the limit η↓ 0, which gives

    lim_η↓ 0ln(1-s+iη/m^2)=
    	ln(1-x)    s<m^2, 
    ln(1-x)-2iπ    s>m^2,

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 (<ref>), instead,
the physical sheet coincides with the conventional one for all x≠1.




 §.§ 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),_g g → d d̅ H.
In what follows, we restrict the discussion to a single quark family with a diagonal CKM matrix, g_=g_^* and g_=g_^*.
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.
[
Our choice corresponds to the NDR treatment of γ^5 (see e.g. <cit.>).
Note however, that the amplitude is finite diagram by diagram and traces over γ^5 enter only in the interference of the AV-part with e.g. the QCD background.
Since for the interference there is no explicit ϵ-dependence anymore, the traces can be treated as four-dimensional objects, without the need of imposing additional constraints.
]





Within the purely axial amplitude A^(0,2),_g g → d d̅ H both EW couplings are ∝ g_γ^μγ^5.
The γ chains that appear in the amplitude are shown in Fig. <ref>.
One always needs to do an even number of anticommutations to arrive at γ^5γ^5=1 from which immediately follows that

    A^(0,2),_g g → d d̅ H
    	=|g_|^2/|g_|^2A^(0,2),_g g → d d̅ H.
 
The axial-vector piece A^(0,2),_g g → d d̅ H features the γ chains shown in Fig. <ref>.
These chains represent the cases where only the vertex closest to the outgoing d quark (of momentum p_3) contributes with an axial coupling.
[The case where only the vertex “furthest” to the outgoing d quark contributes with the axial coupling is completely analogous.]
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  part of the amplitude is therefore given by

    ℱ^μν_ =
    	-2 g_/g_γ^5 ℱ^μν_.

We thus conclude that the complete EW amplitude can be determined from its purely vector piece.




 §.§ Ancillary files


The notation employed for the ancillary files is the following:
We write every form factor in the ancillary files as the scalar product

    ℱ_s_1 s_2,lm^μν,ab = (T_s_1 s_2,lm^μν,ab)_i S_i,

where the vector T spans the direct product of colour and Lorentz spaces.
The pairs of indices a, b and l, m are associated with the adjoint and fundamental representations of SU(3) respectively;
μ, ν are the Lorentz indices and s_1, s_2 are the spinor ones.

The interference between two amplitudes A and Ã in this notation then reads

    ℳ = 21/41/(N_c^2-1)^2Re(s̃^*_i B_ij s_j),

where B is the structure matrix obtained
summing over all colours, spins and polarisations:

    B_ij = ∑ε_μ(p_1)ε_ν(p_2)
    ε^*_μ'(p_1)ε^*_ν'(p_2)
    u^s_1'(p_3)v̅^s_2'(p_4)u̅^s_1(p_3)v^s_2(p_4)
    ×(T̃_s_1' s_2'^μ' ν')_i(T_s_1 s_2^μν)_j.

The ancillary files contain the vector T_s_1 s_2,lm^μν,ab and the vector S for the QCD background, the  and the  part of the EW amplitude.
We furthermore provide the structure matrices B for
A^(0,2),,()_g g → d d̅ H𝒜^*^(2,0)_g g → d d̅ H
and
A^(0,2),,()_g g → d d̅ H𝒜^*^(2,0)_g g → d d̅ H,
which are sufficient to reproduce analytically the one-loop mixed QCD-EW matrix element for light quarks (excluding Higgs-strahlung contributions).





§ VALIDATION MATERIAL



In order to facilitate the reproduction of our results, we provide below the numerical result for the matrix element ℳ^(α_s^3α^2)_gg→ Hdd̅ and ℳ^(α_s^3α^2)_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 <ref>, 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 's conventions is irrelevant in this case. For the two kinematic points shown in Table <ref>, we find:




The first two matrix element evaluations given in Table <ref> are exactly those used for obtaining the results of Table <ref>. The next six correspond to simplified setups that are only meant to ease comparisons against independent computations.
More specifically, the matrix element denoted ℳ^(α_s^3α^2,Γ_t,W^±,Z=0,W^±@[<ref>,<ref>,<ref>],VV)_g g → H d d̅ corresponds to the case where:


  * all widths are set to zero (then using on-shell renormalisation conditions)

  * only the diagrams from the classes <ref>, <ref> and <ref> 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  evaluations. We found perfect agreement at the level of the 10^ th digit on average.

The above matrix elements can readily be generated by  (from within  v2.6+) using commands similar[See <https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/MadLoopStandaloneLibrary> for instructions on how to generate the corresponding standalone library for linking against your own code.] to the following which generates the matrix element ℳ^(α_s^3α^2)_gg→ Hdd̅:



   set complex_mass_scheme True

   import loop_qcd_qed_sm

   generate g g > h d d [virt=QCD QED] QED2==4 QCD2==6

   output my_gg_hddx

   launch -f


Note that in order to select only the diagrams of the classes <ref>, <ref> and <ref>,
the following –loop_filter option
[
Also note that the rather long loop filters indicated on the command line can alternatively be specified directly in the user-function user_filter() of the   module loop_diagram_generation.py.
]
can be passed to the following generate command, yielding the matrix element ℳ^(α_s^3α^2,Γ_t,W^±,Z=0,Z@[<ref>,<ref>,<ref>],+++)_g g → H d d̅:



 

–loop_filter=not(23\\ in\\ struct_pdgs\\ or\\ 250\\ in\\ struct_pdgs) 






JHEP


