Virtual corrections to $gg\to ZH$ via a transverse momentum expansion

We compute the next-to-leading virtual QCD corrections to the partonic cross section of the production of a Higgs boson in association with a $Z$ boson in gluon fusion. The calculation is based on the recently introduced method of evaluating the amplitude via an expansion in terms of a small transverse momentum. We generalize the method to the case of different masses in the final state and of a process not symmetric in the forward-backward direction exchange. Our analytic approach gives a very good approximation (better than percent) of the partonic cross section in the center of mass energy region up to $\sim 750 \,\textrm{GeV}$, where at the LHC $\sim 98\%$ of the total hadronic cross section is concentrated.


Introduction
The study of the characteristics of the Higgs boson is one of the primary tasks of the LHC program: the forthcoming Run3 and the High-Luminosity phase will increase the accuracy in the measurement of Higgs production cross sections and decay rates, allowing for a more stringent test of the Standard Model (SM) predictions. One of the main production processes being investigated is the so called Higgs-strahlung process pp → V H, in which a single Higgs boson is emitted together with a weak vector boson (V = Z, W ). The leptonic decays of the weak boson can be exploited as a trigger for measurements of elusive Higgs decays. In particular, the decay H → bb has been observed for the first time by ATLAS and CMS, using an analysis focused precisely on the associated production category [1,2].
In this paper, we are interested in the associated production of a Higgs and a Z boson. The theoretical predictions for pp → ZH are accurate at next-to-next-to-leading-order (NNLO) in QCD, and at next-to-leadingorder (NLO) in the EW interactions [3]. The leading and next-to-leading contributions are connected to the qq-initiated channel, allowing to interpret pp → ZH mainly as a Drell-Yan process [4,5].
The gluon-initiated channel gg → ZH arises for the first time at NNLO in QCD. It is an O(α 2 S ) correction, but the contribution from this process to the hadronic cross section is non-negligible because of the large gluon luminosity at the LHC. It has been shown that the relevance of gg → ZH is even more enhanced in the boosted kinematic regime, to the point of being comparable to the quark-initiated contribution near the tt threshold [6]. The factorization-and renormalization-scale uncertainties related to the gluon-induced process also affect significantly the uncertainty on the total pp → ZH cross section. This issue is specific to the ZH final state, since the gluon-induced channel is absent in pp → W H. The knowledge of the NLO corrections to gg → ZH would reduce the scale uncertainties, facilitating precision studies in the next runs of the LHC. The gg → ZH contribution is relevant also for New Physics (NP) studies, since it is sensible to both sign and magnitude of the top Yukawa coupling, dipole operators [7] and can receive additional contributions from new particles [8].
An improved knowledge of the SM prediction for the gluon-induced contribution is therefore very important both for precision measurements of ZH production within the SM and for testing NP in this channel. The leading order (LO) contribution to the gg → ZH amplitude, given by one-loop diagrams, was computed exactly in refs. [9,10]. At the NLO the virtual correction part contains two-loop multi-scale integrals that constitute, at present, an obstacle to an exact evaluation of the NLO contribution. Specifically, the corrections due to the two-loop box diagrams are still not known analytically. A first computation of the NLO terms was obtained in ref. [11] using an asymptotic expansion in the limit m t → ∞ and m b = 0, and pointed to a K-factor of about 100% with respect to the LO contribution. Soft gluon resummation has been performed in ref. [12] including next-to-leading logarithmic terms, and the result has been matched to the fixed NLO computation of ref. [11]. Finite top-quark-mass effects to gg → ZH have been investigated in ref. [13] using a combination of large-m t expansion (LME) and Padé approximants. In addition, a data-driven method to extract the non-Drell-Yan part of pp → ZH, which is dominated by the gluon-induced contribution, has been proposed in ref. [14], exploiting the known relation between W H and ZH associated production when only the Drell-Yan component of the two processes is considered. A qualitative study focusing on patterns in the differential distribution has been conducted in ref. [15], where 2 → 2 and 2 → 3 LO matrix elements were merged and matched to improve the description of the kinematics.
Very recently, a new analytic computation of the NLO virtual contribution based on a high-energy expansion of the amplitude, supported by Padé approximants, and on an improved LME, has been carried out [16]. The results are in agreement with a new exact numerical study [17], in the energy regions where the expansions are legitimate. Nonetheless, an improvement on the analytic calculation is still desirable, since the heavy-top and the highenergy expansions do not cover well the region 350 GeV √ŝ 750 GeV, where √ŝ is the partonic center of mass energy. It should be remarked that this region provides a significant part of the hadronic cross section at the LHC, about 68%.
In this paper, we present an analytic calculation of the virtual NLO QCD corrections to the gg → ZH process that covers the region √ŝ 750 GeV, which contributes about 98% to the hadronic cross section. The most difficult parts, i.e. the two-loop box diagrams, are computed in terms of a forward kinematics [18] via an expansion in the Z (or Higgs) transverse momentum, p T , while the rest of the virtual corrections is computed exactly. We remark that our calculation is complementary to the results of ref. [16], which covers the region of large transverse momentum of the Z. Furthermore, the merging of the two analyses allows an analytic evaluation of the NLO virtual corrections in gg → ZH in the entire phase space.
The paper is structured as follows: in the next section we introduce our notation and the definitions of the form factors in terms of which we express the amplitude. In section 3, we present the expansion of the amplitude in terms of the Z transverse momentum. Section 4 is devoted to a discussion of the expected range of validity of the evaluation of the amplitude via a p T -expansion, by comparing the exact result for the LO cross section with the p T -approximated one. In section 5 we present an outline of our NLO computation, while the next section contains our NLO results. Finally we present our conclusions. The paper is complemented by two appendices. In appendix A, we report the explicit expressions for the orthogonal projectors we employ in the calculation. We present also the relation between our form factors and the ones used in ref. [16]. In appendix B, we report the exact results for the triangle and the reducible double-triangle contributions.

Definitions
In this section we introduce our definitions for the calculation of the NLO QCD corrections to the associated production of a Higgs and a Z boson from gluon fusion.
The amplitude g µ a (p 1 )g ν b (p 2 ) → Z ρ (p 3 )H(p 4 ) can be written as where G F is the Fermi constant, α S (µ R ) is the strong coupling constant defined at a scale µ R and a µ ( 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 and the Born partonic cross section can be written aŝ . The Feynman diagrams that contribute to the gg → ZH amplitude up to NLO can be separated into triangle, box and double-triangle contributions, the last type appearing for the first time at the NLO level. Examples of LO (NLO) triangle and box categories are shown in fig.1 (a) -(c) ((d) -(f )). Due to the presence of a γ 5 in the axial coupling of the Z boson to the fermions in the loop, the projectors P µνρ i are proportional to the Levi-Civita total antisymmetric tensor αβγδ (see appendix A), whose treatment in dimensional regularization is, as well known, delicate and will be discussed in section 5.
In our calculation we treat all the quarks but the top as massless. As a consequence, the contribution to the amplitude of the first two generations vanishes. Concerning the third generation, the contribution of the bottom is present in the triangle diagrams with the exchange of a Z boson ( fig.1(b), (e)) and in the double-triangle diagrams ( fig.1(g)). A nice observation in ref. [11]  allows to compute easily the full (top+bottom) triangle contribution. As noticed in that reference, the triangle contribution with a Z exchange contains a ggZ * subamplitude which in the Landau gauge can be related to the decay of a massive vector boson with mass √ŝ into two massless ones, a process that is forbidden by the Landau-Yang theorem [19,20]. As a consequence, the full triangle contribution can be obtained from the top triangle diagrams with the exchange of the unphysical scalar G 0 , with the propagator of the G 0 evaluated in the Landau gauge. This part of the top triangle diagrams can be obtained from the decay amplitude of a pseudoscalar boson into two gluons which is known in the literature in the full mass dependence up to NLO terms [21,22].
Given the above observation, our calculation of the NLO corrections to the gg → ZH amplitude focuses on the analytic evaluation of the doubletriangle ( fig.1(g)) and two-loop box contributions ( fig.1(f )). The former contribution is evaluated exactly. The latter is evaluated via two different expansions: i) via a LME, following ref. [23], up to and including O(1/m 6 t ) terms, which is expected to work below the 2 m t threshold; ii) via an expansion in terms of the Z transverse momentum, following ref. [18], whose details are presented in the next section.

Expansion in the transverse momentum
The transverse momentum of the Z boson can be written in terms of the Mandelstam variables as From eq.(6), together with the relation between the Mandelstam variables, one finds where ∆ m = (m 2 H − m 2 Z )/2. Eq. (7) implies p 2 T /ŝ < 1 that, together with the kinematical constraints m 2 H /ŝ < 1 and m 2 Z /ŝ < 1, allows the expansion of the amplitude in terms of these three ratios.
A direct expansion in p T is not possible at amplitude level, since p T itself does not appear in the amplitudes. However, as we argued in ref. [18], the expansion in p 2 T /ŝ 1 is equivalent to an expansion in terms of the ratio of the reduced Mandelstam variables t /s 1 or u /s 1, depending whether we are considering the process to be in a forward or backward kinematics. The s , t and u variables are defined as and satisfy The cross section of a 2 → 2 process can always be expanded into a forward and backward contribution. Looking at the dependence of σ upon t , u we can write and t m is the value of t at which t = u = (−s +∆ m )/2. The two terms in the second line of eq.(10) represent the expansion in the forward and backward kinematics, respectively. If the amplitude is symmetric under t ↔ u exchange then (11) so that the expansion in the forward kinematics actually covers the entire phase space.
In the case of gg → ZH the process itself is not symmetric under the t ↔ u exchange. However, as can be seen from the explicit expressions of the projectors in appendix A, it can be written as a sum of symmetric and antisymmetric form factors. To perform only the expansion in the forward kinematics one can proceed in the following way. On the symmetric form factors the expansion can be directly performed. For the antisymmetric ones, it is sufficient first to extract the overall antisymmetric factor (t −û) just by multiplying the form factor by 1/(t −û), written as 1/(2s − 4t − 2∆ m ), then perform the expansion in the forward kinematics and finally multiply back by (t −û).
As discussed in ref. [18], to implement the p T -expansion at the level of Feynman diagrams it is convenient to introduce the vector r µ = p µ 1 + p µ 3 , which satisfies and therefore can be also written as where From eq.(6) one obtains that implies that the expansion in small p T (the minus sign case in eq.(15)) can be realized at the level of Feynman diagrams, by expanding the propagators in terms of the vector r µ around r µ ∼ 0 or, equivalently, p µ 3 ∼ −p µ 1 , see eq. (13).
The outcome of the evaluation of the gg → ZH amplitude via a p Texpansion is expressed in terms of a series of Master Integrals (MIs) that are functions ofŝ and m 2 t only, and whose coefficients can be organized in terms of powers of ratios of small over large parameters where p 2 T , m 2 H and m 2 Z are identified as the small parameters while m 2 t andŝ as the large ones. Thus, the range of validity of the expansion depends on the condition that p 2 T can be treated as a "small parameter" with respect to m 2 t because all the other ratios, small over large, are always smaller than 1.

LO Comparison
In order to investigate the range of validity of the evaluation of the gg → ZH amplitude via a p T -expansion, we compare the exact result for the LO partonic cross section [9,10] with the result obtained via our p T -expansion. The latter is expressed in terms of the same four MIs that enter into the analogous calculation of the gg → HH LO amplitude [18], or where are the Passarino-Veltman functions [24], with n the dimension of spacetime and µ the 't Hooft mass.
As an illustration of our LO result we present the explicit expressions for one symmetric, A 2 , and one antisymmetric, A 6 , form factor including the first correction in the ratio of small over large parameters which will be referred to as 1 O(p 2 T ). We divide the result into triangle ( ) and box ( ) contribution or and  In the lower part of fig.2 the ratio of the exact result over the p T -expanded one is shown. From this ratio one can see that the O(p 0 T ) contribution covers well the ZH invariant mass region M ZH 2m t , corresponding to the range of validity of an expansion in the large top quark mass. Furthermore, when the contributions up to O(p 4 T ) are taken into account a remarkable agreement with the exact result is found up to M ZH 750 GeV. This agreement is extended to sligthly higher values of M ZH when the O(p 6 T ) contribution is included, a finding in close analogy to the result for di-Higgs production [18]. Similar conclusions can be drawn from table 1, where it is shown that the partonic cross section at O(p 4 T ) agrees with the full result for M ZH 600 GeV on the permille level and the agreement further improves when O(p 6 T ) terms are included. As a final remark for this section, we notice that, from the comparison with the LO exact result, the p T -expanded evaluation of the amplitude is expected to provide an accurate result up to M ZH ∼ 700 − 750 GeV that corresponds, from eq.(7), to p T 300 − 350 GeV ≈ 2 m t .

Outline of the NLO Computation
In this section we discuss our evaluation of the three different types of diagrams that appear in the virtual corrections to the gg → ZH amplitude at the NLO. The triangle contribution ( fig.1(d), (e)) was evaluated using the observation of ref. [11], i.e. we adapted the result of ref. [22] for the decay of a pseudoscalar boson into two gluons to our case. This contribution is evaluated exactly and explicit expressions for the form factors are presented in appendix B. We notice that if we interpret the exact result in terms of our counting of the expansion in p T , the p T -expansion of the triangle contribution stops at O(p 2 T ). Given the reducible structure of the double-triangle diagrams ( fig.1(g)), an exact result for the double-triangle contribution can be derived in terms of products of one-loop Passarino-Veltman functions [24]. Explicit expressions for this contribution are presented in appendix B. Although we write the amplitude using a different tensorial structure with respect to ref. [16] we checked, using the relations between the two tensorial structures reported in appendix A, that our result is in agreement with the one presented in ref. [13].
The box contribution ( fig.1(f )) was computed evaluating the two-loop multi-scale Feynman integrals via two different expansions: a LME up to and including O(1/m 6 t ) terms, and an expansion in the transverse momentum up to and including O(p 4 T ) terms. The former expansion was used as "control" expansion of the latter. Indeed, the p T -expanded result actually "contains" the LME one. The LME differs from the expansion in p T by the fact that s is treated as a small parameter with respect to m 2 t , and not on the same footing as in the latter case. This implies that if the p T -expanded result is further expanded in terms of theŝ/m 2 t ratio the LME result has to be recovered. This way, we were able to reproduce, at the analytic level, our LME result.
We conclude this section outlining some technical details concerning our computation. We generated the amplitudes using FeynArts [25] and con-tracted them with the projectors as defined in appendix A using FeynCalc [26,27] and in-house Mathematica routines. We used dimensional regularization and the rule for the contraction of two epsilon tensors written in terms of the determinant of n-dimensional metric tensors. This is not a consistent procedure and needs to be corrected. A correction term should be added [28] to the form factors computed as described above, A In order to check eq.(24), following ref. [29] we bypassed the problem of the treatment of γ 5 in dimensional regularization computing the amplitude via a LME working in 4 dimension, employing the Background Field Method (BFM) [30] and using as regularization scheme the Pauli-Villars method. This result was compared with the LME evaluation of A (1,ndr) i , finding that the difference between the two evaluations was indeed given by the second term on the right-hand-side of eq. (24).
After the contraction of the epsilon tensors the diagrams were expanded as described in section 3. They were reduced to MIs using FIRE [31] and LiteRed [32]. The resulting MIs were exactly the same as previously found for di-Higgs production [18]. Nearly all of them are expressed in terms of generalised harmonic polylogarithms with the exception of two elliptic integrals [33,34]. The top quark mass was renormalized in the onshell scheme 2 and the IR poles were subtracted as in ref. [35].

NLO results
We now present our numerical results for the virtual corrections. We have implemented our results into a FORTRAN programme. For the evaluation of the generalised harmonic polylogarithms we use the code handyG [36], while the elliptic integrals are evaluated using the routines of ref. [34]. In order to facilitate the comparison of our results with the ones presented in the literature, we define the finite part of the virtual corrections as in ref. [16] 3 and in the numerical evaluation of eq.(25) we fixed µ R = √ŝ .
First, both the triangle and box LME contributions to A (1) i up to O(1/m 6 t ) terms were checked, at the analytic level, against the results of refs. [13,16] finding perfect agreement. Then, the p T -expanded results for low M ZH were confronted numerically with the LME ones, finding a good numerical agreement. We recall that, at the same order in the expansion, the p T -expanded terms are more accurate than the LME ones, although computationally more demanding.
In ref. [17] a numerical evaluation of eq.(25) was presented. In that reference the exact NLO amplitude was reduced to a set of MIs that were evaluated numerically using the code pySecDec [37,38]. Table 3 of that reference presents the numerical results 4 for various points in the phase space. For the four points in that table lying within the range of validity of our expansion we find a difference with respect to our results of less than 1 permille in 3 cases and reaching the 1 permille level in one case, similarly to what we find at LO. It should be noticed that small differences on the permille level can be explained not only by the different approaches (exact vs. p T -expanded) but also by the fact that in ref. [17] the m 2 Z /m 2 t and m 2 H /m 2 t ratios were approximated by a ratio of two integer numbers.
In order to present our results we define a virtual part of the partonic cross section from the finite part of the virtual corrections in eq.(25) by (26) and show it in fig.3. The dashed lines in the plot show the different orders in our expansion. For all parts of the matrix elements we use the best results available, i.e. both A (0) and the double-triangle contribution are evaluated exactly, while for A (1) we use the various orders in the p T -expansion. For comparison, we show the results where A (1) is replaced by the one computed in LME up to O(1/m 6 t ) (full black line), which as mentioned before is valid up to M ZH < 2m t . We see that within the validity of the LME our results agree well with it. Furthermore, we show the results in the infinite top mass limit reweighted by the full amplitudes squared (full red line), corresponding to the approach of ref. [11], keeping though the double triangle contribution in full top mass dependence. Differently from the LME line, the m t → ∞ reweighted one shows a behaviour, for M ZH 400 GeV, similar to the behaviour of the p T lines. Still, the difference between the reweighted result and the p T -expanded ones is significant. The p T -expanded results show very good convergence. The zero order in our expansion agrees extremely well with the higher orders in the expansion, and all the three results are very close up to M ZH ∼ 500 GeV.
Finally, we note that the evaluation of V f in requires a running time per phase space point less than one second. In addition, the integration over thê t variable in eq.

Conclusion
In this paper, we computed the two-loop NLO virtual corrections to the gg → ZH process. Among the two-loop Feynman diagrams contributing to the process, the ones belonging to the triangle and double-triangle topology were computed exactly. The ones belonging to the box topology, which contain multiscale integrals, were evaluated via an expansion in the Z transverse momentum. This novel approach of computing a process in the forward kinematics was originally proposed in ref. [18] for double Higgs production where the particles in the final state have the same mass. In this paper, we extended the method to the more general case of two different masses in the final state and to a process whose amplitude is not symmetric under the t ↔û exchange.
The result of the evaluation of the box contribution is expressed, both at one-and two-loop level, in terms of the same set of MIs that was found in ref. [18] for double Higgs production. The two-loop MIs can be all expressed in terms of generalised harmonic polylogarithms with the exception of two elliptic integrals.
As we have shown explicitly at the LO, the range of validity of our computation covers values of the invariant mass M ZH 750 GeV corresponding to 98.5% of the phase space at LHC energies. We showed that few terms in our expansion were sufficient to obtain an incredible good agreement with the numerical evalution of V f in presented in ref. [17], at the level of a permille or less difference between our analytic result and the numerical one.
The advantage of our analytic approach compared to the numerical calculation is also in the computing time. With an average evaluation time of half a second per phase space point, an inclusion into a Monte Carlo programme is realistic. Due to the flexibility of our analytic results, an application to beyond-the-Standard Model is certainly possible.
Finally, we remark that our calculation complements nicely the results obtained in ref. [16] using a high-energy expansion, that according to the authors provides precise results for p T 200 GeV. The merging of the two analyses is going to provide a result that covers the whole phase space, can be easily implemented into a Monte Carlo code and presents the flexibility of an analytic calculation.

A Orthogonal Projectors in gg → ZH
In this appendix we present the explicit expressions of the projectors P µνρ i appearing in eq. (2). The projectors are all normalized to 1. They are: where we defined q µ t = (p µ 3 − t s p µ 2 ) and q ν u = (p ν 3 − u s p ν 1 ) and we used the shorthand notation µνρp 2 ≡ µνρσ p σ 2 . Using these projectors we obtained the relations between the form factors Instead, for the triangle diagrams, we obtain , (B.10) , (B.14) where the K (2l) t function is defined in eq.(4.11) of ref. [22]. We do not show the explicit results for the p T -expansion of the two-loop box diagrams, since the analytic expressions are very lengthy, even for the lowest order term of the expansion.