Two-loop helicity amplitudes for $gg \to ZZ$ with full top-quark mass effects

We calculate the two-loop QCD corrections to $gg \to ZZ$ involving a closed top-quark loop. We present a new method to systematically construct linear combinations of Feynman integrals with a convergent parametric representation, where we also allow for irreducible numerators, higher powers of propagators, dimensionally shifted integrals, and subsector integrals. The amplitude is expressed in terms of such finite integrals by employing syzygies derived with linear algebra and finite field techniques. Evaluating the amplitude using numerical integration, we find agreement with previous expansions in asymptotic limits and provide ab initio results also for intermediate partonic energies and non-central scattering at higher energies.


Introduction
Z boson pair production is an essential process at the Large Hadron Collider (LHC). Besides its immediate relevance as a signal process for precision physics [1][2][3][4][5], it is a significant background to on-shell and off-shell Higgs production for the four-lepton final state [6][7][8][9]. Continuum Z pair production significantly contributes to off-shell Higgs production (∼ 10%) through interference effects [10,11]. This is, in particular, important for indirect Higgs width constraints as proposed in [12,13]. The primary production channel for vector bosons at the LHC is quark-antiquark annihilation, which starts at tree level and is known to next-to-next-to-leading order (NNLO) QCD [14][15][16][17][18][19][20]. The gluon fusion channel is loopinduced and starts formally at NNLO for the process pp → ZZ. Nevertheless, it accounts 2 Setup of the calculation

Form factors and helicity amplitudes
We consider Z pair production in gluon fusion, g(p 1 ) + g(p 2 ) − → Z(p 3 ) + Z(p 4 ) . (2.1) Here, p 1 , p 2 are incoming and p 3 , p 4 are outgoing momenta, so that p 1 + p 2 = p 3 + p 4 and that is, we consider the Z-bosons to be on-shell. Our Mandelstam variables are 3) The amplitude can be represented as using polarization vectors λ i (p i ), for which we will also use the abbreviation i ≡ λ i (p i ).
Using Lorentz invariance, the amplitude can be decomposed in terms of 138 parity-even tensor structures [26]: M µνρσ (p 1 , p 2 , p 3 , p 4 ) = a 1 g µν g ρσ + a 2 g µρ g νσ + a 3 g µσ g νρ ( a 1,ij g µν p ρ i p σ j + a 2,ij g µρ p ν i p σ j + a 3,ij g µσ p ν i p ρ j + a 4,ij g νρ p µ i p σ j + a 5,ij g νσ p µ i p ρ j + a 6,ij g ρσ p µ i p ν j ) a ijkl p µ i p ν j p ρ k p σ l . (2.5) Parity-odd tensor structures involving the epsilon tensor do not need to be taken into account due to Bose symmetry and charge-parity conservation for our process [24]. Since the color structure of the external states is straight-forward, we suppress color indices here and in the following. We can reduce the number of tensors using transversality of the gluon polarization vectors, where L Z ff = (I f 3 − q f sin 2 θ W )/(sin θ W cos θ W ), R Z ff = −q f sin θ W / cos θ W , e is the positron charge, and q f is the electric charge of the fermion in terms of e. The vector and axial components are given in terms of the weak mixing angle θ W by v t = 1 2 − 4 3 sin 2 θ W and a t = − 1 2 , respectively. The couplings of the two Z bosons to the fermion line can in principle generate vector-vector (v 2 t ), vector-axial (v t a t ), and axial-axial (a 2 t ) contributions to the amplitude. However, due to Bose symmetry and charge-parity conservation for this process, the vector-axial contribution should vanish identically [24]. This also explains the absence of any terms with the Levi-Civita tensor in (2.5) since such terms would violate parity and are hence forbidden. For the massless quark case, the vector-vector and the axial-axial contributions are identical; after including quark masses, they differ by terms proportional to the quark mass.
We find a total of 166 diagrams containing at least one top-quark propagator. Out of these, 49 diagrams have a single gluon coupled to a closed fermion loop, and hence they vanish due to colour conservation. The remaining diagrams can be divided into two classes shown in figure 1.
Class A : Both Z bosons couple to the same fermion line. To appropriately handle γ 5 in d dimensions, we use the anti-commuting γ 5 scheme described in [43,44]. Since cyclicity of trace is not preserved in this scheme, a reading point prescription is employed to ensure that all traces are read from the same point. However, for a closed fermion loop with an even number of γ 5 matrices, it is trivial to eliminate γ 5 using the anti-commutation relations; this greatly simplifies the implementation of the anti-commuting γ 5 scheme.
Class B : The Z bosons couple to different closed fermion lines. For these diagrams, the vector-vector contribution can be shown to vanish due to Furry's theorem, while the vectoraxial piece is identically zero because of charge-parity conservation. The axial-axial piece, however, vanishes only after summing over a degenerate SU (2) L doublet. Since the third generation of quarks is not degenerate, this cancellation is incomplete and we see a finite remainder from the top-bottom mass splitting. These diagrams have a single γ 5 in each loop which leads to a non-trivial structure and requires careful application of the reading point prescription. Since these diagrams are effectively one-loop, we treat them separately. Exact results for these diagrams were previously presented in [34] and we find full agreement. We will not mention them any further and do not include these contributions in the results presented below.
After generating the diagrams in class A, we employ Reduze 2 to map them to the 4 different integral families shown in table 1. We find 13 top-level topologies (trivalent graphs), of which 7 are irreducible (figure 2) and 6 are reducible (figure 3). We use FORM   [45][46][47] to apply the Feynman rules and generate the amplitude, employing the 't Hooft-Feynman gauge (ξ = 1) for internal gluons. Before applying any symmetries, we find a total of 29247 integrals with up to 4 irreducible scalar products in the numerator. This number can be reduced to 4504 using symmetry relations between the integrals. We see further simplifications after inserting the symmetry relations in the amplitude; due to cancellations only 1584 integrals survive in the form factors. This is a significant improvement over the original number of integrals and underlines the importance of using symmetry relations and working with amplitudes instead of individual diagrams.
Using Reduze 2 [48][49][50][51], we perform a numerical reduction by substituting numbers for kinematics and find, for the diagrams with a single fermion loop, 85 irreducible topologies with the worst sector having 6 master integrals with 6 lines. In total, we obtain 264 master integrals for class A, out of which 172 are not related by any crossing. Our symbolic reduction is discussed in the following section.

Linear relations from syzygies
A general L-loop scalar Feynman integral with N propagators can be represented by where k 1 , ..., k L are the loop momenta, q i are the propogator momenta (linear combinations of loop and external momenta), m i are the masses of the propagators, ν i are (integer) exponents of the propagators, and d = 4−2 . Here, we allow also for non-positive powers ν i of the propagators, i.e. we consider a family of integrals with possible irreducible numerators. The total derivative of an integral in dimensional regularisation vanishes; this allows us to write linear relations between different integrals [52] where v µ could be any linear combination of loop and external momenta. We can eliminate most of the integrals in the amplitude using these relations with the remaining integrals usually referred to as master or basis integrals. This procedure can be systematically used to reduce any integral appearing in the amplitude due to an algorithm by S. Laporta [53]. Many public codes based on this algorithm are available for this purpose [48,[54][55][56][57]. Conventionally, the vector v µ is chosen as a single loop or external momentum; different such choices yield a set of simple equations as a starting point. It is easy to see that the derivatives in (3.2) generate higher powers ν i of the propagators, often referred to as "dots". Such auxiliary integrals with a large number of dots are usually not required for the amplitude and lead to relatively large linear systems that are computationally expensive to reduce.
A method was proposed in [58] to avoid these higher powers of propagators by constructing suitable generating vectors v µ from syzygies. This method involves the computation of a Gröbner basis to obtain the syzygies. A linear algebra based approach was presented in [59], albeit the syzygies can only be obtained to a specified degree using this method. Subsequent work [60][61][62] refined syzygy based constructions in the momentum space representation as well as in Baikov's representation [63]. Syzygies can also be used to derive linear relations [64][65][66] in the Lee-Pomeransky representation [67].
The L-loop Feynman integral in (3.1) can be written in Baikov's representation as where the Jacobian of the variable transformation involves the determinant P , the Baikov polynomial, N 0 is a normalization factor, and E is the number of linearly independent external momenta. The integration-by-parts identities in Baikov's representation are given by where f 1 , . . . , f N are arbitrary polynomials in the Baikov parameters z 1 , . . . , z N , and the kinematic invariants. In the above equation, terms that appear with 1/P lead to dimensionally shifted integrals. Since these integrals don't appear in the amplitude, it may be desirable to avoid them to prevent an unnecessary proliferation of auxiliary quantities in the system. This can be achieved by imposing the constraint Here, we introduced a new polynomial f N +1 in the Baikov parameters. Note that P and its derivatives are known polynomials for the problem, see e.g. [68] for details. A constraint of this type on the vector of polynomials (f 1 , . . . , f N +1 ) is known as a syzygy in algebraic geometry. Explicit solutions to this equation were pointed out in [68] and can easily be written down. The resulting f i are linear polynomials in the Baikov variables z k and the kinematic invariants. It must be noted that these f i generate integration-by-parts relations which cover [68] those derived in the conventional momentum-space approach (3.2).
To enforce the absence of doubled propagators, one requires that for all i with ν i ≥ 1, the f i are proportional to z i to cancel the 1/z i in the relation, While it is straight-forward to fulfil both constraints (3.5) and (3.6) separately, a simultaneous solution requires a non-trivial calculation.

Constructing syzygies with linear algebra
Formally, finding vectors of polynomials (f i ) which are simultaneous solutions of both (3.5) and (3.6) corresponds to the determination of the intersection of two syzygy modules [69]. In practice, computer algebra packages implement algorithms to solve this task. For performance reasons we decided to develop a custom syzygy solver based on linear algebra and finite field arithmetic [70,71]. Note that if polynomials (f i ) satisfy the syzygy constraint in (3.5), then (z k f i ) for any k also satisfy it.
In algorithm 1, we provide a description of our method which converts the intersection problem up to a specific degree of the syzygies to row reduction of a matrix. Here, we treat the kinematic invariants as indeterminates of the polynomial ring, such that the matrices M n have entries which are rational numbers. Alternatively, one can treat the invariants as part of the coefficient field. This decreases the number of columns of the matrices M n , Algorithm 1 Syzygies for linear relations without dimension shifts or dots Input: Syzygies of degree 1 solving (3.5), maximal required degree n max . Output: Syzygies S 1 , . . . , S nmax up to degree n max solving (3.5) and (3.6).
1: Start with syzygies of degree n = 1. Let I 1 be a complete set of solutions (f i ) to the no-dimension-shift constraint (3.5), which are linear in the Baikov parameters z k . These can directly be written down [68]. Abbreviating the momenta squared with variables z N +1 , . . . , the vectors in I 1 are of homogeneous degree 1 in the variables z k .
2: At degree n, form a matrix M n , where each element of (f i ) ∈ I n corresponds to a row. The columns enumerate both the component i of (f i ) and the power products of z i in them; the entries of the matrix are the coefficients. A column is called admissable, if it satisfies the no-doubled propagator constraint (3.6), and non-admissable otherwise. All admissable columns are ordered to the right of the non-admissable columns.
3: Perform a row reduction of M n . In the row reduced form, select all rows, which have an admissable pivot column and form the corresponding syzygies S n from them. S n forms a complete set of linear combinations of the syzygies in I n , which satisfy (3.6) for all of their terms, and are therefore our solutions at degree n.

4:
If n is the user-defined maximal degree, stop and return the solutions S 1 , . . . , S n . Otherwise, proceed.

5:
For each vector of polynomials (f i ) ∈ I n and each z k , form the vector of polynomials (z k f i ). This gives the set I n+1 , which are solutions of (3.5) of degree n + 1 in the z k but not necessarily solutions of (3.6).
6: Replace n → n + 1 and go to step 2.
but the entries are then rational functions of the kinematic invariants. Since in the second approach the kinematic invariants do not count towards the degree n in our algorithm, a lower maximal value may be sufficient for the integral reduction problem at hand compared to the first approach. It is useful to use e.g. an overall mass dimension squared as a homogenizing variable z N +1 for the last component of the syzygy vectors in this setup.
The row reduction of the matrix M n eliminates redundancies between the syzygies at degree n. In our approach, we generate templates for the generation of linear relations between Feynman integrals from the syzygies. We allow the templates to be applied to seed integrals with specific integer propagator powers and perform a subsequent row reduction on the resulting identities, similar to the traditional Laporta algorithm. In this approach, we find it useful to filter out syzygies that are just a lower degree syzygy multiplied with an overall power product in the z k . This is achieved by determining reducible monomials using the row reduced form of an auxiliary matrix for the syzygies induced by lower degree syzygies [72].
For our current process, we generated the required syzygies and performed the subsequent Laporta step with an in-house linear solver, Finred, based on finite field arithmetic and rational reconstruction. To simplify the linear relations further, we set m t = 1 and use a numerical value for the Z-boson mass as a ratio over top-quark mass, m 2 Z /m 2 t = 5/18. This amounts to factoring out powers of m 2 t corresponding to the mass dimension of the respective form factor. In this way, we successfully reduced all of the Feynman integrals in our calculation to master integrals. The reductions proved to be rather challenging nevertheless and required significant computational resources. This is evident from the fact that the reduction tables exceeded 200 GB in size, with rational functions of degrees of up to 190 in the kinematic variables appearing in the reduction tables. The non-planar topologies, unsurprisingly, were the most difficult and accounted for almost all of the computation time and disk space. An interesting point to note is that within the planar topologies, figures 2a, 2c, and 2b, with adjacent gluons are significantly simpler than 2d with the gluons at the opposite vertices.

Inserting the reductions into the amplitude
After having generated the reduction identities, the next task is to insert them into the unreduced amplitude. The reduction identities for this process are very complicated with a size of over 200 GB. As such, this task in itself is a major challenge. We used several tools and techniques to make this more manageable. We first calculate the reduction identities to the conventional Laporta basis and perform multivariate partial-fractioning of the reduction tables based on polynomial reductions with respect to a Gröbner basis [73,74]. We implement this using the public code Singular [75] 1 with a polynomial ordering that prefers polynomials with lower degrees in kinematic variables and smaller coefficients, and are able to drastically reduce the size of the reduction tables. We found it useful to first perform the partial fractioning for the d dependent denominators and then partial fraction the kinematic denominators. We note that this procedure can also be used even in the presence of denominators which depend both on d and the kinematic variables.
We use custom FORM scripts to insert the reduction identities into the amplitudes, and again perform multivariate partial-fractioning on the reduced amplitudes to arrive at a simpler representation in terms of kinematic variables and the irreducible denominators. We see a drastic level of compression at this step; after partial-fractioning, the total size of amplitudes reduces from ∼300 GB to ∼600 MB.
Next, we perform a change of basis to express our amplitudes in terms of finite integrals. We explain this choice of basis in more detail in section 4. After partial fractioning the basis change identities, we insert them into the reduced amplitude to arrive at the final reduced amplitude in terms of our finite basis. This step was computationally expensive, with more than a week of run-time and intermediate expressions with sizes in the terabytes before partial fractioning. Once we have the amplitudes in the new choice of basis integrals, we perform partial-fractioning to simplify them. Note that the form factors in the conventional Laporta basis contain many denominator factors that are polynomials in both d and kinematics; we find that all such denominators no longer appear for our choice of finite master integrals.
As a last simplification measure we expand the form factors around d = 4. Our projectors introduce a spurious pole of order 1/ 5 , which cancels after reduction. Since we are calculating an NLO amplitude, UV and IR subtractions will involve at most 1/ and 1/ 2 poles, respectively. The reduced bare form factors should therefore not have any pole worse than 1/ 2 , which, however, is not completely manifest when using our symbolic master integrals. However, the change of basis to finite integrals removes the 1/ 4 poles at the algebraic level. In section 5 we describe in detail how all the poles show the expected behaviour with high numerical precision.
In the final representation, we are able to bring down the size of the worst coefficients to less than 1 MB. We create a C++ library for fast evaluation of the integral coefficients, either with exact rational arithmetic or with arbitrary precision floating point arithmetic using the GMP library. Even though the expressions are still sizable, we can evaluate all coefficients for a generic point in phase space within half a minute using rational arithmetic or within 3s using floating point arithmetic with a target precision of 15 digits on a single CPU core.

Dimension shifts and dots
To evaluate the master integrals, a powerful approach is to use differential equations to find analytic solutions [77][78][79][80][81][82]. This approach was used to calculate the master integrals in terms of multiple polylogarithms for the 2-loop massless corrections to diboson production in [16,18,83,84]. Due to the massive top-quark loop in the corrections considered here, we expect the presence of functions beyond multiple polylogarithms, which makes the evaluations of master integrals considerably more challenging. While there has been significant progress concerning the analytic evaluation of Feynman integrals beyond polylogarithms [85][86][87][88][89][90], integrals of the type considered here remain a challenge. An alternative is the use of expansions to solve the differential equations numerically [41,[91][92][93][94]. Here, we use a purely numerical approach to integrate the master integrals, namely sector decomposition [95][96][97][98]; see also [99,100] for recent applications.
A naive integration-by-parts reduction using Laporta's algorithm with a generic ordering criterion leads to a conventional basis of master integrals. This basis is rather difficult to evaluate numerically since the integrals are often divergent and numerically unstable, and as such is inadequate for our purpose. We instead choose a different basis of master integrals that is finite in the limit d → 4. It was observed in [99,101] that using a basis of finite integrals is highly beneficial, leading to a numerically more stable behaviour. Additionally, finite integrals often require fewer orders in the expansion which, coupled with better numerical stability, improves the overall performance significantly.
One possible approach to constructing finite integrals is to use dimensionally shifted integrals [102], possibly with doubled (or higher powers of) propagators. It is always possible to construct a basis in this way [103,104] and also straightforward in practice using e.g. the finite integral finder in Reduze 2. Examples of such integrals are shown in figure  4. While it is convenient to find such finite integrals with dimension shifts and dots, they require computation of additional reduction identities beyond those required for the amplitude. For example, reductions for integrals with 2 additional dots are required for the dimension shift of two-loop integrals, and typically such integrals do not directly appear in the amplitude calculation. It may therefore seem interesting to consider alternative choices of finite integrals. Here, we explore a different approach by constructing finite integrals through linear combinations of divergent integrals based on the Feynman parametric representation [105]. In such linear combinations, non-integrable divergences of individual integrals cancel at the integrand level. This results in a single generalised Feynman parameter integral that is finite. We briefly describe the algorithm in the following subsection.

Constructing finite linear combinations
Consider a general L-loop integral in d dimensions with N distinct propagators in the momentum space representation, with integer exponents ν j ∈ Z. If all indices ν j are positive, one can use (see e.g. [106,107]) to derive the Feynman parametric representation of this integral, We can include inverse propagators (numerators) with ν j < 0 by employing the identity [108, 109] Let N + be the set of all positive ν j , N − the set of all negative ν j , and r = j∈N + ν j . Then, an integral with positive or negative indices can be written as Our goal is to combine different integrals sharing a common parent topology into one merged parametric representation. We therefore wish to base our Feynman parametric integral on the resulting U and F polynomials for the parent sector. For integrals belonging to subtopologies of the parent sector, this can be achieved by taking derivatives with respect to the Feynman parameters corresponding to the pinched lines without setting them to zero, for ν j ≤ 0.
Note that we allow the pinched lines to appear as numerators i.e. ν j ≤ 0 for j ∈ N ∆t . The Symanzik polynomials U and F are calculated by taking all indices N into account.
With the prerequisites in place, we can now formulate algorithm 2 to construct linear combinations of integrals, which have a convergent Feynman parametric representation for = 0.

Algorithm 2 Finite Feynman integrals
Input: Dimensionally regularized multiloop integrals with a common parent sector, possibly involving higher powers of propagators, irreducible numerators, or dimension shifts. Output: Linear combinations of the input integrals which are finite, i.e. they have a convergent Feynman parametric representation for = 0.
1: From the n s input or "seed" integrals, form a general linear combination where I i are the seed integrals and a i are the unknown coefficients. The a i are assumed to depend on the kinematic invariants and the dimensional regulator .
2: Using (4.7), write the Feynman parametric representation for each seed integral and bring their linear combination over a common denominator such that where N T is the set of distinct propagators in the parent sector, ν 0 is the effective number of propagators, and d 0 ∈ Z the effective number of space-time dimensions to be expanded around. The numerator P is a homogeneous polynomial in the Feynman parameters, where the coefficients c j are polynomials in a i , the kinematic variables, and , and M j (x 1 , ..., x Np ) are monomials in Feynman parameters. Note that the numerator polynomial P in general depends on and it is crucial to keep this dependence to produce correct results. It is sufficient, however, to set = 0 in the exponents of the U and F polynomials for the convergence analysis in the following two steps.
3: Check the scaling behaviour of the integrand near an integration boundary using the prescription outlined in [104,110]. 4: Make sure a convergent integration of (4.9) is not prevented by a rapid growth of the integrand near the boundary. This can be achieved by requiring the coefficients of the offending monomials in the numerator to vanish, which provides constraints on the a i .   Figure 6. Integrals appearing in (4.12). I 1,2 is the corner integral of the topology under consideration. I 2,2 is a second integral in the topology, but with an extra numerator (k 2 − m 2 t ) where k is equal to the difference of the momenta of the edges marked by the thin dashed lines. Integrals I 3,2 , I 4,2 , I 5,2 , I 6,2 , I 7,2 belong to subtopologies. All integrals are defined in d = 4 − 2 dimensions.
As an example, we applied our algorithm to a set of seed integrals including those shown in figure 5 and obtained the finite linear combination (4.11) Allowing for seed integrals with higher numerator rank the algorithm finds, amongst others, the finite linear combination  very similar. In fact, it is straightforward to see that for any finite linear combination, an additional numerator can be added while keeping the integral IR finite. Through power counting one can see, that the additional numerator does not introduce a UV divergence in our present example. However, linear combinations obtained simply by augmenting the existing integrals with additional numerators aren't the only possibilities at higher numerator rank. Indeed, we observe that generally the number of finite linear combinations increases with the numerator rank.

Numerical performance
One can try to express the amplitude in terms of finite linear combinations, which are defined in 4 − 2 dimensions and have at most additional numerators. In practice, we found it useful to consider integrals with "dots" and dimension-shifts as well, primarily for the following reasons: • It can happen that already the corner integral of a sector has a UV divergence, which can not be cured by a subsector subtraction. Obviously, a numerator insertion is not going to help. One could try to use a supersector instead, but this can have other disadvantages such as an unnecessary increase of analytic complexity.
• Choosing integrals with higher numerator ranks leads to extreme proliferation in the number of terms in the numerator polynomial, often leading to rather large pySecDec libraries that are difficult to compile on GPUs. Our efforts to condense the numerators to a more manageable size resulted in the appearance of spurious poles that often worsened numerical stability.
• In a slightly different approach, integrals with both numerators and dots can be combined to form finite combinations. These integrals, however, have higher powers of the F polynomial in the denominator. In our experiments, this led to significantly worse numerical performance in the physical region, where contour deformation is required.
A comparison of numerical performance for different divergent and finite integrals for the first few orders in expansion is shown in table 2. It is clear that the finite integrals perform significantly better. The finite integral in figure 4c has the lowest exponent for 1/F, and unsurprisingly shows the best numerical performance. We can also see that the finite linear combination in (4.11) is on par with the dimension-shifted integrals, which demonstrates its viability. One interesting point to note is that both linear combinations have integrands with 1/F 3 compared to 1/F 2 for the dimension-shifted finite integral in figure 4d while having similar performance.
We observe the best numerical performance for a combination of both approaches: finite linear combinations and dimension-shifted integrals. In addition, we choose our finite basis of master integrals so that the d-dependence of the denominators appearing in the reduction identities factors out, using the code of [113] (see also [114]). In other words, there are no irreducible denominator factors that are polynomials in both kinematics and d for this choice of master integrals. The definitions of the finite master integrals used in our calculation in terms of divergent integrals are provided in an ancillary file.

UV renormalisation and IR subtraction
We expand the bare form factors A i perturbatively according to where α s,0 is the bare QCD coupling. Since the LO process already starts at one loop, the two-loop process is effectively an NLO correction. We first perform UV renormalisation of α s in the 5-flavour MS scheme, n f = 5, with the top-quark contribution to the gluon self energy subtracted at zero momentum [115] using where S = (4π) e −γ E , γ E ≈ 0.577 is Euler's constant, µ R is the renormalisation scale, and µ 0 is the 't Hooft scale introduced in the dimensionally regularized bare amplitude. The renormalisation constant Z αs is given by where We renormalise the top-quark mass in the on-shell scheme. The renormalised top-quark mass is related to the bare mass according to In practice, we find it convenient to account for the top-quark mass renormalisation by inserting counterterm vertices in the 1-loop diagrams. Finally, we take into acount the gluon wave function renormalisation by multiplying the amplitude with Z 1/2 G for each external gluon, where the gluon renormalisation constant is defined as This gives us the renormalised form factors The IR structure of NLO amplitudes was first predicted by Catani in [116]. Here, we perform IR subtraction using the "q T scheme" described in [117] with where δ We present all of our results for µ 2 R = s.

Checks
We perform the following checks to establish the correctness of our results: (i) We verify our 1-loop amplitude against the literature, specifically the form factors provided in [36]. This is essential to make sure we match conventions and to facilitate comparisons for the 2-loop result.
(ii) We explicitly check that the form factors satisfy the identities in (2.13). We see an exact algebraic identity at the level of reduced amplitude with symbolic kinematics.
(iii) We also verify, numerically for a phase space point, that the relations in (2.14) are satisfied.
(iv) We check all the finite integrals by numerically evaluating them and comparing them against their explicit definitions in terms of divergent integrals for a phase space point.
(v) We observe algebraic pole cancellation for the leading poles, see section 3.3 for details. We calculate our amplitude for a Euclidean point, using Reduze 2 to generate numerical reductions for the Euclidean point. We verify that spurious 1/ 3 poles vanish after integration (15 digits), and that the 1/ 2 and 1/ poles match Catani's IR formula [116] (9 digits for the double pole, 7 digits for the single pole) as shown in the first table of appendix A.
(vi) We verify that the poles match Catani's IR formula [116] for a point in the physical region as shown in the second table of appendix A. For each phase-space point we compute, we also automatically check that the poles match the IR formula within our numerical uncertainty.
(vii) We evaluate the amplitude using an alternate finite basis and compare it against our result from the primary basis. This acts as a very strong check of our calculation since it validates the basis change, the definitions of our finite integrals, and the reliability of their numerical evaluation. We find agreement between the two bases within the expected numerical error, typically within a few percent for the form factors. It must be noted that this alternate basis is numerically a lot less stable and unsuitable for large scale evaluation runs.
(viii) We compare the axial-axial piece of the amplitude evaluated using Kreimer's anticommuting γ 5 scheme [43,44] with a separate amplitude calculation utilising Larin's γ 5 scheme [118,119]. For the latter calculation we avoid the appearance of γ 5 by expressing all axial-currents in terms of Levi-Civita symbols. Metric tensors obtained from contracting two Levi-Civita symbols are treated as d-dimensional. Finally, a finite renormalisation is applied for each non-singlet axial current as required to restore the Ward identities. We emphasize that a verbatim application of the scheme as described in [119] is motivated (ignoring e.g. higher order terms in the symmetry restoration constant), because of the finiteness of our one-loop amplitudes. Performing two independent amplitude calculations utilising different schemes for the treatment of γ 5 provides a strong check of our amplitude calculation. We find agreement between the two calculations for a physical phase space point within numerical precision.
(ix) We check that our result reproduces the large top-mass expansion [32][33][34] below the top-quark threshold and the small top-mass expansion [36] above; a detailed comparison is presented in the next section.

Results
Here, we present the results of our calculation and compare them against several approximations available in the literature. In particular, we perform comparisons against the large top-mass expansion [32][33][34] as well as the small top-mass power series and Padé expansions [36]. Let us define the quantities relevant for presentation of our results. We work in the helicity basis defined by (2.16). Concretely, we can write the UV renormalised and IR subtracted helicity amplitudes with incoming helicities λ 1 , λ 2 and outgoing helicities λ 3 , λ 4 as The amplitudes are expanded as To prepare the summation over polarisations, we consider contributions to the squared 1loop helicity amplitudes V (1) , and to the interference between the 1-loop and 2-loop helicity amplitudes V (2) defined as Note that for our numerical results, we include only the pure top-quark contributions of class A computed here, both in the amplitudes and in the interference terms. There are 36 helicity amplitudes in total for the gg → ZZ process, fulfilling various relations; see section 2. In order to further condense the presentation of our results, we average over the helicities of the incoming gluons and define the quantities where λ 1 , λ 2 ∈ {+, −}, and λ 3 , λ 4 ∈ {+, −, 0}. In the results shown, we choose our electroweak couplings as where the Fermi constant G F and Z boson mass m Z are fixed according to [120]. In our calculation, we fix m 2 Z /m 2 t = 5/18; inserting the value of m Z from (6.6) implies that m t = 173.016 GeV and m W = 80.296 GeV. The weak mixing angle is fixed according to sin(θ W ) = 1 − m 2 W /m 2 Z . Note, however, that the only mass ratio fixed in the computationally expensive part of our calculation is that of m 2 Z /m 2 t ; all other mass values and couplings can straightforwardly be varied in our code. All results are presented at renormalisation scale µ 2 = s.
For numerical evaluation of the integrals appearing in our amplitude we apply sector decomposition and integrate using the quasi-Monte Carlo (QMC) algorithm first applied to sector decomposed Feynman integrals in [112], as implemented in the program pySecDec [97,111]. For a review of QMC methods from a mathematical perspective see, for example, [121]. We separately evaluate terms appearing in the form factors of our amplitude according to their colour factor (C F or C A ) and whether they form part of the vector-vector (v 2 t ) or axialaxial (a 2 t ) contribution. For each phase-space point we aim to obtain percent level or better precision for each of the A i form factors, for each colour structure and for the vector-vector and axial-axial pieces separately. To present our results, we then rotate to the helicity basis defined in section 2. In order to improve the efficiency of this approach, the target precision of each integral is set according to its contribution to the uncertainty on the form factors using a variant of the algorithm presented in [99]. For most phase-space points, the time required to obtain this precision varies between 90 minutes and 24 hours on 2 Nvidia Tesla V100 GPUs. This time is completely dominated by the numerical integration of the master integrals; the time to evaluate the coefficients is basically negligible in this context, see section 3. The result of each integral for a given phase-space point is shared between all form factors, colour structures and vector/axial pieces. We observe that requiring percent level precision on all of the form factors individually typically results in most of them being obtained to per mille or better precision. The resulting precision obtained for the interference terms V (2) λ 3 λ 4 is per mille or better as well. We expect that a further performance improvement can be achieved by optimising the sampling of the integrals according only to their contribution to the numerical error of the interferences rather than the individual unphysical form factors. Table 3 shows our numerical results for the independent helicity amplitudes for a physical point in phase space. The same phase space point is used for the second table in Appendix A, where also the corresponding (γ 5 scheme dependent) values for the form factors A 1 , . . . , A 20 are shown. We would like to emphasize that all of the plots below actually show error bars for our numerical results; however, the errors are too small to be visible in the plots. top mass expansion, we see that the Padé approximation provides a drastic improvement with respect to the power series approach: while the power series data is visible within the plotted range only for the two highest energies and diverges substantially from the exact result for smaller values of √ s, the Padé approximation agrees very well with our result down to much lower energies. We note that the quantitative level of agreement between the exact results and the approximations depends greatly on the details of the scheme in which the comparison is performed. For example, to convert the finite 2-loop interference term V (2) from the "q T scheme" used in this article (see section 5.1) to Catani's original convention [116] which has also been used in eq. 13 of [36], at scale µ 2 R = s we must subtract π 2 C A V (1) (for the real part). The size of the shift resulting from the difference between the two subtraction schemes is comparable to the 2-loop interference terms themselves and can therefore significantly alter the shape of the corrections, their overall size and the level of agreement with the expansions. Furthermore, the 2-loop curves presented in the following roughly follow the form of the 1-loop curves in the "q T scheme", while this is generally not the case in Catani's original scheme. We give explicit examples for these effects in the appendix B.
In figure 8 we show a comparison of our calculation to the expansions as a function of cos(θ) for a fixed value of √ s. It is apparent that the large top-mass expansion is very stable with respect to variation in cos(θ) (Top Left Panel). The small top-mass power series expansion, on the contrary, diverges rapidly away from cos(θ) = 0. The breakdown of the small top-mass approximation away from cos(θ) ≈ 0 can be understood from the fact that the expansion is performed in the limit m 2 Z m 2 t s, |t|, |u|. In particular, for scattering angles θ ≈ 0 or θ ≈ π, the parameters |t| and |u| are not guaranteed to be large compared to m 2 t . The Padé improved expansion substantially cures this problem: the agreement with our exact result is close to perfect for the high energy samples (Bottom Panel) and good within a few percent for the intermediate energy samples as long as the scattering is relatively central (Top Right Panel).
λ 3 λ 4 , as well as interference terms at two-loops, V λ 3 λ 4 , for different outgoing helicities compared against the expansion results (only at two-loops). We find good agreement with the expansions in the relevant regions. We observe that the Padé approximation does not agree with the full result equally well for all helicities. Indeed, for the dominant V (2) 00 helicity configuration the approximation works well rather close to the top-quark threshold. However, for the suppressed V (2) +− and V (2) +0 configurations the approximation begins to visibly deteriorate for √ s/m t 3. It is interesting to observe that the mode with longitudinal polarisation for both the Z bosons dominates both V (1) We also see a rapid increase in V (1) 00 and V (2) 00 past the √ s = 2m t threshold, where the top quarks can be produced on-shell.
In figures 10, 11 and 12 we show our results for the polarized interference terms as a function of cos(θ) and compare them against expansion results for different fixed values of √ s. Note that in many plots the one-loop and two-loop results agree so perfectly when scaled accordingly, that the green points are exactly on top of the black points. In figure  10, we observe that the large top-mass expansion approximates the exact result very well below the 2m t threshold also as a function cos(θ). For the intermediate energy considered in figure 11, the Padé result for small top-quark mass agrees overall rather well with the exact results. As expected from the previous discussion, some deviations are visible for nondominant final state helicities or non-central scattering. Further, we note an asymmetry in the cos(θ) dependence, which is a consequence of the asymmetric expansion used to  construct the Padé approximation. In figure 8, this asymmetry is absent by construction since in this unpolarised case the Padé approximation is calculated for a fixed value of √ s and p T and used for both the forward and backward directions. For the high energy in figure 12, we see excellent agreement between the angular dependence of the Padé result and that of the exact result. As visible from the figure, this is a significant improvement over the angular dependence of the power series approach to the small mass expansion for less central scattering angles.

Conclusions
In this paper, we have presented a calculation of the two-loop top-quark corrections for the process gg → ZZ. Maintaining exact dependence on the top-quark mass, we calculated the helicity amplitudes in terms of finite integrals, which we evaluated using numerical quadrature. For reduction to master integrals, we employed finite field techniques and syzygies which avoid the introduction of squared propagators ("dots"). We presented a new computational method to find these syzygies with linear algebra. We considered finite linear combinations of dimensionally regularized Feynman integrals and presented a novel algorithm to systematically construct them. These linear combinations possess convergent parametric integral representations for = 0 and are formed from building blocks which may involve irreducible numerators, higher powers of propagators, dimensionally shifted integrals, and subsector integrals. The resulting parametric integrand can be expanded in allowing for direct numerical integration. We employed pySecDec and found that such finite linear combinations can significantly improve the numerical performance of the quadrature, similar to what was observed in the case of dimensionally shifted integrals with additional dots. The new approach allows us to stay in a range of integrals which may be considered more natural in the context of the amplitudes themselves. We emphasize that our method is fully automated and works for arbitrary loop order and number of external legs.
Our results for the two-loop amplitudes show good agreement with the large m t expansion and the small m t expansions in the regions, where they are expected to be valid. At moderate energies and for non-central scattering at higher energies, we find that the small m t power series expansion differs substantially from our result. In contrast, the Padé improved results [36] give a very good approximation to our results for a much larger region of phase space.
We observed that the quantitative and even qualitative behavior of the 2-loop finite remainders is rather sensitive to the choice of infrared subtraction terms. In particular, admixtures of 1-loop contributions may actually dominate the overall behavior of the 2-loop remainders and smooth 2-loop threshold effects. As a consequence, the level of agreement between the available approximations and our exact result depends significantly on the choice of the subtraction scheme.
Our amplitudes provide the last major building block required to include the full topquark mass effects in the next-to-leading order cross section for ZZ production in gluon fusion.

A Numerical checks
In this appendix, we present details for the numerical pole checks for our basis of finite integrals in Kreimer's anti-commuting γ 5 scheme. For the UV renormalised two-loop form factors prior to IR subtraction, we observe analytical pole cancellation through to order 1/ 4 and very precise numerical cancellations at order 1/ 3 . For the 1/ 2 and 1/ poles, the following table shows our results compared against predicted IR poles (5.9, 5.10) as well as the 0 term (before IR subtraction) for the Euclidean point s/m 2 t = −191, t/m 2 t = −337, m 2 Z /m 2 t = −853, m t = 1. The digits in parentheses for the 0 term denote the uncertainty in the last digit. +3.055600405 · 10 −4 +1.158849558 · 10 −3 + 1.919890357 · 10 −3 i −4.35036(1) · 10 −3 + 6.0546699(5) · 10 −3 i Pred.
−4.235989677 · 10 −8 −1.659621000 · 10 −7 − 2.661550810 · 10 −7 i In the following table, we compare the 1/ 2 and 1/ poles of the 2-loop form factors against predicted IR poles (5.9, 5.10) as well as provide the 0 term (before IR subtraction) for a point in the physical region with s/m 2 t = 142/17, t/m 2 t = −125/22, m 2 Z /m 2 t = 5/18, m t = 1. We only note here that we observe an improved agreement for the physical combinations of these form factors. The digits in parentheses for the 0 term denote the uncertainty in the last digit.

B Subtraction scheme dependence
In the main text of this article we define the finite remainders for the amplitudes according to the "q T scheme" [117]; see section 5.1. As described in section 6, the choice of scheme can significantly affect the results and the level of agreement with the available approximations.
In this appendix we demonstrate this effect explicitly by presenting a selection of our results using an alternative definition of the finite remainders in Catani's original scheme [116]. At the level of form factors, the finite remainders in Catani's original scheme are obtained from their "q T scheme" analog in (5.11) according to where see also eqs. 4.9 and 4.10 in [26]. The transformation of the helicity amplitudes follows the same pattern. For the interference terms considered in eqs. (6.3) and (6.4) an additional factor of 2 needs to be taken into account and the term due to iπβ 0 does not contribute. As can be seen in figures 13, 14, 15, and 16, the 2-loop corrections can show a rather different qualitative behaviour than the 1-loop corrections in Catani's original scheme. This is in contrast to the corresponding results in the "q T scheme" in figures 9, 10, 11, and 12. Moreover, the relative agreement between the expansion results and our exact calculation depends greatly on the choice of scheme for the finite remainder; it is significantly better in the "q T scheme" than in Catani's original scheme. In order to assess which relative error reflects better the resulting relative error on physical observables, the corresponding real radiation contributions would need to be taken into account in the respective scheme as well.  Figure 16. The cos(θ) dependence of 1-loop and 2-loop interferences for polarised ZZ production in gluon fusion at √ s/m t = 4.703. The small top-quark mass expansion (to order m 32 t ) and Padé improved expansion [36] are shown for comparison. Here we reproduce the top left and bottom right panels of figure 12 using Catani's original subtraction scheme [116].