Virtual corrections to gg → ZH in the high-energy and large-mt limits

We compute the next-to-leading order virtual corrections to the partonic cross-section of the process gg → ZH, in the high-energy and large-mt limits. We use Padé approximants to increase the radius of convergence of the high-energy expansion in mt2/s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {m}_t^2/s $$\end{document}, mt2/t\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {m}_t^2/t $$\end{document} and mt2/u\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {m}_t^2/u $$\end{document} and show that precise results can be obtained down to energies which are fairly close to the top quark pair threshold. We present results both for the form factors and the next-to-leading order virtual cross-section.


Introduction
At the CERN Large Hadron Collider (LHC), gluon fusion processes play an important role due to the large gluon luminosities at high collision energies.As a consequence one often observes that gluon fusion-induced processes provide a numerically large contribution to the theory predictions of production cross-sections.This is true even for processes for which the leading-order (LO) contribution consists of one-loop diagrams.A prime example of such a process is inclusive Higgs boson production, where the gluon-fusion channel is about an order of magnitude larger than all other contributions.
In this paper we consider the associated production of a Z and a Higgs boson, pp → ZH, often called "Higgs Strahlung".This process was of particular importance in the observation of the Higgs boson's decay into bottom quarks at ATLAS [1] and CMS [2].At LO it is mediated by a tree-level process in which a quark and anti-quark annihilate.For this channel, corrections up to next-to-next-to-next-to-leading (N 3 LO) order are available [3] (see also Ref. [4] and references therein) and are included in the program vh@nnlo [5,6].Electroweak and QCD corrections are also included in the program HAWK [7].
Associated ZH production can also occur via the loop-induced gluon fusion process.Although formally of next-to-next-to-leading order (NNLO) with respect to pp → ZH, it provides a sizeable contribution, in particular in the boosted Higgs boson regime in which the transverse momentum of the Higgs boson is large [8,9].Furthermore, the process gg → ZH provides sizeable contributions to the uncertainties of ZH production with a subsequent decay of the Higgs boson in a pair of bottom quarks, e.g., Ref. [10].Being loop-induced, gg → ZH is sensitive to physics beyond the Standard Model.In Ref [11] it was suggested that the gluon-initiated component of pp → ZH can be extracted by comparing to W H production, which only has a Drell-Yan-like component.It is thus important to consider NLO QCD corrections to gg → ZH, requiring the computation of two-loop box-type Feynman diagrams with two different final-state masses (m Z and m H ) and the massive top quark propagating in the loops.
An exact LO (one-loop) calculation was performed in [12].At NLO only approximations in the large m t limit are known [13,14].In this work we consider the two-loop NLO virtual corrections in the high-energy limit, expanding the two-loop master integrals for m 2 t s, t, where s and t are the Mandelstam variables.Furthermore, we also provide analytic results for the form factors in the large-m t limit. 1 A similar approach has been applied to the related process of Higgs boson pair production, gg → HH, where a comparison to exact numerical calculations [15] could be performed and good agreement was found, even for relatively small values for the Higgs boson transverse momentum [16].In Ref. [17] the high-energy expansion was successfully applied to the top quark contribution of the twoloop diagrams contributing to gg → ZZ.
The remainder of the paper is organized as follows: In Section 2 we introduce our notation, and give our definitions for the form factors and the virtual finite cross-section.
In Section 3 we consider the quality of our approximations by comparing with the exact LO expressions.In Section 4, we briefly discuss the two-loop form factors in the large-m t limit and in Section 5 we discuss the form factors in the high-energy limit, and investigate the behaviour of Padé approximants constructed using this expansion.In Section 6 we study the NLO virtual finite cross-section and apply our Padé scheme to extend the approxmation to a larger kinematic region.Finally in Section 7 we conclude our findings.Auxiliary material can be found in the Appendix.In Appendix A we present analytic results for the one-particle reducible double-triangle contribution and in Appendix B we briefly discuss out treatment of γ 5 and the application of projectors to obtain the form factors.In Appendix C we present the relations between the form factors and helicity amplitudes for the process gg → ZH.The expansions of the form factors are provided in an analytic, computer readable form in the ancillary files of this paper [18].

Notation and technicalities
We consider the amplitude for the process g(p 1 )g(p 2 ) → Z(p 3 )H(p 4 ) where all momenta are assumed to be incoming.This leads to the following definitions for the Mandelstam variables, with Additionally, We also introduce the transverse momentum (p T ) and the scattering angle (θ) of the final-state bosons, which are related to the Mandelstam variables as follows, where We denote the polarization vectors of the gluons and the Z boson by ε λ 1 ,µ (p 1 ), ε λ 2 ,ν (p 2 ) and ε λ 3 ,ρ (p 3 ), in terms of which the amplitude can be written as Due to charge-conjugation invariance, the vector coupling of the Z boson to the quarks in the loop does not contribute.The axial-vector coupling is proportional to the third component of the isospin and thus mass-degenerate quark doublets also do not contribute.Since we consider the lightest five quarks to be massless, only contributions from the topbottom doublet remain.Furthermore, each individual term of A µνρ is proportional to the totally anti-symmetric ε tensor from the axial-vector coupling.
At LO and NLO both triangle-and box-type diagrams have to be considered.Examples of these diagram classes are shown in Fig. 1.In the box-type diagrams the Higgs boson couples directly to the quark loop, so diagrams involving the bottom quark are suppressed by their Yukawa coupling with respect to diagrams involving the top quark.This is not the case for the triangle-type diagrams; here contributions from both the top and bottom quark loops must be considered.At NLO there is also the contribution from the one-particle reducible double-triangle contribution, shown as first diagram in the second row of Fig. 1.Note that in the numerical results discussed in the main part of this paper these contributions are excluded, however, we present exact analytical results in Appendix A. The contribution from the reducible double-triangle diagrams to the (differential) partonic cross-section is implemented in the computer program which comes together with Ref. [14].
In total one can form 60 tensor structures from the indices µ, ν, ρ, the independent momenta p 1 , p 2 and p 3 and an ε tensor.Details of the computation of these 60 structures and our treatment of γ 5 in d = 4 − 2 dimensions are given in Appendix B. After applying transversality conditions (p i • ε λ i (p i ) = 0), gauge invariance w.r.t. the gluons (p 1µ A µνρ = p 2ν A µνρ = 0) and Bose symmetry (A µνρ (p 1 , p 2 , p 3 ) = A νµρ (p 2 , p 1 , p 3 )), 14 tensor structures remain which can be grouped such that one has to introduce 7 form factors.
We follow the decomposition of Ref. [12] and write where a, b are colour indices.Note that while the decomposition is the same, the form factors in Ref. [12] have dimension 1/GeV 2 whereas we pull out an overall factor of 1/s such that our form factors are dimensionless.Only F 1 receives contributions from the triangle-type diagrams discussed above.We represent the expansion coefficients of the form factors in the strong coupling constants with the following notation, At this point it is convenient to define new form factors which are linear combinations of those of Eq. ( 6), It is easy to see that F + 2 (t, u) drops out in the squared amplitude and thus does not contribute to physical quantities.It is furthermore convenient to introduce leaving six physically relevant functions: where only F + 12 (t, u) has contributions from triangle-type diagrams.F − k (t, u) (with k = 12, 2, 3) and F 4 (t, u) are anti-symmetric w.r.t. the exchange of the arguments t,u, and F + k (t, u) (with k = 12, 3) are symmetric.At leading order F − 3 (t, u) = 0, however it has non-zero contributions starting at NLO.
For the computation of the one-and two-loop Feynman diagrams (some examples are shown in Fig. 1) in the high-energy limit, we proceed as follows: After producing the amplitude we perform a Taylor expansion in the Z and Higgs boson masses (since t ), leaving one-and two-loop integrals which depend only on s, t and m t .Using integration-by-parts reduction techniques with the programs FIRE 6 [19] and LiteRed [20], these integrals can be reduced to the same basis of 161 two-loop master integrals as in Refs.[21,22], in which they were computed as an expansion in the high-energy limit to order m 32 t .Inserting these expansions into the amplitude yields its high-energy approximation.We use FORM 4.2 [23] for most stages of the computation.The calculation is performed in the covariant R ξ gauge and we allow for a general electroweak gauge parameter ξ Z which appears in the propagators of the Z boson and Goldstone boson χ.Thus only the triangle-type diagrams depend on ξ Z , and this dependence cancels upon summing the Z and χ contributions.
We also investigate the large-m t limit of the form factors.This expansion is straightforward and proceeds in analogy to [17,22].The programs q2e and exp [24,25] are used to produce an asymptotic expansion for m t q 1 , q 2 , q 3 , again performed using FORM, yielding products of massive vacuum integrals and massless three-point integrals.Results for the gg → ZH amplitude, expanded to order 1/m 8 t , have been previously published in [14]; here we provide one additional term in this expansion, to order 1/m 10 t , and provide results at the level of the form factors.
The renormalization of the ultra-violet (UV) divergences proceeds in a similar way as for the processes gg → HH [22] and gg → ZZ [17].In particular, we work in the six-flavour theory and renormalize the top quark mass on shell and the strong coupling α s in the MS scheme.In addition, our treatment of γ 5 (see Appendix B for more details) requires that we apply additional finite renormalization to the axial and pseudo-scalar currents [26], The subtraction of the infra-red poles proceeds according to Ref. [27], using the conventions of Refs.[17,22].The subtraction has the form [27] with β 0 = 11C A /12 − T n f /3 and F (1),IR is UV renormalized but still infra-red (IR) divergent.F (1) is finite.An explicit expression for K g is given in Eq. (2.3) of Ref. [22].
After the second equality sign in Eq. ( 11) we make the µ dependence explicit.Note that below we only need F (1) to construct the NLO cross-section.
In analogy to the process gg → HH [28] we define the NLO virtual finite cross-section for gg → ZH as where the form factors entering the amplitude Ãµνρ sub are the IR-subtracted finite form factors F (1) of Eq. ( 11).The superscripts "(0)" and "(1)" after the square brackets in Eq. ( 12) indicate that we take the coefficients of (α s /π) 0 and (α s /π) 1 , respectively, of the squared amplitude, in accordance with Eq. ( 7).For the discussion in Section 6 it is convenient to introduce the α s -independent quantity

Comparison at leading order
In this section we compare our high-energy expansion with the exact LO result.We implement this by using FeynArts [29] to generate the amplitude and FormCalc [30] to compute it, performing a Passarino-Veltman reduction to three-and four-point one-loop scalar integrals.Schouten identities allow us to write the result in terms of the tensor structures and form factors of Eqs. ( 6) and (8).We use Package-X [31] to evaluate the Passarino-Veltman functions with high precision and to produce an analytic expansion in the limit of a large top quark mass.We have verified that our implementation of the exact LO result reproduces the cross-sections provided in the literature [13,14,32].Additionally we have compared the large-m t expansion derived from this result with a direct expansion of the amplitude as described above.
In Fig. 2 we show the squared amplitude at LO, as a function of √ s, for a fixed scattering angle θ = π/2.The solid black lines correspond to the exact result.The coloured lines correspond to different expansion depths in m Z and m H , as detailed in the plot legend.All of the latter include high-energy expansion terms up to m 32 t .The red curve demonstrates that after including quadratic terms both in m H and m Z , the deviation from the exact result is below 1% (for √ s ∼ > 1000 GeV).The m 2 H corrections appear to be more important than the m 2 Z corrections.The inclusion of the quartic terms (shown as green and pink curves, best visible in Fig. 2(c)) improves the accuracy of the expansion, leading to an almost negligible deviation from the exact expression.These terms are much  harder to compute, however, so in the NLO results we will restrict the approximation to the quadratic corrections only.The solid blue curve and associated uncertainty band in Fig. 2(a) shows the result of a procedure to improve the expansions based on Padé approximants, which is discussed in more detail in Section 5.Here it is based on the expansion to quadratic order in m H and m Z , confirming that our computation of the NLO amplitude only to this order is sufficient.
In Fig. 2 and in the following sections, we use the following parameter values: 4 NLO form factors: large-m t limit In this section we discuss the large-m t expansion of the form factors at NLO.While the expansion of F + 12 starts at m 0 t , we find that the other form factors of Eq. ( 10) exhibit some cancellation in the leading contributions.In particular, F − 12 , F − 2 , F + 3 and F4 start only at 1/m 4 t , and F − 3 starts at 1/m 6 t .We now present the leading terms of the large-m t expansion to establish our notation.For brevity, we restrict ourselves to F + 12 .Expressions for the expansion of all form factors to 1/m 10 t are provided in the ancillary files of this paper [18].Our result reads where C A = 3 and C F = 4/3 are the quadratic Casimir invariants of SU (3).
We have verified that after constructing Ṽfin in Eq. ( 12) we find agreement with the results of Ref. [14] up to order 1/m 8 t .

NLO form factors: high-energy limit
We now turn to the high-energy expansion of the NLO form factors.The analytic expressions are large, so we refrain from showing them here but we provide analytic expressions for all form factors in the ancillary files of this paper [18].For illustration we discuss in the following the results for F + 3 .In Fig. 3(a) we plot F + 3 as a function of √ s.We can see that the high-energy expansions of both the real (blue solid lines) and imaginary parts (green dashed lines) converge well for √ s ∼ > 800 GeV, which is in analogy to gg → HH [22] and gg → ZZ [17].Note that the lighter-coloured curves include fewer m t expansion terms; the darkest lines show the expansion up to m 32 t .As in previous publications on gg → HH and gg → ZZ, we make use of Padé approximants to improve the description of the high-energy expansion.The methodology used follows that of Section 4 of Ref. [17] and so is not described in detail, but is only summarized briefly here.The expansion is used to produce 28 different Padé approximants, which are combined to produce a central value and error estimate for the approximation.In Fig. 3 Padé-based approximation reproduces the expansion.However, it also produces reliable results for smaller values of √ s, as can be expected from the comparison with the LO result shown in Fig. 2(a).
In Fig. 3(b) we show F + 3 for the fixed values of p T = 200, . . ., 800 GeV.The blue (dashed and dotted) curves correspond to the high-energy expansions and the green (dashed) curves to the Padé results.For all Padé curves we also show the corresponding uncertainty band, which for p T = 200 GeV is relatively large but for p T = 250 GeV the uncertainty band is already quite small; it is completely negligible for higher values of p T .Note that the high-energy expansions are only shown for p T ≥ 400 GeV; for lower p T values the curves lie far outside of the plot range.For p T ∼ > 450 GeV the expansions converge and are very close to the Padé results.For p T = 400 GeV, while the expansions initially agree with each other and the Padé close to √ s = 800 GeV, they diverge for larger values of √ s.We recall here that the high-energy expansion is an expansion in m 2 t /s, m 2 t /t and m 2 t /u.For a fixed value of p T , increasing √ s can lead to values of t or u which are not large enough for convergence.
In Section 6 we will apply the Padé procedure to the virtual finite cross-section, in order to compare our results with a state-of-the-art numerical evaluation at NLO [33].

NLO virtual finite cross-section
Our starting point is Ṽfin in Eq. (12). the LO form factors we use the exact results and for the two-loop form factors we use the high-energy expansion.This allows us to write V fin in Eq. ( 13) as an expansion in m t up to order m 32 t . 2 At this point we can apply the Padé approximation procedure as described in Section 5.In Fig. 4 we show V fin for several fixed values of p T as a function of √ s.For p T = 400 GeV and larger we also show two high-energy expansion curves, which include terms up to m 30 t and m 32 t .These curves agree well with each other and with the Padé approximation which they produce.For lower values of p T , these curves do not agree with each other, and are not visible within the plot range.
For p T = 150 GeV Padé procedure produces stable results with an uncertainty of about 10%.For p T = 200 GeV the uncertainties are notably smaller and for higher p T values they are negligible.We advocate to use results based on our approach for p T ∼ > 200 GeV.For p T ≈ 150 GeV the Padé approach provides important results for cross-checks.For lower values of p T other methods have to be used for the calculation of gg → ZH, see Ref. [33].
We have compared our results with those of Ref. [33], which are obtained numerically, but without making any expansions.In the high-energy region, for 104 kinematic points with p T ≥ 200 GeV we agree to within 2-3%, and for 42 points with 150 ≤ p T < 200 GeV we agree to within 10% and our values are consistent to within the uncertainty of the Padé procedure.We also construct V fin using the large-m t expanded NLO form factors of Section 4 and find that for 120 points with √ s ≤ 284 GeV, we agree to within 1%.For larger values of √ s, as expected, the large-m t expansion starts to diverge significantly from the numerical results as one approaches the top quark threshold at √ s ≈ 346 GeV.

Conclusions
In this paper we consider two-loop NLO corrections to the Higgs Strahlung process gg → ZH.The corresponding Feynman integrals involve five dimensionful parameters (s, t, m t , m Z and m H ) which makes an analytic calculation impossible with the current state-ofthe-art techniques.Numerical computations are also very challenging, however they have recently been achieved in Ref. [33].In Section 6 we have discussed the cross-check of our results against these.
In our approach we use the hierarchy between the top quark, Higgs and Z boson masses and perform an expansion for m 2 t m 2 Z , m 2 H , effectively eliminating the dependence on m Z and m H from the integrals.We show at one-loop order that the expansion converges quickly, which allows us to truncate the expansion at NLO after the quadratic terms.As far as the remaining scales are concerned, we concentrate on the high-energy region where s, t m 2 t .We expand our master integrals in this limit such that we obtain results for the form factors including expansion terms up to m 32 t .This allows us to construct, for each phase-space point (e.g., a particular pair of √ s, p T values) a set of Padé approximants, which considerably extend the region of convergence of our expansion.Our approach provides both central values and uncertainty estimates.
We provide results for all form factors involved in the gg → ZH amplitude, both at one-and two-loop order.Furthermore, we provide relations between the form factors and helicity amplitudes in Appendix C. The main emphasis of the paper is the IR subtracted virtual corrections to the partonic cross-section, which can be combined with other (e.g.Drell-Yan-like) corrections to gg → ZH.Our method provides precise results for p T ∼ > 200 GeV with almost negligible uncertainties.For lower values of p T our uncer-tainties increase, in particular for smaller values of √ s.The analytic results for the form factors obtained in this paper can be obtained in electronic form from Ref. [18].
The contribution of the double-triangle diagrams to the squared matrix element have been computed in Ref. [14] and implemented in the corresponding computer program, see Appendix of Ref. [14].

B Projectors and γ 5
Each LO and NLO diagram contributing to the gg → ZH amplitude contains a quark with either an axial-vector coupling to a Z boson, or a pseudo-scalar coupling to a Goldstone boson.The γ 5 matrix present in these couplings is not defined in the d = 4 − 2 spacetime dimensions of dimensional regularization.The couplings are re-written in terms of anti-symmetric tensors, according to [26], as which allow one to compute the loop integrals in two ways.The first is to ensure the ε tensors are not contracted until one can safely work in four dimensions, by solving tensor loop integrals, performing UV renormalization and IR subtraction (as detailed in Eq. ( 11)) and then finally contracting with the ε tensors.
An alternative approach, which avoids the need to compute tensor integrals, is to project the amplitude onto the 60 possible rank-three Lorentz structures which can be formed from the three independent momenta and the ε tensor.Since each term of the projectors onto these structures contains an ε tensor, the result of their contraction with the r.h.s. of Eq. ( 18) can be treated correctly in d dimensions.The 60 projectors act on the amplitude as to produce 60 form factors F i , which can be reduced to a minimal set (see the discussion around Eq. ( 6)).Each projector P µ 1 µ 2 µ 3 i can be written generically as and we contract this general form with our amplitude.In the final result, we specify each of the 60 sets of coefficients {C i,j } in order to obtain the 60 form factors F i .
We furthermore specify the (contravariant) external momenta and polarization vectors as follows: where ε 0 denotes the longitudinal components of polarization vectors.Recall that all external momenta are defined as incoming, see Eq. (1).We have chosen the convention for the polarisation vector of p 2 , following Ref.[35], such that ε + (p 2 ) → ε − (p 1 ) in the center-of-momentum frame.Furthermore, the polarization vectors satisfy which means that we have fixed the gauge for the external gauge bosons.With the above choice, some of the tensor structures in A µνρ ab (p 1 , p 2 , p 3 ), which are proportional to either p ν 1 or p µ 2 , are irrelevant because This reduces the number of Lorentz structures in Eq. ( 6) from 14 to 8. We note that, as we will see in Eq. ( 28), the dependence on F + 2 drops out in the helicity amplitude.

Figure 2 :
Figure 2: (a) LO squared amplitude and (b) ratio of expansions to the exact result.Plot (c) is a zoomed-in version of (b).Note that for better readability in (a) no quartic corrections are shown; the black curve in (a) refers to the exact result and the blue curve (and associated uncertainty band) is the result obtained from Padé approximation.

Figure 3 :
Figure 3: Results for F + 3 as a function of √ s for fixed θ = π/2 (a) and for fixed p T (b).In (a) dashed and solid lines correspond to the imaginary and real parts, respectively.The red curves in (a) represent the Padé results.In (b) only the real part is shown, and the Padé results are shown as green dashed lines.The high-energy expansions up to order m 30 t and m 32 t are shown in blue.The widely-separated pair of curves correspond to p T = 400 GeV.

32 tFigure 4 :
Figure 4: V fin as a function of √ s for eight values of p T , at the scale µ 2 = s.The uncertainty estimate of the Padé procedure is displayed as light-coloured bands.For p T ≥ 350 GeV, we also show two high-energy expansion curves including terms up to order m 30 t and m 32 t .