ZH production in gluon fusion at NLO in QCD

We present fully differential next-to-leading order results for Higgs production in association with a $Z$ boson in gluon fusion. Our two-loop virtual contributions are evaluated numerically using sector decomposition, including full top-quark mass effects, and supplemented at high $p_T$ by an analytic high-energy expansion to order ($m_Z^4, m_H^4, m_t^{32}$). Using the expanded results we also present a study of the top-quark mass scheme uncertainty at large $p_T$.


Introduction
Higgs boson production in association with a Z boson is a particularly interesting process as it probes both the Higgs boson coupling to Z bosons as well as to fermions. Despite its relatively small cross section, the process pp → ZH in combination with the decay channel H → bb was the "discovery channel" for Higgs boson couplings to bottom quarks [1,2], as while inclusive H → bb suffers from large backgrounds, the Z boson offers straightforward triggering. Recently, associated Higgs production was also used to place limits on the Higgs-charm coupling [3,4].
The loop-induced gluon channel formally enters at next-to-next-to-leading order (NNLO), with respect to the pp → ZH process. However, due to the dominance of the gluon parton distribution function (PDF) at the LHC, this channel is sizeable; it contributes about 6% to the total NNLO cross section and becomes significant in the boosted Higgs regime for p H T 150 GeV [5][6][7][8][9]. For more details about calculations of higher-order corrections to the pp → ZH process we refer to Ref. [10]; here we focus on the gluon channel.
Experimental measurements of ZH production [11][12][13][14][15][16][17] still suffer from large statistical uncertainties, however the statistics will improve considerably in LHC Run 3 and at the High Luminosity LHC. On the theoretical side, the scale uncertainties for this process are large. They are dominated by the gg → ZH channel which, being loop induced, enters at its leading order (LO) into the simulation programs [18][19][20][21][22][23][24][25] used by the experimental collaborations. Therefore next-to-leading order (NLO) QCD corrections to the gg → ZH process, calculated at LO in Ref. [26], are important; particularly so in view of providing constraints on anomalous couplings, because of the sensitivity of this process to both the Higgs couplings to fermions as well as to vector bosons. Furthermore, it provides a way to put constraints on the sign of the top quark Yukawa coupling as well as on its CP structure [7,20,21,27], because the cross section has a component where this coupling enters linearly. Recently, it was shown [28] that this process also has the potential to probe anomalous Zbb couplings and thus to shed light on the long-standing discrepancy of the forward-backward asymmetry A b F B measured at LEP, with the Standard Model prediction.
The NLO QCD two-loop amplitude for gg → ZH, including full top-quark mass effects, has been calculated numerically in Ref. [29]. The two-loop amplitude has also been calculated based on high-energy (also sometimes called "small-mass") and largem t expansions, supplemented with Padé approximants to improve the description beyond the high-energy radius of convergence [30]. Results for the virtual corrections based on a transverse-momentum expansion are also available [31]. Recently, the virtual corrections to both gg → HH and gg → ZH, based on a combination of transversemomentum expansion and high-energy expansion, have been presented in Ref. [32]. The total cross section and invariant-mass distribution for gg → ZH at NLO QCD has been calculated in Ref. [33], based on a small-(m Z , m H ) expansion [34] of the two-loop amplitude, retaining the full top-quark mass dependence. An expansion of the virtual and real NLO contributions to gg → ZH for large m t has been computed in [35].
In this work we present full NLO QCD results for the gg → ZH process where the two-loop amplitude is based on a combination of the numerical results of Ref. [29] and an extended version of the results from the high-energy expansion of Ref. [30], thereby providing reliable and accurate results in all kinematic regions, in particular in the boosted Higgs regime which is particularly sensitive to new physics effects. A similar combination has already been carried out successfully for Higgs boson pair production [36]. In addition, we consider two renormalisation schemes for the top quark mass, the on-shell (OS) scheme and the modified Minimal Subtraction (MS) scheme, and investigate how the scheme dependence impacts the phenomenological results.
This paper is organised as follows. In Section 2 we describe the calculation of the NLO corrections and the combination procedure. Section 3 contains results for the total cross section at different center-of-mass energies as well as ZH invariant-mass and transverse-momentum distributions. The second part of Section 3 is dedicated to the Figure 1: Representative Feynman diagrams for the virtual correction to the ggZH amplitude. We neglect the masses of all quarks except the top quark; therefore, due to the Yukawa couplings, only the top quark gives a non-zero contribution for the box diagrams. All quark flavours contribute for the triangle diagrams, however the contribution from each massless generation is zero due to a cancellation between the up-type and down-type quarks. We calculate in the Feynman gauge and so also include the set of diagrams where the Z-boson propagators are replaced by Goldstone bosons. discussion of the top-quark mass scheme dependence, before we conclude in Section 4.

Setup of the calculation
In this section we summarise the computation of the individual contributions to the cross section at NLO and describe the combination of the virtual corrections computed in [30] and [29].

Virtual two-loop contributions
The calculation of the renormalised and infrared (IR) subtracted two-loop amplitude, called V, is described in detail in Refs. [30] and [29]. For completeness we repeat in the following the most important steps.
In Ref. [30] the amplitude of the process gg → ZH has been written as a linear combination of six form factors. At one-loop order it is straightforward to obtain exact results for the form factors. At two-loops expansions for large and small top-quark masses were performed. In this work only the high-energy expansion, for which m 2 H , m 2 Z m 2 t s, |t|, is of relevance. In Ref. [30] an expansion up to order (m 2 Z , m 2 H , m 32 t ) was computed. In this work we extend that result up to quartic order (m 4 Z , m 4 H , m 32 t ) (including also the "mixed" quartic term m 2 Z m 2 H ) and show that including these quartic terms improves the agreement with numerical results. LiteRed [37] is used to expand the integrals appearing in the amplitude, followed by an integration-byparts (IBP) reduction to master integrals using FIRE [38]. The reduction relations are substituted into the amplitude and simplified using FORM [39]. In Ref. [30] it was shown that Padé approximants for the expansion in m t significantly improve the description beyond the radius of convergence of the naive expansion and reliable results can be obtained for p T ≥ 150 GeV (see Fig. 2); we apply the same procedure here. Note that in this approach one can construct the high-energy expansion of V as a function of m H , m Z and m t which allows for a variation of these parameters. Furthermore, it is straightforward to perform a scheme change and convert the top quark mass from the OS to the MS scheme. The subsequent construction of the new Padé approximants requires only negligible CPU time.
In Ref. [29], the two-loop amplitudes have been calculated via a projection onto a basis of linear polarisation states as suggested in Ref. [40]. For the IBP reduction of the resulting form factors to master integrals we used Kira [41][42][43] in combination with the rational function interpolation library FireFly [44,45]. For the numerical integration, using a quasi-finite basis [46] of master integrals is beneficial. This basis is related to the default basis by dimension shifts and higher powers of propagators (dots). For the derivation of the set of dimensional recurrence relations which connects integrals in D + 2n and D dimensions, we used LiteRed [37] and Reduze [47].
To evaluate the master integrals, we applied sector decomposition as implemented in the program pySecDec [48,49], using a quasi-Monte Carlo algorithm [49,50] for the numerical integration. In particular, we made use of one of the new features of pySecDec [51] to integrate a weighted sum of integrals such that the number of sampling points used for each integral is dynamically set, according to its contribution to the total uncertainty of the amplitude.
We renormalise the strong coupling in the MS scheme with 5 active quark flavours. The top quark mass is renormalised in either the OS scheme or the MS scheme, as indicated. Expanding each renormalised form factor A i=1,...,n in powers of the strong coupling, a s = α s /(4π), we may write and then obtain IR-finite amplitudes using the Catani-Seymour subtraction operator [52] The squared 2 → 2 amplitude, in the helicity basis, can be written as  where the squared Born amplitude (B) and the Born-virtual interference (V), are given by The sum/average, denoted by , runs over all helicities and averages over the incoming spin and colour indices.

Combination of the two approaches
In Fig. 2 we show the difference between the two calculations of the two-loop virtual amplitudes, relative to the LO amplitude. This difference is independent of the IR subtraction scheme used for removing the IR singularities. The plot shows that the Padé-improved high-energy expansion converges on the full result for p T 150 GeV. However, using only the quadratic terms in m 2 Z and m 2 H in the expansion, a difference at the two permill level remains even at large p T . After including the quartic terms in the expansion, most of the points with p T > 200 GeV agree within the numerical uncertainty at the 2 · 10 −5 level (with a few outliers with large numerical uncertainty reaching up to the 2 · 10 −3 level). At low p T the differences increase, reaching up to 0.15% at 200 GeV and up to 2.8% at 150 GeV.
Since the two results are consistent for sufficiently large p T , in the following we use the results based on the high-energy expansion for contributions with p T > 150 GeV and use the numerical evaluation of the amplitude only for p T < 150 GeV. The two contributions are then combined at the histogram level. Our results are therefore valid in all phase-space regions, but avoid the costly numerical evaluation in large parts of the phase space.

Phase-space sampling
The integration of V over the phase space is achieved by a reweighting procedure based on the Born events. Specifically, we use the events of a LO calculation and apply the accept-reject method to obtain a list of sampling points for the virtual contribution distributed according to the probability density function ∼ L gg, where L gg,0 is the gluon-gluon luminosity as defined in Ref. [53] and B 0 is the LO matrix element as used in the LO calculation. The factor dPS is the Jacobian of the phasespace integration and the function f (p T , m ZH ) can be used to enhance the number of sampling points in specific regions. Choosing, e.g., f = f (m ZH ) ∝ dσ B /dm ZH leads to sampling points which are uniformly distributed in m ZH , thus enhancing the number of events in the tail of the distribution, whereas f (p T , m ZH ) = 1 results in sampling points distributed according to the fully differential LO cross section.
The virtual contribution to the cross section is then given by where L gg and V are the new gluon-gluon luminosity and virtual matrix elements, whereas L gg,0 , B 0 and σ B,0 are obtained from the original LO calculation, with the total LO cross section σ B,0 .
In the following results we use three different sets of sampling points, optimized for the total cross section, as well as the m ZH and p T distributions. While the exact form of f (p T , m ZH ) is not important, a good choice can be obtained with a fit using a Padé ansatz. In the region p T , m ZH ≥ 2 TeV, where a uniform sampling of points is not needed, we keep f constant. In total, we use 1294 numerically evaluated points distributed according to the differential LO cross section. For the m ZH and p T distributions, we combine these results with sets containing an additional 6000 points, optimised for the corresponding distribution, evaluated using the Padé-improved highenergy expansion.

Computation of the real radiation contributions
The real radiation matrix elements are calculated using the one-loop amplitude generator GoSam [54,55] together with an in-house C++ code, similar to the one used in Refs. [53,56], where the IR singularities are subtracted in the Catani-Seymour scheme [52], supplemented by a dipole phase-space cut parameter α cut [57]. We have checked that our implementation of the dipoles reproduces the matrix element in the soft and collinear limits and that our results are independent of α cut for 0.2 ≤ α cut ≤ 1.
To check the numerical precision of our real matrix elements we use several rotation tests (i.e., we perform azimuthal rotations about the beam axis and recompute the phase-space point). We first compute the matrix element at a given phase-space point and a rotated phase-space point in double precision. If the results do not agree to 10 digits, we compute the phase-space point in quadruple precision and check if it agrees with the double-precision evaluations to 7 digits. If the results do not agree we compute a rotated point in quadruple precision and check that the quadruple-precision results agree to 10 digits; a vanishingly small fraction of points failed this test and were discarded. Using the above procedure, we did not find it necessary to apply a technical cut to our 2 → 3 phase space.
In Fig. 3 we show examples of the Feynman diagrams included in our real radiation. We include all diagrams appearing in the ggZHg and qqZHg amplitudes (as well as their crossings) which contain a closed fermion loop and have either a Z-boson or Goldstone boson coupled to that loop. We consider n f = 5 massless quarks and a massive top quark running in the fermion loops. Due to the presence of the Yukawa coupling in the upper diagrams, the massless quarks give a non-zero contribution only for the lower row of diagrams. For the diagrams in the second and third column of the lower row, the contribution of each massless generation is zero due to a cancellation between the up-type and down-type quarks. We calculate in the Feynman gauge and so also include the set of diagrams in which the Z-boson propagators are replaced by Goldstone bosons.
In Fig. 4 we show examples of Feynman diagrams which are not included in this work; these diagrams do not have a Z boson or Goldstone boson coupled to the closed fermion loop and so belong to the Drell-Yan class of diagrams, which we do not consider. The class of diagrams with the Higgs boson coupled to a closed quark loop, represented by the left figure, is UV/IR finite and separately gauge invariant; it was considered in detail in Ref. [58].

Total and differential cross sections
Our results are based on a center-of-mass energy of √ s = 14 TeV unless stated otherwise. We use the NNPDF31 nlo pdfas parton distribution functions [59] with masses determined by the ratios m Results for the total cross section with full top-quark mass dependence at three different center-of-mass energies, including scale uncertainties resulting from the 7-point scale variation, µ R,F = ξ R,F m ZH with ξ R,F = (2, 2), (2, 1), (1, 2), (1, 1), (  Differential results for the invariant mass m ZH = (p Z + p H ) 2 of the Z-Higgs system are shown in Fig. 5 for the central scale choices m ZH and H T , with where the sum runs over all final state massless partons k. For the fully-inclusive case (left), the K-factor is relatively flat with a value of about two, except at very low invariant masses where threshold corrections are significant. The kink in the distribution at m ZH 350 GeV is related to the tt-production threshold. Only a small reduction of the scale uncertainty is observed going from LO to NLO. Note that the quark-gluon channel for this process first opens up at the NLO level. The cuts p T,H ≥ 140 GeV, p T,Z ≥ 150 GeV (Fig. 5 (right)) somewhat decrease the K-factor.
The Z-boson transverse momentum distributions at LO and NLO are shown in Fig. 6. In the left plot we observe a K-factor which rises with increasing p T,Z , reaching a value of almost 5 at p T,Z = 1 TeV, it is only slightly tamed by the cuts on p T,H and p T,Z (right plot). Fig. 7 shows the Higgs-boson transverse momentum distributions with and without p T cuts. In the inclusive case (left) an extreme rise of the K-factor with increasing p T,H , up to values of about 20 towards p T,H = 1 TeV, is observed. The cuts p T,H ≥ 140 GeV, p T,Z ≥ 150 GeV decrease this K-factor by a factor of about 3 at large p T,H values. The cuts have such a large effect on the K-factor of this distribution as they remove configurations with a hard jet recoiling against a relatively hard Higgs while the Z boson is soft, this configuration dominates the tail of the distribution but is not present at   LO. This behaviour was already reported in Ref. [20] and traced back to diagrams with t-channel gluon exchange, it was further studied in Ref. [60]. The reason why the rise of the K-factor is more pronounced in the p T,H case than in the p T,Z case can be related to the coupling structure of the Z and Higgs bosons to top quarks. In the diagrams where both the Higgs and the Z boson are radiated from a top quark loop, the probability to radiate a "soft" Z boson while the Higgs boson recoils against a hard jet is related to  the soft Eikonal factor p µ /(p·p Z ), where p µ generically denotes the radiator momentum. The probability to radiate a "soft" Higgs boson on the other hand is proportional to m t /(p · p H ). The ratio of these Eikonal factors is p T /m t 1, thus at large transverse momentum p T of the radiator it is more likely that the Z boson is soft and the Higgs boson is hard.

Investigation of different top quark mass renormalisation schemes
We now turn to the discussion of the uncertainties stemming from the use of different top quark mass renormalisation schemes. Such uncertainties have been investigated in detail for the case of Higgs boson pair production in Refs. [60][61][62][63]. For top quark pair production at NNLO, scheme uncertainties have been studied in Ref. [64]. Top quark renormalisation scheme uncertainties also have been investigated for NLO ttH [65] and ttj [66] production at the LHC, as well as for off-shell Higgs production and LO Higgs+jet production [60].
In this section we investigate the top-quark mass renormalisation scheme dependence of ZH production. For this purpose we convert the top quark mass to the MS scheme, which is an appropriate renormalisation scheme in the high-energy region. It is thus sufficient to perform the scheme change in the analytic high-energy expansion of the virtual corrections, where it is straightforward to obtain the corresponding analytic expressions by making the replacement  and to apply the Padé procedure as described in Section 2.1. Afterwards the result is combined with the real-radiation contributions where the MS top quark mass is used in the numerical evaluation.
We make three different choices for the central value of the top quark renormalisation scale, µ t , T,i + k |p T,k | (k sums over massless partons) and vary µ t up and down by a factor of 2 to obtain an uncertainty estimate. For the conversion of the numerical values of the top quark mass between the OS and the MS schemes we proceed as follows: we first convert the top quark OS mass to the MS scheme at the scale µ t = m t , at four-loop accuracy. For our input values, m t = 173.21 GeV and α s (m Z ) = 0.118, this gives m t (m t ) = 163.39 GeV. We then use the renormalisation group equation, at five-loop accuracy with six active quark flavours, to run from µ t = m t to the desired renormalisation scale for m t . For both the numerical scheme conversion and the running we use the Mathematica and C++ codes RunDec and CRunDec [67,68].
At LO the difference between the schemes is purely parametric and is driven both by the top quark mass appearing in the propagators and the Higgs-top Yukawa coupling. At NLO the OS and MS schemes differ by the parametric choice of m t and a shift proportional to the derivative of the LO, i.e., the mass counterterms, which partly compensates the parametric difference.
In Fig. 8 we show predictions at LO (dashed lines) and NLO (solid lines) for µ R = µ F = m ZH with three different choices of the top-quark renormalisation scale, µ t . The red band is generated by varying the scale µ t = m ZH up and down by a factor of 2, keeping µ R = µ F = m ZH fixed. We observe that the scheme choice for the top quark mass has a large impact on the predictions for both the invariant mass distribution (with p T,H ≥ 140 GeV and p T,Z ≥ 150 GeV cuts) and the p T,Z distribution. Focusing on the invariant mass distribution, we observe that at LO for m ZH ∼ 1 TeV the OS result is approximately a factor of 2.9 times the MS result with µ t = m ZH . At NLO the difference between the schemes is somewhat reduced, with the OS result around 1.9 times the MS result with µ t = m ZH , see Table 2. Taking, for example, the difference between the OS and the MS result with µ t = m ZH as a mass scheme uncertainty would result in a +0% −65% uncertainty at LO and a +0% −47% uncertainty at NLO for m ZH = 1 TeV. Alternatively, taking the MS result with µ t = (m ZH /2, m ZH , 2 m ZH ) as an uncertainty gives +26% −21% at LO and +17% −14% at NLO for m ZH = 1 TeV. We observe a similar pattern for large p T,Z , the difference between the OS scheme and the MS with µ t = m ZH scheme at p T,Z ∼ 1 TeV is reduced from a factor of 2.8 at LO to a factor of about 1.9 at NLO.
The K-factor, defined as the ratio of the NLO result in a given scheme to the LO result in the same scheme, is typically larger in the MS scheme than in the OS scheme; this feature is also observed in Higgs pair production [63]. The K-factor of the invariant mass distribution is relatively flat for all scheme choices, with the dynamic scale choices µ t = H T and m ZH yielding K ∼ 2.5−2.7 while the OS scheme has K ∼ 1.6 for m ZH ∼ 1 TeV. The MS scheme with µ t = m t (m t ), where the logarithm appearing in Eq. (3.2) vanishes, has a very similar K-factor to the OS scheme; this differs from the HH case where the two schemes had a similar shape but a different normalisation.
For the p T,Z distribution the pattern of K-factors for the different schemes is broadly the same as for the invariant mass distribution, but in all cases the K-factors rise with p T,Z reaching up to K = 5 for dynamic µ t choices at p T,Z = 1 TeV.
Comparing the results obtained here for gg → ZH to other loop-induced processes, such as off-shell Higgs production, Higgs pair production and Higgs plus jet production, we note that the ZH process has a larger mass scheme dependence at LO. For off-shell Higgs production and Higgs pair production going from LO to NLO approximately halves the uncertainty due to the mass scheme choice; in the ZH case we also observe a reduction in the uncertainty, but by less than a factor of 2.
In the HH case, in the high-energy limit, the triangle contribution is suppressed by a factor of 1/s w.r.t. the box form factors. Here, the leading high-energy behaviour of the box form factors has the form [62,69] where the log[m 2 t ] term in A (1) i is due to the renormalisation of m t , and the overall power of m 2 t comes from the Yukawa couplings. Converting to the MS scheme using Eq. (3.2) results in a logarithm of the form log[µ 2 t /s]. In Ref. [62] it was argued that choosing µ 2 t ∼ s minimizes these logarithms and is thus the preferred central scale choice of the Yukawa couplings. However, in the present ZH case, the structure is different. Firstly, the triangle contribution is not suppressed w.r.t. the box form factors, and secondly logarithms involving m t appear in the box form factors already at leading order. Unlike in the HH case, where the overall power of m 2 t in Eq. (3.3) comes entirely from the top Yukawa couplings, in gg → ZH one of the overall m t factors must come from the top-quark propagators, hence the leading term in the small-mass expansion is already power-suppressed by one power of m t . Similar, power-suppressed, mass logarithms have been studied in the context of single Higgs production, see for example Ref. [70] and references therein. The leading helicity amplitudes for ZH in the high-energy limit have the form

Conclusions
The gg → ZH channel contributes to the pp → ZH process starting at NNLO, accounting for around 6% of the total cross section. However, the gluon-fusion channel suffers from a large scale dependence at LO and is a significant source of theoretical uncertainty for Z boson production in association with a Higgs boson at the LHC [1,2,[12][13][14][15][16].
In this work, we have presented the complete NLO corrections for the loop-induced gluon-fusion channel; they increase the gluon-fusion cross section by about a factor of 2, and reduce the scale dependence. We have investigated the invariant mass distribution and transverse momentum distributions for both the Z boson and Higgs boson. At large transverse momentum, we found that the NLO corrections can be very large, more than 10 times the LO result for p T,H ; the origin of this behaviour can be traced back to extremely large real radiation corrections when a soft Z boson is radiated from a top quark loop [20].
We have also studied the top quark mass scheme uncertainties for this channel, i.e., the difference between results produced with the top quark mass renormalised in the on-shell scheme and the MS scheme. As for other loop-induced processes with a scale above the top quark pair-production threshold [60][61][62][63], we found a large mass scheme uncertainty at LO. At NLO the mass scheme uncertainty is smaller than at LO, but remains at least as large as the usual renormalisation and factorisation scale uncertainties.
The inclusion of the NLO corrections to the gluon-fusion channel is essential for correctly describing ZH production at the LHC and HL-LHC. The size of the NLO corrections, especially for the transverse momentum distributions, and the mass scheme uncertainty motivate a calculation of gg → ZH beyond NLO and the study of this process beyond fixed order.
[14] ATLAS collaboration, Combination of measurements of Higgs boson production in association with a W or Z boson in the bb decay channel with the ATLAS experiment at √ s = 13 TeV, tech. rep., CERN, Geneva, Sep, 2021.
[15] ATLAS collaboration, G. Aad et al., Direct constraint on the Higgs-charm coupling from a search for Higgs boson decays into charm quarks with the ATLAS detector, 2201.11428. [17] CMS collaboration, Combined Higgs boson production and decay measurements with up to 137 fb-1 of proton-proton collision data at sqrts = 13 TeV, tech. rep., CERN, Geneva, 2020.