$gg\to ZZ$: analytic two-loop results for the low- and high-energy regions

We compute next-to-leading order virtual two-loop corrections to the process $gg\to ZZ$ in the low- and high-energy limits, considering the contributions with virtual top quarks. Analytic results for all 20 form factors are presented including expansion terms up to $1/m_t^{12}$ and $m_t^{32}$. We use a Pad\'e approximation procedure to extend the radius of convergence of the high-energy expansion and apply this approach to the fini\ te virtual next-to-leading order corrections.


Introduction
The remainder of the paper is organized as follows. In the next section we introduce our definitions and notation, and describe our methodology for the computation of the highenergy and large-m t expansions. We also discuss how one can obtain helicity amplitudes from our form factors. In Section 3 we compare the expansions to the exact LO result and justify our choices for the expansion depths used at NLO. In Section 4 we describe how one can improve the radius of convergence of the high-energy expansions by making use of Padé approximants. Using this method, in Section 5 we show NLO results for form factors and for the finite virtual corrections to the cross section. For the latter, we consider different values for the transverse momentum of the Z bosons and demonstrate that we can obtain stable predictions for this quantity for transverse momentum values as small as 150 GeV. Our conclusions are presented in Section 6. In the Appendix we provide the explicit results for the relations which can be used to rotate to the orthogonal tensor basis of Section 2.3. Furthermore, numerical results for all LO and NLO form factors and analytic results for some example LO form factors are presented.

Technical details
In Fig. 1 we show one-and two-loop sample Feynman diagrams contributing to process where all momenta p i are incoming. The Mandelstam variables are defined as and fulfil the property For later convenience we also introduce the velocity β and the transverse momentum of the Z bosons as where θ is the scattering angle.
In this paper we consider only the top quark as the virtual particle in the loop. We exclude from our analysis the two-loop contribution which originates from the product of two one-loop triangle diagrams (the so-called anomaly contribution) since this contribution is discussed in detail in Ref. [20], in which exact results are presented.
coupling, and we note that only f (j) 1,tri is non-zero. In the case of massless quark loops, f (j),vt i,box /v 2 t = f (j),at i,box /a 2 t [15]. This property is satisfied by the leading term of our high-energy expansions (m 0 t m 0 Z ) but is violated in higher order terms, including for higher order terms in m Z since m Z < m t .
These form factors are, at this point, divergent in 4 dimensions. We perform the renormalization of the top quark mass, the strong coupling constant and the gluon field to remove the ultra-violet divergences. The remaining divergences are infrared in nature and are removed by the subtraction procedure of Ref. [28], which we outline here.
We construct finite form factors which are defined as where f (1),IR is ultraviolet renormalized but still infrared divergent. K (1) g can be found in Ref. [28] and is given by where γ E is Euler's constant. Note that the poles in the terms proportional to n f from Eq. (14) cancel against the counterterm contribution induced by the α s renormalization. However, finite terms proportional to log(µ 2 /(−s − iδ)) remain, which can be cast in the form , which are independent of µ, contain new information and thus only they will be discussed in Section 5.
We now discuss the work-flow for our calculation of these form factors, as expansions in both the high-energy (Section 2.1) and large-m t (Section 2.2) limits. Analytic expressions for the results of both of these expansions can be found in the ancillary file of this paper [29]. In both cases, the amplitude is generated using qgraf [30]. Each Feynman diagram is then contracted with one of the 138 possible tensor structures discussed above, as a separate computation. This splitting is particularly important for the large-m t expansion of Section 2.2, in order to avoid overly large intermediate expressions.
We additionally reproduce the exact LO result from [15] using the programs FeynArts 3.10 [31] and FormCalc 9.8 [32]. The scalar Passarino-Veltman functions B 0 , C 0 and D 0 are rewritten in terms of polylogarithms with the help of Package-X [33], which allows for a high-precision evaluation within Mathematica. We use this exact LO result to evaluate the performance of our expansions and approximation methods in Section 3.

High-energy expansion
For each contraction we compute the fermion traces and write the result in terms of scalar Feynman integrals, belonging to one of the integral families defined in Refs. [34,35] (there in the context of an NLO calculation of gg → HH in the high-energy limit). We then construct the appropriate linear combinations which are required to obtain the form factors of the 20 tensor structures given in Eq. (9). Up to this point our calculation is exact in all kinematic variables and masses.
Next, we Taylor expand both the scalar Feynman integrals and their coefficients in m Z , using the program LiteRed [36] and in-house FORM [37] routines. For each integral family we perform an integration-by-parts (IBP) reduction to master integrals using version 6 of FIRE [38] and symmetry relations obtained using LiteRed [36]. Since we have performed a Taylor expansion the integrals depend on the kinematic variables and m t , but no longer on m Z ; this makes the IBP reduction much more tractable. For the most complicated family (numbered 91 in Appendix A of Ref. [35]) this takes about 4.5 days 1 on a 3.5 GHz machine with 32 cores.
Inserting the reduction tables into the amplitude and expanding the resulting expressions in m t and took around three weeks on a reasonably sized cluster of computers. Using the results for the master integrals of Refs. [34,35], we produce an expression for the amplitude expanded up to m 32 t and m 4 Z . The coefficients of the expansion terms are functions of s and t, and are written in terms of Harmonic Polylogarithms with a harmonic weight of at most 4, for the numerical evaluation of which we use the package HPL.m [39].
The expansions contain terms with both even and odd powers of m t ; the odd powers come from the expansions of the two-loop non-planar master integrals. Most of the odd powers cancel in the amplitude, however starting from m 3 t , odd m t powers remain in the imaginary part of the non-abelian contribution to the form factors. The situation is analogous to gg → HH [35] where the contributions of odd powers is discussed in detail at the level of master integrals.

Large-m t expansion
We keep the discussion of the large-m t expansion brief, since the methods are largely the same as used in the expansion of Higgs boson pair production, described in detail in Ref. [40]. At the level of the individual Feynman diagrams contracted with one of the 138 possible tensor structures, we apply an asymptotic expansion for m t p 1 , p 2 , p 3 using the program exp [41,42].
This leads to one-and two-loop vacuum integrals with the scale m t multiplied by massless three-point integrals with the scale s. We expand, at one and two loops, to order 1/m 12 t .
After expansion, we compute the appropriate linear combinations of the contractions in order to arrive at the coefficients of the 20 tensor structures of Eq. (8), yielding the largem t expanded expressions for the form factors defined in Eq. (12). For the convenience of the reader we show the leading terms in the 1/m t expansion for some form factors in Appendix D.
In Ref. [19] the amplitude for gg → ZZ has been calculated at LO and NLO up to the first non-vanishing expansion term in 1/m t , which only involves the axial-vector part. We find agreement after fixing two obvious typos 2 . Furthermore, in Ref. [20] analytic results for the gg → ZZ amplitude projected to the triangle contribution are presented as an expansion up to order 1/m 12 t . After performing the same projection we could successfully compare our results for the vector and axial-vector part which constitutes a welcome check for our approach.

Orthogonal tensor basis
The tensors given in Eq. (9) have the advantage of being simple and compact. However, they are not orthogonal; this leads to non-vanishing cross terms when squaring the amplitude. For this reason we construct a new basis T i , using the Gram-Schmidt orthogonalization procedure, which has the property For d = 4 the coefficients c i are given by Note that in four dimensions c 19 and c 20 vanish which means that in the orthogonal basis only 18 form factors contribute to the final results. To obtain Eq. (17) we have made use of the polarization sums already listed in Eq. (10).
The basis change from S i to T i is described in Appendix A. In terms of T i the amplitude in Eq. (8) reads where the factors 1/ √ c i have been introduced such that the coefficients F i are dimensionless. As for f i , the coefficients F i have a decomposition into form factors F should be multiplied by (−1) and in Eq. (7) p µ 1 and p ν 2 should be replaced by p µ 2 and p ν 1 , respectively. and F (j),at i,box as described in Eq. (12). It is in terms of these form factors, of the orthogonal basis of tensor structures, that we write an expression for the differential cross section: where "R" denotes the corrections due to real radiation which we do not consider here. The basis change is computed numerically, upon evaluation of the differential cross section for particular values of the kinematic parameters.

Helicity amplitudes
In this subsection we describe how one can obtain the helicity amplitudes for the process gg → ZZ from the tensor decomposition which we have introduced above. For this purpose it is convenient to explicitly specify the external momenta and to introduce polarization vectors as follows: 3 where ε 0 denotes the longitudinal components of polarization vectors. Recall that all external momenta are defined as incoming and that the polarization vectors are chosen such that they satisfy Eq. (10). The helicity amplitudes M λ 1 ,λ 2 ,λ 3 ,λ 4 are given by Eq. (7). In total there are 2 × 2 × 3 × 3 = 36 helicity amplitudes. However, due to various symmetries only eight of them are independent. First, due to which holds for p, p = p 1 , ..., p 4 , we have which reduces the number of independent amplitudes to 18. Furthermore, there are additional symmetries [15] relating helicity amplitudes with different polarization states and there are symmetry relations due to β → −β, Note that this replacement changes none of the Mandelstam variables s, t, u. Using Eqs. (23) and (24) reduces the number of independent helicity amplitudes by six and four, respectively, and we arrive at eight independent helicity amplitudes.
It turns out that the above symmetries are fulfilled when the form factors satisfies the relations Note that up to this point we do not make use of any approximation. We use the relations (25) as a cross check of our calculations.
The LO results for the eight independent helicity amplitudes are provided in Ref. [15] and we confirm the agreement between them and our results. 4 For the results of the high-energy expansion we have expanded in m Z , making the symmetry relations due to β → −β hard to realize, since β = 1 − 2m 2 Z /s + O(m 4 Z ) and we do not distinguish the origin of m Z terms in the expression. For this reason, in the ancillary file [29] we provide results for the twelve helicity amplitudes

Comparison at leading order
This section is devoted to the discussion of the LO contribution to gg → ZZ with virtual top quarks. We first consider the form factors and helicity amplitudes, and compare the large-m t and high-energy expansions with the exact results. Afterwards we discuss our approach to improve the radius of convergence of our expansions, which is based on Padé approximations. We furthermore investigate the importance of finite Z boson mass corrections. For the numerical evaluation we use the following input values [43] As typical examples for the LO form factors, in Fig. 2 we show the results for F    Figure 3: LO partonic differential cross section for gg → ZZ for θ = π/2. Solid, dashdotted and dotted lines correspond to the exact, high-energy and large-m t results.  We show a selection of expansion depths between m 2 t and m 32 t . One observes that five to six expansion terms are necessary in order to obtain a good approximation of the exact result for √ s 1000 GeV. The deeper expansion depths show agreement down to √ s ≈ 750 GeV, which cannot be further improved even by including terms up to m 32 t . It appears that the simple expansions in m 2 t /s, m 2 t /t and m 2 t /u have a finite radius of convergence, which for θ = π/2 manifests itself around √ s ≈ 750 GeV. This feature can be understood by inspecting the functions which are present in the exact one-loop result. Among others we have identified logarithms and di-logarithms which depend on the quantity which has, in the high-energy limit, a radius of convergence of ut/s = 4m 2 t . For θ = π/2 we have t = u = −s/2 which leads to √ s = 4m t 700 GeV.
Let us next discuss the importance of finite m Z terms. The high-energy approximations shown in Figs we show how the number of expansion terms in m 4 Z affects the quality of the expansion of the LO differential cross section. Curves including terms to m 32 t and m 0 Z , m 2 Z and m 4 Z are shown, normalized to the exact result. For all three curves we observe, as discussed above, a divergent behaviour for √ s 750 GeV. The m 0 Z curve shows a more than 5% deviation from the exact result and including the m 2 Z term leads to a significant improvement, with the deviation reducing to around 1%. Finally, including the m 4 Z term produces a permille level agreement with the exact result, which motivates our computation of the m 4 Z expansion terms of the NLO quantities discussed in later sections of this paper.

Padé-improvement of the high energy expansion
In Section 3 we investigated the behaviour of the expansions and, in particular, noted that the high-energy expansion fails to converge below √ s ≈ 750GeV regardless of how many expansion terms are included. In this section we discuss a method by which we can extend the prediction of the high-energy expansion to smaller values of √ s.
The method is an extended version of the approach used in Ref. [44] in the context of Higgs boson pair production, and we describe it in detail below. It is based on the construction of a number of Padé approximants using the terms of the high-energy expansion, and subsequently combining the approximants to produce a central value and uncertainty estimate for a given phase-space point { √ s, p T }. We describe the procedure in terms of a generic quantity F for which we assume an expansion in m t is available. F also depends on the kinematic quantities s and p T and on m Z . In our practical applications F can be either a form factor, a helicity amplitude or the virtual finite cross section defined in Section 5.
The approximation procedure for F is then as follows: • We write F as an expansion in m t and define where F 0 contains the exact (in m t and m Z ) expressions of the LO contributions. x k for the odd and even powers of m t . We insert numerical values for m t , m Z , s and p T , yielding a polynomial in the variable x.
• Next we construct Padé approximants of F N in the variable x and write F N as a rational function of the form The coefficients a i and b i are determined by comparing the coefficients of x k after expanding the right-hand side of Eq. (30) around the point x = 0. Evaluation of this rational function at x = 1 yields the Padé approximated value of F N .
The numerator and denominator degrees in Eq. (30) are free parameters; one only must ensure that n + m ≤ N/2 such that a sufficient number of expansions terms are available to determine the coefficients a i and b i . We construct many Padé approximations and combine them to obtain a prediction for the central value and the uncertainty of F.
The rational function of Eq. (30) develops poles which, for some Padé approximants, might lie close to the evaluation point x = 1 and yield unphysical results. In the following we describe a weighting scheme which minimizes the influence of such Padé approximants. We call this approach a pole distance reweighted (PDR) Padé approximation.
• For each phase-space point { √ s, p T } we compute, for each Padé approximant, the value at x = 1 and the distance of the nearest pole which we denote by α i and β i , respectively.
• Introduce a weighting function, which reduces the impact of values α i from Padé approximations with poles close to x = 1. We define where the sum runs over all Padé approximants under consideration.
• Use the values α i and ω i,poles to compute the weighted average and weighted standard deviation of the Padé approximants, These form the central value and error estimate of the approximation.
At this point, the procedure is the same as that of Ref. [44], in which expansions up to m 30 t and m 32 t were used to create Padé approximants with 15 ≤ n + m ≤ 16, with the additional restriction to "near-diagonal" approximants which satisfy |n − m| ≤ 2. This results in 5 possible approximants, [7/9], [8/8], [9/7] which were weighted according to the above procedure to produce a central value and error estimate.
In this paper we further refine the method which allows us to loosen the restrictions and thus include more approximants in the computation. We introduce two additional weights into the averaging procedure which a) emphasize the contribution from Padé approximants which are derived from a larger number of expansion terms and b) emphasize the contribution from "near-diagonal" approximants. These weights are defined as follows: a). An [n i /m i ] Padé approximant is weighted by b). An [n i /m i ] Padé approximant is also weighted by As above, the sums run over all Padé approximants under consideration. The weights of Eqs. (31), (34) and (35) are combined according to and used to form a central value and error estimate The approximation of Eq. (32) used in Ref. [44], with the restrictions described above, can be considered to be a special case of the same procedure with the weights of Eqs. (34) and (35) replaced with step functions. In this refined procedure we include a wider set of Padé approximants. We define the quantities N low and N high such that In Section 5 we will study the quality of the approximations due to the choices of In Fig. 5 we demonstrate the effect of including higher order terms in the expansion in m Z in the construction of the Padé approximants for the LO differential cross section dσ/dθ. We show plots for p T = 150 GeV and p T = 200 GeV and, from top to bottom, approximations formed from high-energy expansions which include m 0 Z , m 2 Z and m 4 Z terms. One observes that it is crucial to include corrections at least to order m 2 Z , and that the results are further improved by including additionally the m 4 Z terms. These improvements are in line with the expectations due to the behaviour of the high-energy expansion as demonstrated in Fig. 4. We note that in the p T = 200 GeV plots, the majority of the points have error bars which are too small to be visible. After including the higher order m Z terms, the exact results lie within the error estimates of the approximations, demonstrating that they are realistic.
The bottom-left plot of Fig. 5 we additionally show, in black, a Padé approximation according to the simpler prescription of Eq. (32) using the five Padé approximants of the set Eq. (33). One observes that for small values of √ s the exact result lies outside of the error estimates and that for large values of √ s the errors appear to be overestimated. In our view the purple points provide a more reasonable description of the uncertainty.
It is interesting to have a closer look at the effective expansion parameters entering the high-energy expansion; the final result is expressed as an expansion in m 2 t /s, m 2 t /t and m 2 t /u.

NLO results for gg → ZZ with virtual top quarks
In this Section we apply the approximation procedure of Section 4 to our results for NLO quantities. We begin by considering two example form factors, in the orthogonal basis of Section 2.3: F 1 , the plot suggests that for real and imaginary parts the curves merge smoothly into the large-m t expansion. In the case of F (1) 16 we expect a resonance-like structure (as for F (0) 16 in Fig. 2) which is also indicated by the Padé curves. The thickness of the solid curves reflect our estimate of the uncertainty due to the approximation procedure, as defined in Eq. (37).
Although not the focus of this paper, the form factors also receive contributions from massless quarks running in the loop. We extend the notation of Section 2 and define where for up-type quarks v f and a f are as given in Eq. (6) and for down-type quarks they are given by Although the NLO massless contributions are known (see Refs. [16,26]), they are not relevant for our discussion of the quality of the approximations of the NLO top quarkloop contributions. They can be added easily to the final result since they only interfere with the exact LO expressions.
In terms of these F (0) i , we now define the finite, virtual contribution to the differential cross section (in analogy to gg → HH [44,45]), where C i is defined by In Eq. (41) α s corresponds to the five-flavour strong coupling constant. We introduce the quantity which is discussed in the following. For the renormalization scale we choose µ 2 = s/4.
When evaluating V fin we use exact expressions for the LO form factors and our highenergy expansions for the NLO parts. Exact results are known for the two-loop triangle form factors [22][23][24], however as shown in [35], the high-energy expansions reproduce the exact result almost down to the top threshold, which justifies our use of the expansions to evaluate these contributions also.
In Figs. 7 and 8 we show V fin as a function of √ s for different values of p T . In all cases we show Padé-improved results which are obtained by applying the method from Section 4 to V fin . We start with the numerical evaluation of the form factors f i , keeping the variable x as introduced in Section 4, perform the basis change to F i and then construct V fin as an expansion in x. At this point we apply the procedure outlined in Section 4. We take into account the different sets of Padé approximants listed below Eq. (38). Our best prediction, {10,16}, is shown as black points. For higher values of p T (Fig. 7) we show in addition two curves from the high-energy expansion.   For lower values of p T (Fig. 8) the high-energy expansion curves lie completely outside of the range of the plot axes. Nonetheless stable Padé predictions are observed, even for the low p T value of 150 GeV. For this value one observes that the higher orders in the m t expansion are crucial to obtain estimates with small uncertainties. The bottom-right panel of Fig. 8 shows the same approximations, but excludes the light quark contributions introduced in Eq. (39). This shows in more detail the improvement of the {10,16} approximation with respect to, for example, the {5,9} approximation.
The results presented in Figs. 7 and 8 make predictions for V fin which should eventually be confronted with numerical results as, e.g., announced in Ref. [46]. To this end we provide, in the ancillary files of this paper [29], a simple C++ program which interpolates a pre-evaluated grid of V fin points and uncertainties, evaluated with {10,16} for µ 2 = s/4, using routines from the GNU Scientific Library [47]. It can thus be used to reproduce the black points and uncertainty bars of Figs. 7 and 8.

Conclusions
We compute NLO QCD corrections to the process gg → ZZ induced by virtual top quark loops. We concentrate on the high-energy limit which corresponds to an expansion in the parameters m 2 t /s, m 2 t /t and m 2 t /u. We furthermore expand for small Z boson masses and show, at LO, that three expansion terms are sufficient to obtain per-mille accuracy. Analytic results, including terms up to order m 32 t m 4 Z , are presented for all 20 form factors in the ancillary file [29]. Additionally we include in this file the large-m t expansions of these 20 form factors, up to order 1/m 12 t . Using simple tensor structures as a starting point, we construct an orthogonal basis which is convenient when computing the squared amplitude. Alternatively we also provide LO and NLO results for the helicity amplitudes.
We extend the radius of convergence of the high-energy expansions with the help of Padé approximations. Our method provides both a central value and an uncertainty estimate. This is validated by comparisons to known exact results at LO. The Padé method is applied both to the form factors and the NLO virtual corrections to the differential cross section. In the latter case we include in our predictions also the LO contributions which originate from massless quark loops. In this setup the interference of the one-and two-loop top quark contributions amounts to about 5% as can be seen from Fig. 8.

D Example analytic results for LO expansions
In the ancillary file, we provide analytic expressions for the form factors defined in Eq. (12), as expansions in both the large-m t and high-energy limits. For the convenience of the reader we show here, for LO f 1 and f 3 , the leading terms of the expansions in a typeset form. In the large-m t limit they are given by