Next-to-leading-order QCD corrections to Higgs production in association with a jet

We compute the next-to-leading-order (NLO) QCD corrections to the Higgs p T distribution in Higgs production in association with a jet via gluon fusion at the LHC, with exact dependence on the mass of the quark circulating in the heavy-quark loops. The NLO corrections are presented including the top-quark mass, and for the ﬁrst time, the bottom-quark mass as well. Further, besides the on-shell mass scheme, we consider for the ﬁrst time a running mass renormalisation scheme. The computation is based on amplitudes which are valid for arbitrary heavy-quark masses. © 2023 The Author(s). Published by Elsevier B


INTRODUCTION
Since the discovery of the Higgs boson [1,2], one of the main goals of the Large Hadron Collider (LHC) physics program has been to investigate couplings and quantum numbers of the Higgs boson as accurately as possible.At the LHC, the dominant Higgs production mode is via gluon fusion, with the coupling of the Higgs boson to gluons being mediated by a heavy-quark loop.This gives the opportunity to test the Standard Model (SM) and to look for possible deviations from it, which would be footprints of New Physics (NP) [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19].
A promising observable to probe possible NP effects is the Higgs p T distribution [20][21][22][23][24], which allows one to analyse how the production rate depends on the heavyquark loop at different p T values.The Higgs p T distribution is known at leading order in α s for arbitrary quark masses in the heavy-quark loop [25][26][27], and at NLO for the top-quark mass [28,29] in the on-shell (OS) mass scheme.Approximate results are known at NLO [30], which include the bottom-quark mass [31][32][33] as well as the top-bottom interference [34], and beyond NLO in the high-energy limit [35].Further, the mixed QCDelectroweak contributions to the Higgs p T distribution are known at leading order [36].The Higgs Effective Field Theory (HEFT) approach provides a good approximation to the Higgs p T distribution when p T < m t [26].In HEFT, the Higgs p T distribution is known at NNLO in α s [37][38][39][40].Conversely, when p T > m t the shape of the Higgs p T distribution in HEFT overshoots the exact Higgs p T distribution [29].
In this letter, we present the Higgs p T distribution at NLO in α s , with top-quark mass dependence, and for the first time including also the exact bottom-quark mass contribution.Further, since the top-and bottom-quark masses are treated as dynamical parameters and not as fixed numbers, besides the OS mass scheme, we compute the Higgs p T distribution using for the first time a dynamical mass renormalisation scheme, the MS scheme.

CALCULATION
Our computation of the Higgs p T distribution is based on amplitudes which are valid for arbitrary quark masses circulating in the heavy-quark loops.
We consider the production of a Higgs boson in association with a jet in proton-proton collisions, p + p → H + j + X.The corresponding cross section can be found by convoluting the partonic cross section with the respective parton distribution functions (PDF), for the partonic channels: gg → Hg, q q → Hg and q(q)g → Hq(q).At LO in the coupling constants, this involves the one-loop 2 → 2 amplitudes.The coupling of the Higgs boson to light quarks is suppressed by the light-quark mass.Therefore, only the contributions in which the Higgs boson is coupled to a heavy quark are important.In this letter we consider the contribution of the top and the bottom quarks.
At NLO in QCD, we have to consider the O(α S ) corrections to the LO 2 → 2 amplitudes, together with the one-loop 2 → 3 amplitudes, in which the Higgs boson and the jet are produced with an additional parton in the final state: gg → Hgg, q q → Hgg, and q(q)g → Hq(q)g.The calculation of the two-loop 2 → 2 amplitudes, which has been performed assembling the relevant master integrals [41][42][43], is briefly described in the following subsection and will be the subject of a forthcoming publication [44].The calculation of the one-loop 2 → 3 amplitudes, which are known analytically [45,46] and numerically through public automated tools, has been performed using the results of refs.[46,47] reported in the program MCFM-9.1 [48].We have cross-checked the one-loop 2 → 3 amplitudes with automated tools (MG5 aMC@NLO [27,49] and GoSam [50]) and found perfect agreement for every phase-space point we compared, including points where an unresolved limit is approached.
After renormalisation of the ultra-violet (UV) divergences, each of the two sets of amplitudes still contains infra-red (IR) divergences.In order to regulate such divergences, we used the dipole subtraction scheme [51].For the numerical evaluation of the 2 → 3 amplitudes and the subtraction terms we have used the program MCFM-9.1 [48].Although the integration of the 2 → 3 amplitudes, even when generated with the automated tools, is not time intensive, we preferred to use the analytic result of refs.[46,47] which saved about a factor of a hundred in the integration time of the gluonic one-loop 2 → 3 amplitudes.

TWO-LOOP 2 → 2 AMPLITUDES
The two-loop 2 → 2 amplitudes for the processes under consideration can be written in terms of form factors, exploiting their Lorentz and Dirac structures, where every form factor can be expressed as a combination of scalar integrals that are divergent in four space-time dimensions.We regularise both UV and IR divergences in dimensional regularisation.Since the scalar integrals are not all independent, we use Integration-by-Parts Identities [52][53][54], implemented in the computer programs FIRE [55,56] and Kira [57,58], in order to reduce them to a set of independent integrals, called Master Integrals (MIs).We calculate the MIs using the differential equations method [59][60][61][62][63][64][65][66].The MIs are expressed as Laurent series in the dimensional parameter = (4−d)/2, where d is the dimension of space-time.The system of first-order linear differential equations satisfied by the MIs is solved at every order in using DiffExp [67], a Mathematica implementation of the series-expansion method [68].The evaluation of all the 447 MIs involved in the calculation is described in detail in [41][42][43].For a number of points in the physical phase space region we compared the values of the full set of MIs with results obtained using the numerical package AMFlow [69,70], always finding agreement with the requested full precision (16 digits).
The bare two-loop 2 → 2 amplitudes need UV renormalisation.We employ two different schemes.In the first, we renormalise the external fields on-shell, the strong coupling constant in a mixed scheme in which the light-flavour contribution is renormalised in MS while the heavy quark is renormalised at zero momentum, and the Yukawa coupling and the heavy-quark mass in OS.In the second scheme, we renormalise both the Yukawa coupling and the heavy-quark mass in MS.
Once the UV divergences are renormalised, the twoloop 2 → 2 amplitudes still contain residual IR divergences that appear as poles in 1/ .The structure of these IR poles is universal and is described by factorization formulae [51].
We examined the behaviour of the two-loop 2 → 2 amplitudes in the soft and collinear limits of one unresolved parton.To do so, we generated sequences of phase space points tending to the desired limit and compared the value of the amplitudes to that predicted by the corresponding one-loop IR factorization formulae [71][72][73][74][75][76].We found that the amplitudes approach the limit formulae with the rate expected from the cancellation of the leading singularity order-by-order in .We observed the expected behaviour of the full amplitudes independently of the value of internal quark mass as well as for the interference of two massive internal quarks.
As a further check, we verified that in the limit of very large transverse momentum, the full amplitude is in reasonable agreement with the approximated results computed in ref. [77].

THE CROSS SECTION
The results for the hadronic cross sections have been obtained by combining the integration of the one-loop 2 → 2 (Born) and 2 → 3 (real) amplitudes, including the subtractions and the initial-state collinear remnants, with the integration of the two-loop 2 → 2 (virtual) amplitudes.To efficiently integrate the virtual part we have first produced a grid using the virtual correction in the HEFT.To increase the precision in the computation of the tail of the Higgs p t distribution, we have biased the weight with an appropriate exponential factor and generated a second grid.Subsequently we have used these grids to integrate the amplitudes in the full theory.In particular for every choice of scales we have run about 14k points on the cross section grid and 16k on the biased one.When we use the MS mass renormalisation, the value of the internal quark masses changes dynamically in each phase space point so that no form factors or integrals can be recycled from the evaluation of another point.The run time of the evaluation of the two-loop matrix element for one kinematic configuration and one choice of masses circulating in the loops varies greatly depending on the phase-space point and the existing precomputed points that DiffExp can transport from.However, we can qualitatively report that it ranges from 5 minutes to one hour, with a median around 15 minutes.We implemented our two-loop amplitudes through a MadGraph5 aMC@NLO plugin, similar to what was done in ref. [36], where the amplitude appeared as an effective vertex in a UFO [78] model.In that representation, the form factors are understood as dynamical couplings which will be re-evaluated for each phase-space point by the tree-level Fortran matrix element code via an interface to Mathematica.The main advantage of such an implementation is that it greatly facilitates the distribution and reproducibility of our two-loop amplitude as it can now readily be generated like any tree-level process in a version of MadGraph5 aMC@NLO equipped with our plugin.It moreover offers additional flexibility regarding the selection of the type of contribution the user is interested in.This includes 1) the possibility of choosing between OS and MS mass renormalisation schemes, 2) selecting particular interference terms and which massive quark flavours are included in each fermion loop, 3) potentially attaching Higgs decay structures to the production process and finally 4) computing renormalisation scale variations through a reweighting procedure necessitating only a single evaluation of the form factors.
For the case of OS top-quark mass renormalisation, we validated this implementation by aligning our setup with that chosen in refs.[28,29] and found great agreement both at the inclusive and differential level.*

RESULTS
In this section we show a selection of predictions obtained using our computation.For our simulation we use the value G F = 1.16639 • 10 −5 GeV −2 for the Fermi constant and the NNPDF40 nlo as 01180 [79] PDFs set and α S .Jets are reconstructed using the anti-kt algorithm [80] with resolution variable R = 0.4.We select events of proton-proton collisions at √ s = 13 TeV where the Higgs boson is produced in association with a jet with transverse momentum larger than p j1 T > 20 GeV.The masses of the Higgs boson and of the internal quarks are set to m H = 125.25 GeV, and m OS t = 172.5 GeV when considering on-shell top mass renormalisation, whereas we compute GeV in the case of MS top-quark mass renormalisation, using leading logarithmic evolution.We also included the effect of having an internal massive bottom quark attached to the Higgs boson.For a consistent use of the decoupling scheme with five-flavour evolution of the strong coupling and PDFs we have chosen to retain the bottom mass Table I.Cross section for Higgs boson production in association with a jet with transverse momentum larger than p j 1 T > 20 GeV, with top-and bottom-quarks circulating in the heavy-quark loop, in the MS scheme (upper line); with top-quark circulating in the loop, in the MS scheme (middle line) and in the OS scheme (lower line), at LO and NLO accuracy.
only in the loop where it is coupled to the Higgs boson and set it to zero in the other pure QCD contributions.The bottom-quark mass and the corresponding Yukawa coupling are computed evolving 18 GeV (also with leading logarithmic accuracy).Our choice for the central renormalisation (µ R ) and factorization (µ F ) scales is with the sum running over all partons in the final state.
As for the scale variations, we take the envelope of the seven scale choices obtained by varying the central scales by a factor of two in both directions excluding the maximal distance.
Our results for the semi-inclusive Higgs production cross section within our jet acceptance cut are given in Table I, where we note that the relative scale uncertainty goes down from about 30% at LO to about 14% at NLO.Further, by comparing the first two lines, we see that the top-bottom interference yields a negative contribution at LO, and instead a positive contribution at NLO which cancels the offset between the cross section with and without top-bottom interference.
The bottom quark is expected to affect the production rate only at small and intermediate values of the Higgs p T distribution.Approximated results for the effect of the bottom quark at intermediate p T values have been previously reported in [34].In fig. 1, we show the results of our exact computation for the Higgs p T distribution in 20GeV-wide bins at intermediate p T values, at LO (left panel) and NLO (right panel) accuracy, with topand bottom-quarks circulating in the heavy-quark loops in MS (black curve); with top-quark only, in MS (blue curve) and in OS (red curve).Vertical bars represent the statistical error from the Monte Carlo integration.In the left panel, the first bin is empty, since the LO kinematic contraint, p T = p j T , forces the Higgs p T to be not lower than 20GeV; in the second bin, the contribution of the top-bottom interference is negative, and fades away from the third on.In the right panel, the NLO contribution of the top-bottom interference is negative in the first bin, positive in the second, negative again, but smaller, in the third, and dies out from the fourth on.Note that the seven-point scale variation (not shown) provides a much larger uncertainty than the difference among the bars that are shown in the figure.Also note that at LO, the impact of the change of top-quark mass renormalisation scheme is almost indiscernible in the second bin of the left panel.However, at NLO this impact is far greater, and overall 15 times larger than at LO in our semi-inclusive cross-section of tab.I.We leave further investigation of this enhanced sensitivity to the mass and Yukawa renormalisation scheme at NLO to future work.While the scale uncertainty is expected to be reduced in a resummed calculation of the Higgs p T , we anticipate that the difference among the different predictions will persist.Further investigation of the Higgs p T distribution in the low-energy limit is beyond the aim of the present work and will be presented elsewhere.In fig.2, we plot the Higgs p T distribution in 25GeVwide bins, with top-and bottom-quarks circulating in the heavy-quark loops in MS, at LO (green curve) and NLO (red curve) accuracy.The scale uncertainty bands at LO (yellow band) and NLO (purple band) accuracy are obtained by taking the envelope of seven-point scale variations.Not to clutter the plot, we refrain from showing the same distribution with the top-quark only, either in MS or in OS, opting for highlighting their behaviour in the next figures.In fig.3, we plot the ratio of the Higgs p T distribution at NLO over the same at LO in 50GeV-wide bins, with top-and bottom-quarks circulating in the heavy-quark loops, in MS (upper panel); with top-quark only, in MS (middle panel) and in OS (lower  We note that except for the first bin all ratios have a numerical value greater than or equal 2. In particular, the curves of the upper and middle panels have a similar, rather flat, shape with a numerical value which is larger than 2 on most of the p T range, while the curve of the lower panel wiggles about the value 2. In figs. 4 and 5, we plot the ratio of the Higgs p T distribution, with top-and bottom-quarks circulating in the heavy-quark loops, over the distribution with the topquark only, both in MS (upper panel); the ratio of the distribution with the top-quark in OS over the one with top-and bottom-quarks in MS (middle panel); the ratio of the distribution with the top-quark in OS over the one with the top-quark in MS (lower panel).The scale uncertainty bands as reported in the y-labels are given by the ratio of the bands of the distributions over their and 5, we note that the Higgs p T distribution with the top-quark only in MS falls off faster than the same distribution in OS, as p T increases, the more so at LO than at NLO accuracy.This can be understood by the fact that the top-quark mass has an OS fixed value, while the running top-quark mass decreases as the p T values, and thus the renormalisation and factorisation scales, eq. ( 1), increase.This is at the origin of the difference between the upper/middle and the lower panels of fig. 3 at high p T values.

CONCLUSIONS
Building on two-loop amplitudes for Higgs + three partons [44], which are valid for arbitrary quark masses circulating in the heavy-quark loops, we have computed for the first time the NLO QCD corrections to the Higgs p T distribution in Higgs + jet production via gluon fusion, with top and bottom quarks circulating in the heavyquark loops.The exact mass dependence on the top and bottom quarks has been included using, for the first time in this context, a running mass renormalisation scheme, the MS scheme.We have also provided predictions for the Higgs p T distribution with only the top quark, in MS and OS schemes.
We find that within the scale uncertainty the LO contribution of the bottom quark, and thus of the topbottom interference, to the Higgs boson production is almost erased at inclusive level by the NLO corrections.On the other hand, at the low end of the p T distribution, the interference induces a non trivial change of shape.However, for precision studies on the high energy tail and with the current attainable accuracy, the use of the Higgs p T distribution with only the top quark circulating in the heavy-quark loops is fully justified.Finally, we find that the Higgs p T distribution with the top-quark only in MS falls off faster than the same distribution in the OS scheme as p T increases.This would have an obvious impact on any numerical study, requiring then that the choice of mass renormalisation scheme be done with great care.

Figure 1 .
Figure 1.Higgs pT distribution in the intermediate pT range.

Figure 2 .
Figure 2. Higgs pT distribution with top-and bottom-quarks.

Figure 3 .
Figure 3. NLO/LO ratio of the Higgs pT distribution.

Figure 4 .
Figure 4. Ratio of Higgs pT distributions at LO.

Figure 5 .
Figure 5. Ratio of Higgs pT distributions at NLO.