Top quark contribution to two-loop helicity amplitudes for $Z$ boson pair production in gluon fusion

We compute the top quark contribution to the two-loop amplitude for on-shell $Z$ boson pair production in gluon fusion, $gg \to ZZ$. Exact dependence on the top quark mass is retained. For each phase space point the integral reduction is performed numerically and the master integrals are evaluated using the auxiliary mass flow method, allowing fast computation of the amplitude with very high precision.


Introduction
Production of Z boson pairs is an important process at the LHC. The gluon fusion channel, gg → ZZ, is loop-induced. For this reason, it is suppressed by the strong coupling constant α s in comparison to the quark annihilation channel qq → ZZ which enters at tree level. However, the large gluon flux as well as event selection enhance the contribution of the gluon fusion channel to the hadronic cross section [1]. Therefore this production mode is essential for a reliable description of Z boson pair production.
The current status of the amplitude calculations for this process is as follows. The one-loop amplitude was calculated long ago [2,3]. The two-loop amplitude for massless internal quarks is also known [4,5]. However, until very recently, contributions of massive quarks have only been calculated approximately [6,7]. The goal of this paper is to present a calculation of the gg → ZZ two-loop amplitude keeping the dependence on the top quark mass. We note that when this paper was being completed, ref. [8] appeared where the two-loop amplitude with full top quark mass effects was calculated using analytic integral reduction and numerical integration using sector decomposition [9,10].
To compute the gg → ZZ amplitude, we largely follow the method described in our paper on a similar process, gg → W W [11]. The main difference lies in the integral reduction. For the present calculation, we decided to perform integration-by-part reductions individually for each phase space point. As before, the master integrals are evaluated efficiently using a system of ordinary differential equations and the auxiliary mass flow method [12,13].
The organisation of this paper is as follows. In Section 2 we describe the amplitude calculation, discuss the ultraviolet and infrared divergences that arise and summarise the renormalisation procedure. In Section 3 we discuss the evaluation of the master integrals. A benchmark point for the amplitude for physical kinematics and plots for the helicity amplitudes across partonic phase space are presented in Section 4. We conclude in Section 5.

Calculational setup
We consider the two-loop amplitude for Z boson production in gluon fusion This process is mediated by quark loops. We calculate the contribution where the external Z bosons couple directly to top quarks and disregard other quark flavours. The exact dependence on the quark mass m t is retained. We do not consider diagrams with an intermediate γ-, Z-or Higgs-boson. The photon mediated amplitude vanishes identically and calculations for the two latter are available in the literature [14][15][16][17][18]. All external particles are on-shell We define the usual Mandelstam variables which satisfy the relation s + t + u = 2m 2 Z . For four-dimensional external states the amplitude A can be decomposed in terms of 18 tensor structures Their definitions are given in ref. [1] and reproduced in appendix A for convenience. Our goal is to calculate the form factors A I=1,...,18 .

Pole structure
The form factors in eq. (2.4) contain ultraviolet (UV) and infrared (IR) divergences. In order to regulate them we employ dimensional regularisation. We work in space-time with dimension d = 4 − 2 .
Expanding the unrenormalised amplitude, A, in the bare strong coupling constant α s we write We employ the same renormalisation scheme as used in ref. [11] where the gluon field G µ and the top quark mass m t are computed in the on-shell scheme while the strong coupling constant, α s , is renormalised in the MS scheme. The relations between bare and renormalised quantities read where µ is the renormalisation scale and S = (4π) − e γ E . All renormalisation constants are expanded in the strong coupling The renormalised amplitude is related to the unrenormalised one by (2.10) The mass counterterm, C (1) , is calculated as a separate amplitude. The relevant renormalisation factors read [19][20][21][22][23] where γ g = 11 6 C A − 2 3 T F (n l + 1), n l is the number of massless quarks, and T F = 1 2 . After renormalisation the remaining poles are of IR origin. Their structure is predicted by the Catani formula [24]. Since gg → ZZ is loop-induced and has trivial colour structure, the IR poles are particularly simple. By subtracting the IR poles we define a finite remainder, (2.14) In the present case the Catani operator reads The -dependence of the leading order amplitude A (1) ( , µ) is important for computing the finite remainder.
From here on we disregard the trivial contribution proportional to n l , which in our calculation enters only through eqs. (2.11) and (2.15) and contributes a finite term proportional to the leading order amplitude.

Amplitude calculation and integral reduction
Using QGRAF [25] we generate 8 and 138 diagrams for the one-and two-loop amplitudes respectively. Colour and Dirac algebra is performed in FORM [26][27][28]. Symmetries between diagrams are established using REDUZE 2 [29] and we find that all two-loop diagrams can be mapped on to 21 integral families. The two-loop diagrams are classified according to the colour factors C A and C F . Factorisable two-loop diagrams are present in the t-and u channels only. A total of 37 two-loop diagrams vanish due to colour conservation, leaving 101 non-vanishing diagrams. Representative diagrams and their colour classifications are given in Table 1.
Colour factor # of non-vanishing diagrams Table 1: The classification of representative diagrams by colour factors. Curly lines represent gluons, wavy lines are Z bosons and straight lines represent top quarks. All nonplanar diagrams come with colour factor C A .
We decompose the two-loop amplitude according to the colour factors where c 1,2 are colour indices of the incoming gluons. From the discussion of the pole structure we note that A [ 2 ] , which is formed entirely of the factorisable diagrams, is finite. 1 A [C F ] is finite after UV renormalisation and hence IR poles are only found in The Z bosons couple to fermions through vector and axial currents. The gg → ZZ process involves two such vertices and can therefore be split into three parts; vector-vector, vector-axial, and axial-axial. The vector-axial part vanishes due to charge parity conservation and the two remaining parts can be considered separately. For each colour factor we have The weak couplings are given by where 4πα is the absolute value of the electron charge. We note that the vector-vector part of the factorisable diagrams, A vv , vanishes due to Furry's theorem.
In all non-factorisable diagrams both Z bosons are attached to the same quark loop leading to a single Dirac trace involving zero or two γ 5 matrices. These single traces can be evaluated naively in d dimensions by using the anti-commutative property of γ 5 . The factorisable diagrams contain traces involving a single γ 5 . For these diagrams we employ the Larin scheme [30], replacing the axial current by As this class of diagrams consists of products of one-loop anomaly diagrams, no finite renormalisation is required. After projection onto the tensor structures introduced in eq. (2.4), the form factors, A I , are in a form suitable for integration-by-parts (IBP) reduction. At this point we switch to an entirely numerical approach in which masses are chosen to be numbers close to their experimental values, 2 m t = 173 GeV, m Z = 91 GeV. (2.21) We perform numerical reductions to master integrals across partonic phase space using rational values for s and t. Having space-time dimension d as the only free parameter keeps all intermediate expressions compact and manageable. In addition to the computational advantage of working numerically, we also carefully choose the master integral basis to limit the size of the reduction tables and computation time. Even in the purely numerical case, significant simplification is achieved by avoiding reduction coefficients with denominators that are not factorised in kinematic invariants and space-time dimension. We find that by introducing numerator insertions and increasing denominator powers, non-factorisable denominators are eliminated.
The reductions to master integals required for the evaluation of the amplitude can be performed for a single phase space point in approximately 3 hours on a single CPU core using KIRA 2.0 [32]; speed strongly depends on fast read/write access to disk drives. Each reduction requires between 3 and 4 gigabytes of memory. The short runtime and low memory consumption enable straightforward parallelisation for a large number of phase space points.
The approach to integral reduction described here is in variance to our approach in previous work [11] where reductions were performed with s and t as free parameters. We kept the reductions tractable by projecting integrals onto one master integral at a time which, however, in turn made bookkeeping more involved. We find the strategy adopted in the present calculation with straightforward parallelisation and compact univariate reduction tables to be more efficient.

Numerical evaluation
Having expressed the full amplitude through master integrals, we need to evaluate them. The master integrals are defined as follows   To compute the 205 two-loop master integrals we follow the same procedure as described in ref. [11] and employ the auxiliary mass flow method [12,13]. We construct a system of differential equations with respect to m 2 t and solve it starting from the boundary conditions at m 2 t → −i∞ and moving to the physical value m t = 173 GeV. 3 Evaluating all 205 master integrals to 20 digits at a typical phase space point takes less than 1 hour on a single CPU core.
Comparing to the calculation of the gg → W W amplitude, we have more massive propagators in the gg → ZZ case, which leads to fewer regions at the boundary and simpler boundary conditions. Indeed, we find that the computation of the boundary conditions only requires 2 regions: 1. All internal momenta are comparable to m 2 t → −i∞.
2. Some internal momenta that form a closed loop are comparable to m 2 t → −i∞, while the remaining momenta are much smaller than m 2 t and are comparable to other kinematic parameters (i.e. s, t, m 2 Z ).
In figure 2 we show a typical master integral in these two regions. In each region we compute the large mass expansion of the integral. The expansion coefficients can be expressed in terms of the original integral with some of the propagators contracted, together with an overall power of m t . All of these contracted integrals can in turn be expressed by the two one-loop integrals listed in figure 3 through IBP reduction. Figure 2: A typical master integral and its leading regions. Solid and dashed lines describe massive and massless particles respectively. All internal massive particles have masses equal to m t , while all external massive particles have masses equal to m Z .
(a) I 1 (b) I 2 Figure 3: Master integrals for the boundary conditions. Solid and dashed lines correspond to massive and massless particles respectively. See appendix C for their explicit expressions [33].
We cross-check the evaluation of the master integrals against pySecDec [34,35] at a physical phase space point above the top quark threshold and find agreement to the default precision of pySecDec (3-10 digits).
In addition, we also check the self-consistency of the differential equations. Evaluations at two different phase space points should be connected by a system of differential equations with respect to s and t. We start from the phase point eq. (3.2), and move to the following two phase space points using differential equations in s and t: We compare the results with those obtained by directly solving the m 2 t differential equation at (s 1 , t 1 ) and (s 2 , t 2 ). We find that master integrals evaluated in the two different ways agree up to the precision used when solving the differential equations (20 digits).

Helicity amplitudes
Having discussed all the preliminary steps in the previous sections, we turn to the evaluation of the helicity amplitudes. We note that for each phase space point the numerical IBP reduction, evaluation of master integrals, and expansion of the form factors in require approximately 5 hours on a single CPU core.
The form factors are insensitive to the polarisation of external particles, which only enters through the tensor structures defined in eq. (2.4). The tensor structures are evaluated by constructing polarisation vectors for the external particles using spinor-helicity formalism. For the gluons the polarisation vectors are given by For on-shell production of Z bosons there is a total of 2 × 2 × 3 × 3 = 36 helicity amplitudes, but only 8 of them are independent [3,7,8]. For presenting our results, we only consider decays of the Z bosons to two massless fermions. In this case, the polarisation vectors of the Z bosons can be represented by the following currents in the spinor-helicity formalism, * µ  The spinors represent massless fermions in the decays p 3 → p 5 + p 6 and p 4 → p 7 + p 8 . Note that we omit propagators and couplings related to the decays. The equivalence of conjugation and interchange of leptons 5 ↔ 6 and 7 ↔ 8 further reduces the number of helicity configurations to two. We label the two independent helicity amplitudes according to the helicities of particles 1, 2, 5, 7: LLLL and LRLL.
In table 2 we present renormalised helicity amplitudes for the phase space point given in eq. (4.4). Extrapolating the precision loss order by order in in table 2 we find a precision of around 10 digits for the finite part. This precision is expected to be attainable for the bulk of partonic phase space. Closer to the production threshold ( √ s = 182 GeV) we observe strong cancellation and precision drops as a result thereof. Since the precision of the helicity amplitudes is governed by the master integral evaluations, it can be increased to the desired level by reevaluating the master integrals with more digits.
The partonic phase space can be parametrised using the relative velocity of the Z bosons, β, and the scattering angle between p 1 and p 3 , θ, in the centre-of-mass frame. The the kinematic invariants are related to β and θ as follows, Plots for the interference between the finite remainder of the two-loop amplitude eq. (2.14) and the leading order amplitude are given in table 3. We define the interference as and present the plots separately for each colour factor. Note that while we use integer values for the particle masses in the amplitudes A aa , see eq. (2.18), the weak mixing angle is evaluated using physical Z and W masses [31], m Z = 91.1876 GeV, m W = 80.379 GeV. (4.8) We have checked that the form factors satisfy required (crossing) symmetry relations from refs. [5,7]. We have also checked that we obtain the same finite remainders by using the four-dimensional tensor structures in eq. (2.4) and the d-dimensional tensors used in refs. [4,5,7,8]. Using four-dimensional tensors makes intermediate expressions simpler and requires fewer terms in the -expansion of the master integrals.
We compared the results of our calculation with those of refs. [7,8] and found good agreement. For the result presented in ref. [7], we compared the helicity amplitudes at two high-energy points ( √ s = 1290 GeV, cos θ = ±0.2). 4 We cross-checked d-dimensional form factors against the points provided in ref. [8], including finite terms kindly provided by the authors.

Conclusions
We described a numerical calculation of the top quark contribution to the two-loop helicity amplitudes of the process gg → ZZ. The dependence of the results on the top quark mass is accounted for exactly. We have presented results for a single benchmark point as well as plots for the finite remainders of the helicity amplitudes across the phase space. The plots are produced from a grid of amplitude evaluations. For each phase space point integrationby-parts reductions were performed numerically for integer mass values and the master integrals were evaluated using the auxiliary mass flow method. This allows for efficient evaluation of the amplitude to, essentially, arbitrary numerical precision.  Table 2: Evaluation of the renormalised two-loop helicity amplitudes, see eq. (2.10), for the phase space point defined in eq. (4.4). We normalise by the one-loop amplitude A (1) | =0 and show the infrared pole structure for comparison, see eq. (2.14). The renormalisation scale is set to µ = m Z = 91 GeV.
Amplitude LLLL LRLL Table 3: Interference of the finite remainders of the two-loop amplitude and the leading order amplitude as functions of β and cos θ. See eqs. (4.6) and (4.7) for details. The renormalisation scale is set to µ = m Z = 91 GeV.

C Boundary condition of the differential equation
The explicit expressions for the boundary integrals in figure 3 are listed below [33], where q 2 corresponds to the four-momentum squared of the external legs. Note that these boundary integrals are one-loop integrals, thus they enter the boundary condition through products among themselves.