On the NLO QCD Corrections to Gluon-Initiated ZH Production

We compute the QCD corrections at next-to-leading order for the process $gg \rightarrow ZH$, including both the virtual two-loop terms and real-emission contributions. The two-loop box diagrams in the virtual corrections are approximated analytically over the complete phase space, combining the results of an expansion in the limit of small transverse momentum and an expansion in the regime of high energy. We obtain both inclusive and differential results for the cross section. We find that the NLO QCD corrections are of the same size as the LO contribution up to $ZH$ invariant masses close to 1 TeV, but they increase significantly when higher energies are considered, due to a class of real-emission diagrams in which the $Z$ boson is radiated from an open quark line. Finally, we estimate the uncertainty due to the renormalization scheme used for the top-quark mass both on the total and differential cross section.


Introduction
In the last ten years following the discovery of the Higgs boson [1,2] a great experimental effort has been undertaken at the Large Hadron Collider (LHC) in order to measure the properties of this particle with high precision. The increased statistics expected from the next LHC runs and from the High-Luminosity phase will allow for more accurate comparisons with the Standard Model (SM) predictions, and refined theoretical calculations are needed to extract as much information as possible from the forthcoming measurements 1 .
The process in which a Higgs is produced together with a weak vector boson, known as V H associated production (V = W, Z), is of great relevance at the LHC not only as a probe of the couplings between the Higgs and the weak bosons, but also because of its sensitivity to the H → bb decay. Indeed, the large QCD backgrounds affecting searches of H → bb can be reduced more efficiently in V H than in other production modes [4][5][6].
Recently V H production has been considered also to improve the constraints on the charm Yukawa coupling [7].
In current V H analyses the impact of theoretical uncertainties from missing higher-order terms in perturbative calculations depends on the choice of the final-state boson V : while in the W H case these uncertainties are around 1%, and are comparable to the uncertainties from the PDF and α s determination, the corresponding uncertainties for ZH production are slightly larger, close to 3% [8]. In view of the progress in experimental precision it is hence essential to improve the theoretical control over pp → ZH.
There are two partonic channels that contribute to the pp → ZH hadronic cross section. The qq-initiated (Drell-Yan-like) channel gives the leading contribution, and higherorder corrections are well under control as they are known to next-to-next-to-leading order (NNLO) in QCD [9][10][11] and to NLO in the EW interactions [12,13]. The gg-initiated channel has been computed at LO in Refs. [14,15], and it contributes for the first time as a NNLO QCD correction to the hadronic process. However, this sub-leading channel is enhanced by the large gluon luminosity at the LHC, and it provides about 6% of the total cross section for 13 TeV collisions. Differential analyses have pointed to significant differences in the shapes of distributions between the qq and the gg channels [16][17][18][19], showing an increased relative importance of the gg channel in the boosted regime [20]. The gg → ZH process has been considered also as a probe of new physics effects, with examples including anomalous couplings [19][20][21][22] and new degrees of freedom [20,23]. Finally, the gluon-initiated channel has the largest impact on the theoretical uncertainties for pp → ZH, as only the LO contribution is included in the Monte Carlo programs used for experimental studies. Since an equivalent gluon-initiated contribution is not allowed in W H production because of electric charge conservation, gg → ZH is responsible for the larger uncertainties in ZH production compared to W H. The reduction of the current theoretical uncertainties requires the calculation of the NLO QCD corrections to gg → ZH, which are the main subject of this paper.
Since the main contribution to gg → ZH at LO comes from one-loop diagrams involving loops of massive quarks, the most challenging part of the NLO calculation is associated to two-loop multi-scale integrals in the virtual corrections. An estimate of these terms has been obtained in Ref. [24] in the m t → ∞ limit and included in the program vh@nnlo [18,25], while the effects of a large but finite top-quark mass have been considered in Refs. [26,27]. An approximation that retains the effects of the top-quark mass and that is accurate over the complete phase space has been presented in Ref. [28], whereas the very recent calculation of Ref. [29] relies on a combination of the Padé-improved high-energy expansion presented in Ref. [30] and a numerical evaluation using sector decomposition [31]. All of the previous results show that the NLO corrections are comparable in magnitude to the LO contribution 2 , and this fact makes the implementation of the NLO terms in a Monte Carlo code (possibly interfaced with parton-shower generators [32,33]) a priority for a reliable interpretation of the experimental results 3 . In this paper we take a step towards this goal, as we provide a fast and flexible way to compute gg → ZH at NLO in QCD including the effects of a finite top-quark mass .
We obtain a very reliable approximation of the two-loop virtual corrections by merging analytic results that are accurate in two complementary phase-space regions, namely the results of the transverse-momentum expansion of Ref. [35] and the high-energy (HE) expansion of Ref. [30]. In Ref. [36] it has been shown that this approach can provide a fast evaluation of the virtual corrections with an accuracy of 1% or below. In this paper, we present a complete assessment of the NLO QCD corrections to gg → ZH by including the real-emission contributions, which are related to one-loop diagrams with an additional parton in the final state. We use our results to quantify the effects of the gluon-initiated channel on the hadronic cross section, both at the inclusive and differential level, and we compare our findings with previous independent calculations. Additionally, we discuss the theoretical uncertainty related to the choice of renormalization scheme for the top-quark mass. We finally point out a so-far unnoticed feature of the NLO QCD corrections, as we observe that the contribution from the class of 2 → 3 diagrams in which the Z boson is emitted from an open quark line becomes numerically relevant for collisions at very high energies, namely for invariant masses in the range M ZH > 1 TeV. To the best of our knowledge, the impact of these diagrams on ZH production over this energy regime has not been discussed in the literature. This paper is structured as follows: in the next section we set our notation and we describe the ingredients for the calculation of the NLO corrections to gg → ZH; in section 3 we present our results for the inclusive and differential cross section and we assess the top-mass scheme uncertainty. In section 4 we state our conclusions.

The gg → ZH channel at NLO
The cross section for the subprocess gg → ZH + X in the hadronic reaction pp → ZH + X at center-of-mass energy √ s, can be written as where M 2 ZH is the invariant mass of the Z-Higgs system, τ = M 2 ZH /s, µ F is the factorization scale, f a (x, µ 2 F ), the parton density of the colliding proton for the parton of type a, (a = g, q,q) andσ ab is the partonic cross section for the subprocess ab → ZH +X at the partonic center-of-mass energyŝ = x 1 x 2 s. The partonic cross section can be written in terms of the Born (i.e. LO) partonic cross sectionσ (0) as: where, up to NLO terms in QCD, with α s (µ R ) the strong coupling constant defined at the renormalization scale µ R . The LO contribution G The amplitude for g µ a (p 1 )g ν b (p 2 ) → Z ρ (p 3 )H(p 4 ) can be written as where G µ is the Fermi constant and a µ (p 1 ) b ν (p 2 ) ρ (p 3 ) are the polarization vectors of the gluons and the Z boson, respectively. The tensors P µνρ i are a set of orthogonal projectors whose expressions can be found in Ref. [35]. The corresponding form factors A i (ŝ,t,û, m t , m H , m Z ) are functions of the masses of the top quark 4 (m t ), Higgs (m H ) and Z (m Z ) bosons, and of the partonic Mandelstam variableŝ whereŝ +t +û = m 2 Z + m 2 H and we took all the momenta to be incoming. The A i form factors can be expanded up to NLO terms as where the LO amplitude can be written in terms of two contributions, namely the oneloop triangle diagrams (A  contributions are understood as regularized with respect to the ultraviolet (UV) and infrared (IR) singularities via the introduction of a counterterm as in Ref. [27].
The Born partonic cross section can be written aŝ where dΦ is the two-particle Lorentz-invariant phase space. The NLO terms G (1) ab in Eq. (3) include, besides the gg channel, also the 2 → 3 processes gq → ZHq,qg → ZHq and qq → ZHg. The NLO contribution to the gg channel involves the two-loop virtual corrections to gg → ZH discussed above and one-loop real corrections from gg → ZHg. As well known, the individual contributions are IR divergent while their sum is finite. In this work, at the level of cross section, the IR singularities in all channels were treated via the dipole subtraction method [37]. The outcome of this procedure for the gg channel can be summarized as follows where C A = N c (N c being the number of colors), β 0 = 11/6 C A − 1/2 N f C F T R (N f being the number of active flavors, C F = (N 2 c −1)/(2N c ) and T R = 1/2) is the one-loop β-function of the strong coupling in the SM, P gg is the LO Altarelli-Parisi splitting function and where the plus distribution is used. The first line of Eq. (12) displays the two-loop virtual contribution, with In the second line of Eq. (12), the term R gg contains the integration over the three-particle phase space of the gg → ZHg squared amplitude minus two dipole subtractions 5 which cure the IR singularities when the final-state gluon becomes soft or becomes collinear with either of the initial-state gluons. According to our normalization (see Eqs.(2, 3)) R gg is obtained dividing the latter quantity by α s (µ R )σ (0) (zŝ) z/π. The other NLO contributions to G ab , i.e. the gq → ZHq,qg → ZHq and qq → ZHg channels, require a little discussion. The LO contribution in the hadronic reaction pp → ZH + X is given by the tree-level Drell-Yan-like process qq → Z * → ZH. The NLO real contribution, O(α s ), to this channel includes the tree-level process qq → ZHg and the crossed channels. The interference between these tree-level diagrams and their oneloop corrections is usually considered as an NNLO correction, O(α 2 s ), to the Drell-Yan-like contribution and taken into account in the NNLO evaluation of pp → ZH + X. The NLO real corrections to the gg → ZH process are formally an N 3 LO contribution, O(α 3 s ), to the pp → ZH + X reaction. In this paper we identify them as the square of the one-loop diagrams containing a closed fermion loop to which a Higgs or a Z boson, or both particles, are attached in the gq → ZHq,qg → ZHq and qq → ZHg channels. It should be remarked that our definition of the NLO contribution to gg → ZH process differs from the one employed in Refs. [28,29]. In those references, diagrams in which a Z boson is emitted from an open fermion line, which we denote as Z-radiated diagrams for brevity (see subsection 2.2), were not taken into account because they were assigned to uncalculated N 3 LO corrections to the Drell-Yan-like contribution.
The contribution of the qg → ZHq channel, and similarly forqg → ZHq, can be written as: where and R qg is obtained from the integration over the three-particle phase space of the squared amplitude of the process minus one dipole subtraction, which cures the IR singularity when the final quark becomes collinear to the initial quark, normalized to α s (µ R )σ (0) (zŝ) z/π. Finally, the qq → ZHg channel is IR safe. Its contribution can be written as where R qq is given by the square of the one-loop diagrams integrated over the phase space and normalized to α s (µ R )σ (0) (zŝ) z/π. In the following subsections we discuss the method used for the calculation and the implementation of the various contributions to Eq. (3). We evaluated our results using different renormalization schemes for the top-quark mass, namely the on-shell (OS) and the modified minimal subtraction (MS) scheme. In particular, for evaluating the top mass in the MS scheme, we first convert the OS mass to m MS t (µ t = m OS t ) using the three-loop relation [38], and then run it at three-loop order [39] numerically to the indicated scale µ t .

Virtual Corrections
Concerning the two-loop form factors A , as well as the LO contribution, were evaluated in exact top-mass dependence using the results available from Ref. [35]. Instead, the two-loop box integrals associated to A (1, ) i were evaluated combining two analytic approximations corresponding to different kinematical regimes, and here we briefly recall the main features of this approach, described in detail in Ref. [36].
The box integrals depend on five scales, namely m Z , m H , m t and the kinematic variablesŝ andt, where the latter can be traded for the transverse momentum of the final-state particles, p T . In the p T expansion of Refs. [35,40] it is assumed that the scales associated to m Z , m H and to p T are small compared to the scales set byŝ and m t . Under this assumption, the box integrals are expanded in ratios of small over large scales, and the resulting simplified integrals are written as linear combinations of 52 master integrals (MI) using Integration-by-Parts (IBP) identities obtained with LiteRed [41,42]. On the other hand, in the high-energy expansion used in Ref. [30] the two-loop box integrals are first expanded in terms of small m Z and m H , then an IBP reduction is performed on the expanded integrals; the resulting MIs are further expanded in the limit m 2 t ŝ, |t| and expressed in terms of harmonic polylogarithms.
In Ref. [36] it has been shown that, in the forward kinematic regime 6 defined by |t| ≤ |û|, the results of the p T expansion are accurate in the phase-space region |t| 4m 2 t , while the HE expansion is accurate in the complementary region |t| 4m 2 t . If the convergence of both expansions is improved using Padé approximants, then the combination of the results allows to approximate the exact cross section with an accuracy of 1% or below over the whole phase space. This procedure has been implemented in a FORTRAN program, which uses the Padé-improved HE expansion for phase-space points such that |t| > 4m 2 t and |û| > 4m 2 t , and the Padé-improved p T expansion for all the remaining phase-space points. When we compute the results using the MS scheme for the top-mass renormalization, the form factors A  to be shifted by a quantity that is defined in Eqs. (14,15) in Ref. [36]. This shift is applied to the form factors before the construction of the relative Padé approximants.
We use handyG [43] for the evaluation of polylogarithms and the routine of Ref. [44] for the evaluation of the two elliptic integrals occurring in the p T expansion results. As a result, the average timing for evaluating one phase-space point using the Padé-improved HE expansion is around 0.004 s, while using the Padé-improved p T expansion the timing ranges from 0.02 s to 0.09 s. For comparison, the evaluation of a phase-space point for the virtual terms using the results of a small-mass expansion requires on average 2 s [28].

Real Corrections
For the gg → ZHg process, we adopt Recola2 [45,46] to compute the one-loop matrix element, and we cross checked the result with MadGraph5 aMC@NLO [47]. We include all diagrams with massive and massless closed quark loops. Some representative diagrams are shown in Fig. 1. We note that in Recola2, the value of the top-quark mass cannot be changed after process initialization, and hence the process needs to be reinitialized each time if a dynamical top mass is adopted (which is the case if a dynamical scale is chosen  Fig. 2b 2d) the Higgs boson is attached to a closed quark loop, but the Z boson is radiated from an open fermion line. We note that both types of diagrams can interfere with tree-level diagrams, hence produce O(α 2 s ) contributions. Such contributions were studied in detail 7 in Ref. [11] and they were considered as part of the NNLO corrections to pp → ZH. On the other hand, in this paper we compute the square of those diagrams, corresponding to O(α 3 s ) contributions that we consider as NLO corrections to gg → ZH.

Results
In this section, we present our numerical results for a center-of-mass energy √ s = 13 TeV.
We adopt the following input parameters: m OS t = 172.5 GeV, m W = 80.385 GeV, m Z = 91.1876 GeV, m H = 125 GeV, G µ = 1.1663787 × 10 −5 GeV −2 . We adopt the NNPDF31 nnlo as 0118 [48] parton distribution functions in a five flavour scheme. 7 They belong to the classes R I and R II for the top-mediated terms considered in Ref. [11].

Inclusive Cross Section
In Table 1, we show the total cross section at LO and NLO adopting different top-quarkmass renormalization schemes, i.e. OS and MS with different scale choices. We fix the central value of the renormalization and factorization scales to be µ C = M ZH /2. The scale uncertainty is obtained from the envelope of a 7-point variation of the central scale according to (µ R /µ C , µ F /µ C ) = (1, 1), (1, 1 2 ), (1, 2), ( 1 2 , 1 2 ), ( 1 2 , 1), (2, 1), (2, 2). We find that the NLO corrections are large for each choice of the top-mass renormalization scheme, with an approximate K-factor, K = σ N LO /σ LO , of around 2. Moreover, the relative size of the scale uncertainties is essentially the same regardless of the top-mass renormalization scheme. We note that going from LO to NLO the relative size of the scale uncertainties is reduced by a factor of about 2/3. The OS scheme leads to the largest value of the total cross section both at LO and NLO, while in the MS scheme for µ t = M ZH the smallest cross section value is obtained. At LO, the difference between these two schemes amounts to about 23%, while it decreases to 13% at NLO.
We notice that our OS results are about 20% larger at LO and 14% larger at NLO than those of Ref. [29] (see Table 1 therein). This discrepancy is mainly due to the different choice for µ C , and it is related only in a minor way to the additional diagrams included in our calculation and to the different input parameters adopted. To verify this, we have computed our results including the same diagrams and adopting the same input parameters as in Ref. [28] (which is in accordance with Ref. [29]) and we have found an agreement at the level of the Monte Carlo error. Furthermore, when we consider the relative importance of the scale uncertainties, we observe very similar results to Ref. [29].

Differential Distributions
In Fig. 3, we plot the M ZH distribution in both the OS scheme (3(a)) and the MS scheme with µ t = M ZH /2 (3(b)) in the region M ZH ∈ [200, 800] GeV. In both schemes, the K-factor is about 3 in the ZH threshold region, then it decreases as M ZH increases. In the top-pair threshold region (M ZH ∼ 2 m t ), the OS scheme gives a peak with K-factor slightly above 2, while the MS scheme shows a small dip followed by a peak instead. Increasing M ZH to about 800 GeV, the K-factor in the OS scheme decreases to about 1.5, while it remains If M ZH is increased to very large values, as shown in Fig. 4, we observe that the Kfactor starts to increase rapidly. At M ZH = 2.5 TeV, the value of the K-factor can reach ∼ 6 in the OS scheme, and ∼ 10 in the MS scheme. Such behaviour is due to the inclusion of diagrams where the Z boson is radiated from an open quark line, as in Figs. 2b and 2d. For comparison, in Fig. 4 we show also the NLO cross section when these contributions are excluded. Indeed, one finds that in the latter case the K-factor remains rather flat at high M ZH in both schemes.
To further assess the contribution of these Z-radiated diagrams, in Fig. 5(a) we show various pieces of the O(α 3 s ) corrections to the pp → ZH + X cross section: we compare the isolated contribution of the Z-radiated diagrams for the qg and qq initial states (green and yellow lines, respectively) to the contribution from each of the two partonic channels including all the relevant diagrams (orange and red lines for qg and qq, respectively). We can see that in the qg channel when M ZH > 1 TeV, the dominant contribution comes from the square of the Z-radiated diagrams. We ascribe this feature to the contribution from logarithmic terms of EW origin, of the form log 2 (m 2 Z /M 2 ZH ), which become large when the typical scale of the process (M ZH ) and the EW scale (represented by m Z ) are very different, see e.g. Ref. [49]. On the other hand, although the qq channel includes diagrams where the Z boson is radiated from the external quark lines, its size remains negligible with respect to the total cross section, though also there we observe that the Z-radiated diagrams are dominating the respective initial state at high M ZH . This suppression can be mainly attributed to the reduced partonic luminosity with respect to the qg channel. For comparison, in Fig. 5(a) we also report the size of the Drell-Yan type contribution at NNLO (black line), which we obtained using vh@nnlo [18,25] with MCFM [50][51][52]. In the lower panel of Fig. 5(a) we plot the ratio of the O(α 3 s ) corrections computed by us with respect to the NNLO Drell-Yan contribution. We can see that despite being O(α 3 s ), the relative importance of the Z-radiated contribution can reach 2% when M ZH ∼ 2 TeV.
In Fig. 5(b) we compare our results for gg → ZH at LO (green line) and NLO (blue line) with the Drell-Yan type contribution (black line). In the upper panel we show the size of the differential cross section for the various channels, while in the lower panel the ratio of the gluon-fusion with respect to the NNLO Drell-Yan contribution is displayed. We can see that the gluon-fusion contribution peaks around the top-pair threshold, which increases its relative size over the Drell-Yan contribution by about 25% at LO, and about 45% at NLO. The relative size of the gluon-fusion contribution decreases above the top-pair threshold as M ZH increases, and at NLO becomes dominated by the Z-radiated terms for very large values of M ZH . In particular, at 2 TeV the latter constitute more than half of the gluon-fusion contribution. only Z-rad, qg → qZH real, qg → qZH only Z-rad, qq → gZH real, qq → gZH

Change of Renormalization Scheme
To assess the impact of the top-quark-mass renormalization scheme at differential level, we show in Fig. 6 the differential cross section at LO (6(a)) and NLO (6(b)) in various top-mass schemes. The lower panel shows the ratio of the MS results with respect to the OS scheme.
Before discussing the results, we note that the value of the top mass can affect the shape of the distributions in two ways: first, in the high-energy regime, it directly controls the overall size of the LO and NLO contributions via the proportionality of the amplitude to the top-Yukawa coupling 8 ; second, in the low M ZH region the value of m t shifts the position of the peak associated to the top-pair threshold, M ZH = 2m t . The former effect has been already observed by the authors of Ref. [29], and our results are in agreement when the same scale choices (µ t = M ZH and µ t = m MS t ) and the same invariant-mass range   (M ZH > 400 GeV) are considered. On the other hand, since we can use our analytical results to investigate the low M ZH region, in Fig. 6 we are also able to observe the effect of the peak shift on the invariant-mass distribution. In particular, we can see that the ratio (orange line) between µ t = M ZH and the OS scheme is about 2 in the bin M ZH ∈ [300, 325] GeV at LO, and decreases to about 1.75 at NLO. Instead, for µ t = M ZH /4 (green line) this effect is rather small. Indeed, since renormalization-group evolution predicts that m MS t (µ t ) decreases monotonically as µ t increases, we expect larger deviations from the OS results when larger values of µ t are chosen, see also Ref. [54]. In the region M ZH > 350 GeV, the OS scheme provides the largest cross section, while the cross section for the choice µ t = M ZH is the smallest. Across all regions, we can see that going from LO to NLO reduces the difference among different top-mass schemes.
To quantify the size of top-mass renormalization uncertainties, we follow the procedure adopted in Refs. [54,55]. In this approach, after a binning of the M ZH range is chosen, the OS scheme result is taken as the central value, while for each bin the uncertainty is obtained using the minimal and maximal values from a set of results in different renormalization schemes. In particular, for the top mass we considered the set m OS  results under different choices of the bin size. Due to the way in which the overall uncertainty is constructed, the latter becomes of course bigger the smaller the bin size is chosen. This is though mainly due to the top-pair threshold region, where the location of the bin can lead to a bigger or smaller top-mass renormalization uncertainty in dependence on how well the structure at the top-pair threshold 9 is resolved by the binning, while for large invariant masses the OS scheme is always the largest and the scheme using m MS t (M ZH ) is the smallest. In the very-low invariant-mass region (i.e. M ZH 275 GeV) the scheme dependence is small, as this region can be well described by a large mass expansion. We hence expect that if the description of the top-pair threshold region was improved by a Coulomb resummation, the dependence on the binning of the uncertainty would become smaller.
We observe that the uncertainty nearly halves when going from LO to NLO for choosing a bin size ≥ 100 GeV. This is similar to what has been observed for gg → HH in Ref. [55]. When smaller bin sizes are adopted, the uncertainty does not half but still shows a sizable reduction.

Conclusions
In this paper we have presented the evaluation of the QCD corrections at NLO for the gginitiated channel to ZH associated production. The computation of the virtual corrections are based on the approach of Ref. [36], while for this work we have computed the real corrections and included them in a fast and flexible code. For the inclusive cross section, we have found that the NLO QCD corrections are of the same size as the LO contribution, and they increase the cross section by a factor of about 2. Our findings are in accordance with independent results in the literature [28,29]. We have also studied the invariant-mass distribution of the gg → ZH channel at NLO, observing that the perturbative K-factor is not flat across the M ZH range, specifically in the region of threshold ZH production and in the very-high-energy tail (M ZH > 1 TeV). Furthermore, we have shown that in this latter region the contribution from real-emission diagrams in which the Z boson is emitted from an open quark line is the dominant one, and it causes the K-factor to rise up to 10 in the MS scheme. We expect that in W H production an analogous class of W -radiated diagrams will give a similar contribution in the high-energy regime. While the focus of our analysis is on LHC phenomenology, this feature could be of interest for studies at future colliders.
The gg-initiated channel is responsible for the larger theoretical uncertainties in the prediction for ZH production, compared to its W H counterpart. We have shown that the inclusion of the NLO corrections bring only a mild reduction of the scale uncertainties, about a factor 2/3, suggesting that more accurate calculations are necessary in order to describe the gg-initiated channel at a level that is adequate for experimental studies in the future. Finally, the implementation used for our results allowed us to study for the first time the impact of the uncertainty due to the renormalization scheme for the top quark mass over the whole invariant-mass range: we have found that different choices for the top mass scheme can lead to substantially different results, and we suggest that this uncertainty should be included in refined theoretical predictions.