Next-to-leading order corrections for $gg \to ZH$ with top quark mass dependence

In this Letter, we present for the first time a calculation of the complete next-to-leading order corrections to the $gg \to ZH$ process. We use the method of small mass expansion to tackle the most challenging two-loop virtual amplitude, in which the top quark mass dependence is retained throughout the calculations. We show that our method provides reliable numeric results in all kinematic regions, and present phenomenological predictions for the total and differential cross sections at the Large Hadron Collider and its future upgrades. Our results are necessary ingredients towards reducing the theoretical uncertainties of the $pp \to ZH$ cross sections down to the percent-level, and provide important theoretical inputs for future precision experimental collider programs.


Introduction
One of the top priorities of the Large Hadron Collider (LHC) and the High Luminosity LHC (HL-LHC) is to accurately measure the properties of the Higgs boson including its various couplings. In particular, the Yukawa couplings of the Higgs boson to fermions are key to understand the origin of fermion masses as well as to verify the family universality of different fermion generations. At the LHC and the HL-LHC, the main production channel gg → H followed by hadronic decays of the Higgs boson is overwhelmed by large hadronic activities arising from strong interactions. Therefore, to measure the Yukawa couplings of quarks lighter than the top quark, the pp → ZH production channel is the most promising one, where the leptonic decay of the Z boson can be efficiently detected and leads to a significant reduction of the hadronic background [1,2].
Besides, the pp → ZH process also provides an independent handle of the HZZ coupling with respect to the H → ZZ * decay process [3,4]. Therefore, it is of crucial importance to have precise knowledge of its differential cross sections both from the theoretical and the experimental points of view.
At the partonic level, the pp → ZH process proceeds through the subprocess qq → Z * → ZH at the leading order (LO) in standard model (SM) couplings. Since the quantum effects due to strong interactions are often quite large for such processes, it is necessary to calculate higher order corrections in the strong coupling α s within perturbative quantum chromodynamics (QCD), in order to improve the accuracy for the theoretical predictions. The next-to-leading order (NLO) and next-to-next-to-leading order (NNLO) QCD corrections were calculated in [5][6][7] and are implemented in the program package vh@nnlo [8,9]. The NLO electroweak (EW) corrections are given by [10,11]. Combined QCD+EW corrections at NLO are also available [12][13][14]. Despite of these progresses, the theoretical precision is still not enough to match the potential of future LHC and HL-LHC runs [15].
There is one special class of higher order corrections coming from the partonic subprocess gg → ZH mediated by a closed top quark loop. Such contributions enter at order α 2 s in perturbation theory. However, this formal suppression by the coupling constant is compensated partly by the high luminosity of gluons at colliders such as the LHC and the HL-LHC, and partly by the large Yukawa coupling of the top quark. Moreover, the differential cross sections for this subprocess receive additional enhancement by the 2m t threshold effects in the boosted regime, e.g., when the Higgs transverse momentum p T ≥ 150 GeV [16,17]. For the total cross section at the 13 TeV LHC, the one-loop gg → ZH contribution is roughly 50% in size compared to the NLO corrections for the qq → ZH channel, and is much larger than the NNLO corrections for the qq channel. Therefore, it is reasonable to expect that the formally α 3 s contributions in the gg → ZH category are not small compared to the formally α 2 s terms in the qq channel. In fact, the estimate in the heavy top limit [18,19] suggests that, the gg → ZH contribution at order α 3 s is similar in size to that at order α 2 s . In the case of Higgs boson pair production, it was realized that the heavy top limit is not quite reliable, and the finite top mass effects are significant [19]. As a result, the unknown size of the exact gg → ZH contribution at order α 3 s has become a major uncertainty of the theoretical predictions for the pp → ZH cross sections [15].
Last but not least, the loop-induced gg → ZH subprocess is interesting on its own right. It is sensitive to possible new heavy particles running in the loop. From the low energy point of view, it has a different dependence on the HZZ anomalous coupling than the tree-mediated qq channel. It is also sensitive to other anomalous couplings such as Htt and Ztt. This channel therefore provides unique information on the possible new physics beyond the SM.
The above considerations make it highly interesting to calculate the gg → ZH subprocess to the next order in perturbation theory. With a slight abuse of notation (as adopted in the literature), we will refer to the one-loop amplitudesquared contributions as the "LO" cross sections for this process. These have been computed in [20,21]. The "NLO" corrections to this process then consist of two parts: the interference between the two-loop and one-loop amplitudes (the virtual contributions), and the squared-amplitude with one loop and one extra parton radiation (the real contributions). Both contributions are infrared (IR) divergent, while their sum leads to a finite prediction. The bottleneck in such a calculation lies in the two-loop virtual amplitude, which involves multiple physical scales including three masses (m Z , m H and m t ) and two Mandelstam variables.
The presence of multiple physical scales makes the computation of the twoloop amplitude rather challenging. Very recently, a purely numeric study based on sector decomposition [22][23][24][25] has appeared [26]. Approximations in certain kinematic regions have also been worked out, including the high energy expansion [27] and the small transverse momentum expansion [28]. In [29], some of the authors of this work proposed a novel approach to calculate loop integrals with a heavy top quark loop and lighter external particles such as the Higgs boson and weak gauge bosons, which is valid in all phenomenologically relevant kinematic regions. In this work, we apply this method to calculate the NLO virtual corrections to the process gg → ZH. We show that our method gives accurate numeric results for the two-loop amplitudes in the entire physical phase-space. Combining with the real corrections, we for the first time provide the complete NLO predictions for the total and differential cross sections of the gg → ZH channel. We also give the total cross sections for the full pp → ZH process by adding the qq contributions.

Setup and analytic calculations
We consider the partonic process g α a (p 1 ) + g β b (p 2 ) → Z µ (p 3 ) + H(p 4 ), where a and b are color indices, while α, β and µ are Lorentz indices. The amplitude can be written as where the tensor structures A αβµ i are given in [21], whose coefficients F i are functions of the masses m t , m Z , m H and the Mandelstam variablesŝ = ( These variables satisfŷ s+t 1 +û 1 = 0. For later convenience, we also definet =t 1 +m 2 Z andû =û 1 +m 2 H . To calculate the squared-amplitude, we multiply Eq. (1) by its complex conjugate, and sum over the polarization states of the gluon and the Z boson. For simplicity, we define the coefficientsĈ i ≡ 7 j=1 M ij F j , where the matrix elements M ij are given by The coefficientsĈ i can be expanded in terms of the strong coupling constant are the NLO virtual corrections, which involve complicated two-loop Feynman integrals with massive external legs. In the following, we describe the calculation ofĈ (1) i using the method of small mass expansion. We generate the relevant two-loop virtual diagrams and the corresponding amplitudes using FeynArts [30]. The resulting expressions are manipulated with FeynCalc [31] and FORM [32]. The two-loop diagrams can be classified into several categories: diagrams consisting of two one-loop sub-diagrams, two-loop triangle diagrams (involving both top quark and bottom quark loops attached to an offshell Z boson), and two-loop box diagrams. While the first two categories can be easily calculated [19,27,[33][34][35], the two-loop box diagrams are rather challenging. We denote the corresponding contributions to theĈ (1) i coefficients asĈ (1) i,box , and calculate them using the small mass expansion. Namely, we expandĈ (1) i,box as a power series in m 2 Z and m 2 H : It should be noted that in contrast to the case of Higgs boson pair production [29,36], there is an extra factor of 1/m 2 Z due to the polarization sum of the Z boson, as is evident in Eq. (2). For the power counting we regard m 2 Z ∼ m 2 H , and denote the nth term (n ≥ 0) in the expansion as order m 2(n−1) . We will also use the notation O(m 2(n−1) ) to denote the partial sum of the series up to the nth term. In practice, we have performed the expansion up to n = 3, corresponding to O(m 4 ).
Acting on the amplitudes, the partial derivatives ∂ m 2 Z and ∂ m 2 H can be written as partial derivatives with respect to the external momenta: H → 0 after taking the derivatives, the coefficientsĈ (1) i,box can be written as linear combinations of scalar loop integrals with massless external legs. These integrals have already been computed in [29,[36][37][38].
The loop integrals contain both ultraviolet (UV) and infrared (IR) divergences which are regularized in dimensional regularization. The UV divergences are removed via renormalization. We adopt the on-shell scheme for m t and the five-flavor MS scheme for α s . In addition, to handle γ 5 in the amplitudes, we apply the Larin scheme [39] which requires a finite renormalization constant The IR divergences will be eventually cancelled by the real corrections and the renormalization of the parton distribution functions (PDFs). We apply the Catani-Seymour dipole subtraction method [40] to subtract the IR divergences from the virtual corrections. After that, we define the UV-and IR-finite coef-ficientsĈ i,ren are the renormalized coefficients and the IR subtraction operator I is defined as where β 0 = 11C A /3 − 4T F N l /3 with T F = 1/2 and N l = 5, andK g is the finite part and is given byK The contributions of the finite coefficientsĈ (1) i,fin to the squared-amplitude can be summarized into the so-called finite part of the NLO virtual corrections: In the small mass expansion, the computation of V fin is rather efficient, which only takes 2 seconds on average for one phase-space point on a workstation with an Intel Xeon W-2155 CPU. Most of the evaluation time is spent on the computation of scalar master integrals. To further accelerate the phase-space integration (which requires to repeatedly evaluate the amplitude), we have generated a very large grid for the master integrals with 63175 points on the twodimensional plane (−4m 2 t /ŝ, −4m 2 t /t 1 ). Note that this grid can be reused for different masses (including m t , m H and m Z ), as well as for different couplings (HZZ, Htt and Ztt). Using the grid bypasses the most time-consuming part of the numeric evaluation. As a result, computing the amplitude with the grid only takes 0.003 seconds on average per phase-space point.
We now need to add the IR-subtracted real corrections to obtain physical cross sections. We consider squared amplitudes of four partonic subprocesses gg → ZH + g, gq → ZH + q, gq → ZH +q and qq → ZH + g, and perform the dipole subtraction to remove the IR divergences in the phase-space integration. In the above four subprocesses, we only include the diagrams with a closed quark loop. The IR limit of these diagrams has the same topology as that of the two-loop virtual diagrams. This selection also allows us to compare our results with those in the heavy top limit [18]. We compute the amplitudes using the package Gosam [41,42], and integrate over the phase-space using the Vegas algorithm implemented in the Cuba library [43]. To avoid numeric instabilities, we require the transverse momentum of the extra radiation to be larger than δ √ŝ . We have varied δ between 0.05 and 0.0005 and find that the result is stable within the integration precision of 0.3%.

Numeric results
In this section, we present our numeric results for the total cross sections and the M ZH distributions, where M ZH = (p 3 + p 4 ) 2 is the invariant mass  of the Z boson and the Higgs boson. For the input parameters, we take G F = 1.16637 × 10 −5 GeV −2 , m Z = 91.1876 GeV, m H = 125 GeV, m t = 172.5 GeV [44], and use NNPDF3.1 NNLO PDFs [45] with α s (m Z ) = 0.118. We neglect the mass of the bottom quark appearing in the loop. According to [46], we set the default values for the renormalization scale µ r and the factorization scale µ f to be µ def = M ZH , while the scale uncertainties are estimated by varying the two scales simultaneously up and down by a factor of three.
The validity of the small mass expansion has already been demonstrated in the case of Higgs boson pair production [36]. To be more careful, we have performed a similar analysis for ZH production. First of all, we have compared our results for the finite part of the virtual corrections with the results of Ref. [26] obtained using the sector decomposition program pySecDec. It should be noted that the conventions for the finite part are slightly different between Ref. [26] and our work, due to the different choices of the IR subtraction operator I( ) 1 .  For comparison, we have performed the necessary transformation to arrive at their convention (dubbed V fin later). The results of V fin at six phase-space points are listed in Table 1. We find that our results at O(m 4 ) agree quite well with the pySecDec results. In particular, for the first four points where the pySecDec results are accurate enough, the relative errors of our results are much smaller than 0.1%. We emphasize that our results can still be systematically improved by incorporating higher order terms in the small mass expansion. The last two points correspond to the high energy region far above the 2m t threshold, where the pySecDec results show large numeric uncertainties. On the other hand, the small mass expansion is expected to converge fast in the high energy region, and our results should be reliable here. The high energy expansion of Ref. [27] also works in this region, and it will be interesting to compare our results with theirs as a crosscheck.
In Table 1, we also list the O(m 0 ) and O(m 2 ) results of our calculation, which shows the excellent convergence of the small mass expansion. To examine the convergence for a broader range of phase-space points, we show V fin as a function of √ŝ for several representative values of p T in Fig. 1. The blue and red marks represent the O(m 0 ) and O(m 2 ) results, respectively; while the black lines correspond to the O(m 4 ) ones. It can be seen that the O(m 2 ) and O(m 4 ) results almost overlap with each other completely, which demonstrates the reliability of the expansion in the entire phase-space. We expect that the terms at O(m 6 ) are irrelevant for phenomenological applications.
We now combine the finite part of the virtual corrections with the IRsubtracted real corrections, and present our predictions for the total and differential cross sections. We first consider the LHC with a center-of-mass energy of √ s = 13 TeV. We use the program package vh@nnlo [8,9] to calculate the contributions from the qq channel (including QCD and EW corrections). This program also gives the gg → ZH contributions up to the NLO in the heavy top limit, which we use as a reference to compare our results with. The various results for three values of µ r = µ f are listed in Table 2. As expected, the NLO corrections lead to significant enhancement (about 100%) to the gg → ZH cross section. Combining our results with the qq contributions, we arrive at the state-of-the-art fixed-order prediction for the pp → ZH total cross section at the 13 TeV LHC: In the last two columns of Table 2, we show for comparison the results in the heavy top limit given by vh@nnlo. We find that the situation is quite different from the Higgs boson pair production: the finite top mass effects are much milder, which reduces the NLO cross sections in the gg channel only by about 4%. This accidental fact makes it possible that by calculating the O(α 4 s ) contributions in the heavy top limit, one could reduce the residue theoretical uncertainty of the total cross section down to the percent-level. We now turn to the differential cross sections. It is well-known that the heavy top limit is not valid above the 2m t threshold. On the other hand, the small mass expansion provides reliable results for differential cross sections in the entire phase-space. As an example, we show in Fig. 2 the LO and NLO differential cross sections in the gg → ZH channel with respect to the invariant mass M ZH of the Z boson and the Higgs boson at the 13 TeV LHC. The upper plot employs a logarithmic scale for the vertical axis to access the distributions in the broad range 200 GeV M ZH 2500 GeV, while the lower plot shows the ratios to the central values of the LO differential cross sections. It is clear that the sizes of the corrections are kinematics-dependent, and it is not sufficient to use a uniform K-factor to rescale the LO differential cross sections. The NLO corrections are rather large across the whole range, especially around the peak and to the far tail. The significant corrections in the tail region have important implications for new physics searches, since new phenomena are usually most evident in the high energy regime. The boosted region is also relevant when using jet substructure techniques to measure the hadronic decays of the Higgs boson [47].   The corrections in the peak region, on the other hand, are the most important to the total cross section. To see the peak region more clearly, we show the M ZH distributions in a narrower range in Fig. 3, with a linear vertical axis. It can be seen that the total cross section receives its most contributions from regions around the 2m t threshold, where the NLO corrections are significant. The ratio plot also shows a small kink at the 2m t threshold, which comes from the Coulomb-type enhancement in that region entering at NLO. Finally, we envision a possible future high-energy upgrade of the LHC (HE-LHC) operating at a center-of-mass energy of 27 TeV. In Table 3 we list the results for the total cross section at 27 TeV. Again, the NLO corrections are significant, with the top quark mass effect slightly larger than that in the 13 TeV case. The differential cross sections can also be easily computed, which we leave for future investigations.

Conclusion
In this work, we present for the first time a calculation of the complete NLO corrections to the gg → ZH process. We use the method of small mass expansion to tackle the most challenging two-loop virtual amplitude, in which the top quark mass dependence is retained throughout the calculations. We compare our results of the two-loop amplitude with the purely numeric results from pySecDec. We find that at phase-space points where pySecDec is precise enough, the relative deviations of our results are much smaller than 0.1%. We have also demonstrated the excellent convergence of the small mass expansion in the entire phase-space, which makes us confident that the expansion up to O(m 4 ) is sufficient for phenomenological applications.
We employ the dipole subtraction method to combine the virtual corrections with the real radiation contributions, and find that the IR divergences all cancel out. This allows us to give numeric predictions for the total and differential cross sections at the NLO. We add the contributions from the qq channel to obtain the state-of-the-art fixed-order predictions for the total cross sections, which amount to σ pp→ZH = 882.9 +3.5% −2.5% fb at the 13 TeV LHC, and σ pp→ZH = 2.555 +4.0% −2.7% pb at the 27 TeV HE-LHC. We further present our results for a representative differential cross section: the invariant mass distribution of the Z boson and the Higgs boson. We demonstrate that our method can provide reliable predictions for the differential cross sections from the low energy region all the way up to the highly-boosted regime. Our results are necessary ingredients towards reducing the theoretical uncertainties of the pp → ZH cross sections down to the percent-level, and provide important theoretical inputs for future precision experimental programs at the LHC and the HL-LHC. The results for other phenomenologically interesting distributions at the 13 TeV LHC, and the distributions at the 14 TeV LHC/HL-LHC and the 27 TeV HE-LHC will be presented in a forthcoming article.