Two loop correction to interference in $gg \to ZZ$

We present results for the production of a pair of on-shell Z bosons via gluon fusion. This process occurs both through the production and decay of the Higgs boson, and through continuum production where the Z boson couples to a loop of massless quarks or to a massive quark. We calculate the interference of the two processes and its contribution to the cross section up to and including order O(alpha_s^3). The two-loop contributions to the amplitude are all known analytically, except for the continuum production through loops of top quarks of mass m. The latter contribution is important for the invariant mass of the two Z bosons, (as measured by the mass of their leptonic decay products, m_4l), in a regime where m_4l>2m because of the contributions of longitudinal bosons. We examine all the contributions to the virtual amplitude involving top quarks, as expansions about the heavy top quark limit. Comparison with the analytic results, where known, allows us to assess the validity of the heavy quark expansion, and it extensions. We give results for the NLO corrections to this interference, including both real and virtual radiation.


Introduction
The production of four charged leptons is a process of great importance at the LHC. It was one of the discovery channels of the Higgs boson at the LHC. It also provides fundamental tests of the gauge structure of the electroweak theory through the high-energy behaviour. Four charged leptons are predominantly produced by quark anti-quark annihilation; the mediation is by photons or Z bosons dependent on the mass of the four leptons, m 4l .
A smaller contribution, which however grows with energy is provided by gluon-gluon fusion. The Higgs boson is of course produced in this channel; in the Standard Model (SM) this occurs predominantly through the mediation of a loop of top quarks. As pointed -1 -JHEP08(2016)011 Representative diagrams for the ZZ production. In the following we will suppress the Z-decays to leptons. out by Kauer and Passarino [1], despite the narrow width of the Higgs boson, the Higgsmediated diagram gives a significant contribution for m 4l > m H . If we examine the tail of the Higgs-mediated diagrams there are three phenomena occurring: • The opening of the threshold for the production of real on-shell Z bosons, m 4l > 2m Z .
• The region m 4l = 2m, (m is the top quark mass) where the loop diagrams develop an imaginary part.
• The large m 4l region, m 4l > 2m, where the destructive interference between the Higgs-mediated diagrams leading to Z bosons and the continuum production of onshell Z bosons is most important.
A feature of this tail is that it depends on the couplings of the Higgs boson to the initial and final state particles but not on the width of the Higgs boson. Assuming the couplings of the on-and off-peak Higgs-mediated amplitudes are the same, it has been proposed to use this property to derive upper bounds on the width of the Higgs boson [2]. Note that models with different on-and off-peak couplings can be constructed [3].
In the following we shall refer to the production of the bosons V 1 , V 2 . Gluon-gluon fusion first contributes to the cross section for electroweak gauge boson production pp → V 1 V 2 as shown in figure 1(c)-(e) at O(α 2 S ), which is the next-to-next-to-leading-order (NNLO) with respect to the leading-order (LO) QCD process shown in figure 1(a); no two-loop gg → V 1 V 2 amplitudes participate in this order in perturbation theory.
In the context of the Higgs boson width, however, the interference between the Higgsmediated Z boson pair-production and the Standard Model continuum at next-to-leadingorder (NLO) QCD already requires knowledge of the one-and two-loop gg → (H →)V 1 V 2 amplitudes. The requirement for more precise estimates to the Higgs boson width were -2 -JHEP08(2016)011 emphasised in [4][5][6]. Signal-background interference effects beyond the leading order have been considered in ref. [7] for the process gg → H → W + W − for the case of a heavy Higgs boson.
In this work we will limit ourselves to the Z boson pair final state, due to its importance at the LHC. At LO [8] and NLO [9][10][11][12] the amplitudes for single Higgs boson production have been known for quite some time. At LO, the amplitude for the SM continuum gg → ZZ process occurs via massless and massive fermion loops and results are available in each case [13][14][15][16].
The situation, however, is different for the NLO continuum process, although vast progress in terms of two-loop amplitudes has been made [17][18][19][20][21][22]. Recently two-loop gg → ZZ amplitudes 1 via massless quarks became available [21,22]. The complete computation of two-loop amplitudes with massive internal quark loops, on the other hand, is commonly assumed to be just beyond present technical capabilities. Although the contribution of the top quark loops to these diagrams is smaller than the contribution of the light quarks in the region just above the Z-pair threshold, in the high m 4l region the amplitude is dominated by the contributions of longitudinal Z bosons that couple to the top quark loops. Recently a first heavy top quark approximation for the two-loop gg → ZZ amplitude with internal top quarks was published [6]. In that work only the leading term in the s/m 2 expansion was considered. In that approximation, the vector-coupling of the Z boson to the top quark does not contribute. In addition an approximate treatment of this process at higher orders, based on soft gluon resummation, was presented in ref. [23].
In the present work we will push this analysis further. We start by presenting our results for the LO and NLO Higgs-mediated ZZ production in terms of the s/m 2 expansion in section 2, despite the fact that the full result is known. This part is required for the later interference with the SM continuum. Furthermore, it is well suited to introduce our notation in section 2.1 and to assess the validity of the approximation methods with respect to the exact known (N)LO amplitudes in section 2.2.
The results for the LO and virtual NLO contributions to the SM continuum with massive quark loops will be given in section 3 as a large-mass expansion (LME) with terms up to (s/m 2 ) 6 . We will limit our discussion to the interference between the Higgsmediated term and the continuum term. Similar to [6] we will consider on-shell Z bosons in the final state. A theoretical predictions for off-shell Z bosons would be optimal, but in order to reduce the number of scales in the problem, we restrict ourselves to on-shell Z bosons. Since we are primarily interested in the high-mass behaviour this is an appropriate approximation. A limited number of scales is beneficial when we consider the extension of our approach to a full calculation. In section 4 we summarize our treatment of the real radiation contribution, which makes use of results already presented in ref. [16].
The results of our calculation, including loops of both massless and massive quarks, will be presented in section 5. We will compare the effects of the NLO corrections to the interference contribution with the corresponding corrections to the Higgs diagrams alone.

JHEP08(2016)011
In addition, we will discuss the impact of our results on analyses of the off-shell region that aim to bound the Higgs boson width.
All expansion results from section 2 and section 3.4.1 are provided via ancillary files on arXiv as FORM and Mathematica readable code.
2 Higgs production in gluon-gluon fusion and decay to ZZ In this section we give a detailed discussion of single Higgs boson production at LO and NLO QCD and its subsequent decay to a pair of on-shell Z bosons. As mentioned earlier the LO and NLO amplitudes for single Higgs boson production have been known for a long time; either approximate results in terms of Taylor expansions in the inverse of the top quark mass s/m 2 [8,12,[24][25][26][27][28] or results keeping the exact top mass dependence [12,29].
It is understood that, whenever feasible and available, the exact results for LO and NLO amplitudes are used. However, we are mainly interested in approximations to the interference contributions Re A LO B (N)LO , where A denotes the Higgs-mediated and B the SM continuum amplitude. Since no exact results are available for B NLO we will use the, so-called, large-mass expansion [30] as an approximation of the SM continuum. Hence, for consistency, we also perform the expansion of the Higgs-mediated amplitude A to high powers in s/m 2 . Expansion of the two-loop Higgs-mediated amplitude A NLO and its comparison to available results from the literature provides moreover a helpful check of our expansion routines due to the general structure of the LME.
Furthermore, the large-mass expansion in powers of s/m 2 is formally only valid below the threshold of top quark pair-production, as m is assumed to be much larger than any other scale in the problem, e.g. s m 2 . As extensively discussed in literature the naive LME can be drastically improved at (and even far above) threshold by taking the next mass threshold into account, see ref. [30] and references within, or by rescaling the approximated NLO result by the exact LO result, see e.g. refs. [31,32]. We will address this issue in section 2.2.3 and try to draw conclusions for the SM continuum.

Preliminaries
The amplitudes for single Higgs boson production are illustrated in figure 2 for the one-loop and two-loop case. The largest contribution is due to the internal massive top quark loop; in the following we will ignore the contribution of other quarks for the Higgs production process. The gg → H amplitude, with color (Lorentz) indices A, B(α, β) for the initial state gluons, can be written as such that the reduced matrix element A 0 (α 0 S , m 0 , µ, ) is dimensionless and can be expressed as a function of µ 2 /s and r t = m 2 /s. The bare on-shell amplitudes admit the perturbative expansion where we introduced the parameter from dimensional regularisation in d = 4 − 2 spacetime dimensions and µ to keep the amplitudes dimensionless. The calculation is performed in Conventional Dimensional Regularisation (CDR) and the following definition of the d-dimensional loop integral measure is used in accordance with the M S-scheme, to avoid the proliferation of unnecessary γ E − log(4π) terms. The ultraviolet (UV) renormalised amplitudes are given by where Z g denotes the on-shell gluon renormalisation constant. The Htt vertex is renormalised, according to [33], by g 0 H = Z m g H with g H being the Yukawa coupling for the top quark. The bare top quark mass is related to the renormalised mass, m, by m 0 = Z m m. The necessary on-shell renormalisation constants are given by -5 -

JHEP08(2016)011
with T F = 1/2. See appendix A of [34] and references therein for more information. The mass renormalisation enters as an overall factor in eq. (2.5) because of the renormalisation of the Yukawa coupling, and also implicitly in the relationship between the bare and renormalised mass. We will always present mass-renormalised results in the following. The strong coupling constant is renormalised in the M S-scheme according to with [34] Z (n f ) where n f = 6 denotes the number of fermions and β (n f ) 0 the coefficient of the beta function.
The explicit scale dependence of the renormalised strong coupling constant α (n f ) S (µ) is dropped in the following to simplify our notation. All of our quantities are computed in five-flavour (n l = 5) QCD. Hence, we decouple the top quark from the QCD running via with n l the number of light quarks. After UV renormalisation the two-loop amplitude still contains divergences of infrared origin. The structure of these divergences is, however, completely understood at two-loop level. The finite remainder is defined by infrared (IR) renormalisation The infrared renormalisation matrixẐ gg is taken from [34][35][36] and reads for the gluongluon initial state with colourless final state in terms of the renormalised strong coupling constant (2.14) In the end we are interested in the amplitude for the process (2.17) where we have defined an overall normalisation factor, (2.18) From this it is straightforward to square the amplitude to obtain the result for the Higgsmediated diagrams alone. The sum over the polarisations of the gluons and the Z bosons of momentum p can be performed as usual with the projection operators, Using these projectors we get the subsidiary result (2.20) Including also the sum over colors yields the matrix element squared for the signal in this channel, (The statistical factor for identical Z bosons is not included).
where we use the notation r Z = m 2 Z /s and N A = N 2 c − 1 = 8.

Large-mass expansion and improvements
Using the aforementioned conventions we can compute the leading-and next-to-leadingorder amplitude A (1,2) (m, µ, ) for single Higgs boson production. Although we always work with the loop measure S = exp( γ E )(4π) − we factor out in the results presented below to keep factors of π 2 implicit. The dimensional dependent factor c Γ denotes the somewhat more natural loop measure, because it cancels exactly the Γ(1 + ) factor obtained by the loop integration.
The exactly known leading-order result in d-dimensions (d = 4 − 2 ) yields [8,11,27,37,38] where s = (p 1 + p 2 ) 2 . The definitions of the integrals B 0 and C 0 are given in appendix A. The essential idea of the large-mass expansion based on the method of expansion by regions [30] is that the integration domain is divided into different regions where the loop momenta are soft, k i ∼ p i m or hard, p i k i ∼ m. The external momenta p i m are always assumed to be small. In the expansion of one-loop integrals only the region of a hard loop momentum k 1 ∼ m exists, because all propagators are associated with the large mass m. As a result the one-loop expansion consists only of a naive Taylor expansion and its result is given in terms of simple massive one-loop vacuum integrals.
The two-loop integral expansion is more involved since the hard as well as the soft region must be considered. The first region results, with the help of [39], in scalar massive two-loop vacuum integrals. The soft region produces a product of massive one-loop vacuum integrals and massless one-loop bubble and triangle integrals. All occurring integrals are well known and, although, the intermediate expressions become huge, the final results are remarkably simple, as can be seen below. We use our own fully automatic in-house software to perform the large-mass expansion, relying extensively on the features of FORM [40] and Mathematica. For a similar approach to Higgs boson pair-production, see e.g. [41].
Using the large-mass expansion for the B 0 and C 0 integral, given in section A, the corresponding expansion of the full result for A (1) (m, µ, ) in d dimensions is Similarly the two-loop result can be expressed in terms of the leading-order amplitude Ā (1) (m, µ, ) = (S c Γ (µ 2 /m 2 ) ) −1 A (1) (m, µ, ) and with only mass renormalisation in- The first terms of eq. (2.24) and eq. (2.25) fully agree with available results in the literature [27,28]. Especially the NLO corrections presented in [27] cover terms in the expansion up to O 1/r 2 t , 2 and we find full agreement with our results for the amplitudes as well as the cross sections. The analytic results for the exact LO and NLO amplitude A, keeping the full top mass dependence, can be taken from [11,42]. 2 The NLO results for the virtual amplitude have also been checked by our own independent program, using GiNaC [43] to evaluate the harmonic polylogarithms. This serves as a further independent check of the mass expansion results in eq. (2.24) and eq. (2.25). This agreement will be illustrated in section 2.2.3.
The radius of convergence of the large-mass expansion is given by s/(4m 2 ) 1. The polynomial growth leads to an extremely good convergence below and close to threshold of top quark pair-production, as shown later.

Rescaling with exact leading-order result
Above threshold, however, naively no convergence with respect to the exact result can be expected. At least two procedures exist which lead to major improvements in terms of convergence of the expanded result even above threshold. 3 We recall these procedures in this subsection and the next.
A well known method of extending the naive large-mass expansion of the NLO cross section beyond its range of validity relies on factoring out the LO cross section with exact top mass dependence, (2.26) The numerator and denominator are expanded to the same order in 1/r t . It was argued for single Higgs boson production in [31] and for Higgs boson pair-production in [32] that varying N in the above formula allows to check for additional power corrections. Including sufficient orders in the expansion should lead to stable approximations σ NLO imp,N . 2 The overall sign of the NLO term differs between the published paper [11] and the thesis of Beerli [42].
We believe that the sign in the latter is correct, which is also supported by the comparison with the NLO results using the large-mass expansion [27,28]. 3 The region above threshold could also be approximated by fitting a suitable ansatz to the high-energy limit [31,44,45]. This, however, would require additional knowledge of the high-energy behaviour and is beyond the scope of this work.

JHEP08(2016)011
The method relies on the expansion of numerator and denominator in eq. (2.26) and evidently, requires the knowledge of all of the ingredients in terms of series expansions. Although this requirement usually does not pose any problem per se it might turn out to be disadvantageous in certain cases. In our particular case at hand, we require the SM continuum as well as the Higgs-mediated amplitude as large-mass expansions. Certainly the Higgs-mediated amplitude is well known at LO and NLO including its full top mass dependence. Any approximation of this amplitude poses a potential threat of introducing unnecessary uncertainties. We will discuss this point further in section 3.5 and see that the method introduced in the next section provides a way to circumvent this issue.

Conformal mapping and Padé approximants
Having sufficiently many terms in the 1/m expansion at hand allows for a more powerful resummation method, the Padé approximation [30,[46][47][48][49]. The univariate Padé approximant [n/m] to a given Maclaurin series with a non-zero radius of convergence z 0 is defined via the rational function such that its Taylor expansion reproduces the first n+m coefficients of f (z); the coefficients b i and c i are uniquely defined by this expansion. The advantage of Padé approximants over other techniques, e.g. Chebyshev approximation, lies in the fact that they can provide genuinely new information about the underlying function f (z), see [49] for more information.
The downside of Padé approximants is their uncontrollability. In general, there is no way to tell how accurate the approximation is, nor how far the range z can be extended. Computing the Padé approximants [n/n] or [n/n ± 1] for different orders n allows, at least, checking the stability of the approximation. We will refer to [n/n] as diagonal and to [n/n ± 1] as non-diagonal Padé approximants in the following.
Although the Padé approximation can be directly applied to eq. (2.27), it is advantageous to apply a conformal mapping [46] first. The amplitudes at hand, gg(→ H) → ZZ, with z = s/m 2 develop a branch cut starting from z 0 = 4 and extending to +∞ due to the top quark pair-production threshold. Applying the mapping, eq. (2.29), transforms the entire complex plane into the unit circle of the w-plane, such that the upper (lower) side of the cut corresponds to the upper (lower) semicircle and the origin of the original z-plane is left unchanged. The initial power series can now be transformed into a new series in w [30] f and, subsequently, its Padé approximants computed. We will illustrate these features using the example of single Higgs boson production in the next section.

Comparison of LME with full result
Let us briefly compare the results from the large-mass expansions, eq. (2.24) and eq. (2.25), and their, previously discussed, improvements to the known LO and NLO QCD result with full top mass dependence [9][10][11][12]. We include the subsequent H → ZZ decay, as given in eq. (2.16), perform the UV+IR renormalisation and compute the phase space integral over eq. (2.21) including all corresponding phase space factors and coupling constants. The NLO contribution so obtained is not physical, since we neglect the real-radiation contribution for now. Considering the obtained finite parts of the LO and virtual NLO corrections alone, on the other hand, allow to better verify the validity of our approximations. To be specific, we set (2.32) We utilise the CT 10nlo PDF set [50] within LHAPDF [51] to determine α S (µ f ) and use the input parameters where S and s denote the hadronic and partonic center-of-mass energy, respectively. The orange curves in figure 3 depict the large-mass expansion results of eq. (2.32) for the LO and the NLO case, where each 4 finite remainder F (1,2) A is expanded up to 1/m 20 . A minimum cut √ s ≥ 2m Z has been imposed and the threshold for top quark pair-production is given by s/m 2 = 4. The relative deviation of the approximated results with respect to the exact result are shown in the bottom plots. The large-mass expansion describes the exact LO and virtual NLO results up to the top threshold very well, with only 5% deviation at LO and 7% at NLO at s = 4m 2 . 4 The ambiguity between expanding the product F consists only of power corrections which are numerically negligible. We checked that the difference in ∆σ/σ at threshold of both approaches is 1%. The same arguments hold for the series expansions including the conformal mapping.  The second choice of improving the naive LME is given by the rescaling from eq. (2.26). The results are shown in figure 4. The exact virtual NLO result is again shown in black. The rescaled LMEs are indicated by the shaded grey area and its envelope is given by the expansions σ NLO imp,1 (orange) and σ NLO imp,10 (blue). Although the heavy-quark approximation σ NLO imp,1 gives a reasonable estimate of the exact result above threshold it fails to describe the threshold behaviour and peak structure of the exact result. At threshold the deviation is 10%. Taking higher orders in the expansion into account improves the threshold prescription, with 3% deviation for σ NLO imp,10 at threshold, but worsens the trend for higher energies. In both cases we find more than 20% deviation for s/m 2 > 20.
We end this section by drawing our conclusions from the results presented. We see that, at least in the single-scale Higgs boson production and having a sufficient number of terms in the LME at hand, applying the conformal mapping (and the Padé approximation) yields excellent prescriptions of the exact results. The conformal mapping is imperative, whereas the additional Padé approximants give only small improvements in terms of uncertainty reduction and stability of the approximations. We conclude that we should favour these approximations over the rescaling method.
-12 - One important point to notice, however, is that the kinematics change when moving from the single Higgs boson production to the SM Z boson pair-production. 5 Therefore, the results discussed here may not necessarily transfer easily. Still, the comparisons within this section should give an idea of the validity of the improved large-mass expansions. We will discuss analogous considerations for the Z boson pair-production in section 3.5. available only recently [17][18][19][20][21][22]. Due to the complexity of the computation and present technical limitations no full two-loop correction to the amplitudes with massive internal quarks is presently known. The authors of [6] made the first attempt in approximating the virtual NLO corrections with internal top quarks. Their results, however, includes only the first term of the 1/m expansion. At this order contributions from the vector coupling of the Z bosons to the quarks are neglected completely. This is not necessarily troubling since the vector coupling contribution is a f /v f ∼ 2.5 times smaller than the axial coupling contribution. However, to fully incorporate the physics of the Z boson interactions and to give an estimate of power corrections s/m 2 we compute the virtual two-loop corrections up to O r −7 t . We keep the Z bosons on-shell, sum over their polarisations and project onto the tensor structure of the gg → H → ZZ amplitude (eq. (2.17)) since we are only interested in the interference of both.
This section is structured as follows: in section 3.1 we give our definitions of the SM ZZ amplitude, as far as the conventions differ from section 2.1. The leading-order and next-to-leading-order results are presented in section 3.3 and section 3.4, respectively. The latter is divided into two parts; the first consists of diagrams where both Z bosons couple to one fermion line and the second handles anomaly style diagrams where a single Z boson is connected to one fermion string.

Preliminaries
The on-shell Z boson pair-production in gluon-gluon fusion via the heavy top quark loop can be completely expressed in terms of kinematical invariants or equivalently, using the on-shellness condition, by the rescaled variables The SM continuum amplitudes B 0,AB µναβ (α 0 S , m 0 , µ, ) admit the same perturbative expansion as given in eq. (2.3) for the Higgs-mediated process. The bare amplitudes are renormalized in accordance with eqs. (2.5)-(2.14), omitting the superfluous Higgs vertex renormalisation. As mentioned earlier we project onto the tensor and color structure of the Higgs-mediated amplitude (eq. (2.17)) with where N A = N 2 c − 1 = 8 and P Z,αβ (p) from eq. (2.19). We shall consider a single quark of flavor f to be circulating in the quark loop. The Standard Model coupling of this fermion to a Z boson is given by, The superposition of vector and axial coupling allows to write the scattering amplitude as where we factored out the normalisation factor from eq. (2.18). The mixed coupling structure v f a f vanishes due to charge parity conservation. With the amplitudes outlined above it is straightforward to compute the interference.
Writing eq. (3.7) in this way establishes that A(α

Projected exact result at one loop
The leading-order amplitude for the SM continuum production of two Z bosons is known exactly in d = 4 − 2 dimensions. The usual normalisation factor eq. (2.22) is chosen. We split the result, according to eq. (3.6), into vector-vector (V V ) and axial-axial (AA) contribution. Table 1. Scalar integrals occurring in full LO SM continuum ZZ production. The momenta are defined as q 1 = p 1 , q 2 = p 2 , q 3 = −p 3 and q ij = q i + q j . The scalar integrals are defined in appendix A.
The notation for the scalar integrals B, C and D is given in table 1. We re-introduced factors of s in eq. (3.8) and eq. (3.9) to indicate the correct dimensionality of the expressions. We note that, in contrast to the case where the Z bosons are off-shell and their decays included, these formulae for the interference take a very simple form. Eq. (3.8), (3.9) extend the results of ref. [16] to include the terms of order 1 and 2 .

Large-mass expansion at one loop
Equivalently, eq. (3.8) and eq. (3.9) can be expressed by means of the large-mass expansion.
The result for the vector-vector part yields The result for the axial-axial part is The leading term in the vector-vector expansion is sub-dominant with respect to the axialaxial part. The reason for this difference has been given in [6].

Large-mass expansion at two loops
The two-loop SM continuum amplitude consists in total of 93 + 16 non-zero diagrams. 93 diagrams belong to topologies where both Z bosons couple to the same fermion string, as illustrated in figure 5. Due to momentum conservation and assuming an anti-commuting γ 5 in d-dimensions, no γ 5 contribution arises in the fermion traces of the respective diagrams. The large-mass expansion results for the vector-vector and axial-axial part of these diagrams are shown in section 3.4.1.
The remaining 16 anomaly style diagrams belong to the topology shown in figure 6, where the Z bosons couple to distinct fermion lines. These diagrams must, in principle, be handled with care when using dimensional regularisation due to the non-conservation of the axial-current. Furthermore, contributions from each quark-doublet have to be considered simultaneously. Only the sum over one quark-doublet leads to a gauge anomaly free theory.

JHEP08(2016)011
In case of massless quark doublets all contributions vanish and we only have to consider the third-generation quark doublet, i.e. top and bottom quarks. Results for these diagrams are presented in section 3.4.2.

Non-anomalous diagrams
In this section we give explicit formulae for the large-mass expansions for the sum of the 93 anomaly free diagrams. Including again only mass renormalisation, setting XX (r t , µ, ) and log(−r t ) = log m 2 /(−s − i ) , we can write the divergent two-loop V V part as  And the AA part as The leading term in eq. (3.13) can be compared to the projected results of [6]. We find agreement with their formula. 6 We also performed a consistency check of the renormalisation scale dependence of the presented two-loop expansions by means of the technique given in appendix B.

Anomalous diagrams
The two-loop gg → ZZ amplitude contains, in addition, two topologies which consist of products of one-loop sub-diagrams. On the one hand diagrams containing gluon selfenergy contributions vanish due to color conservation. The diagrams in figure 6, on the other hand, give a finite mass dependent contribution as long as both Z bosons couple to distinct fermion loops. These diagrams are proportional only to the axial coupling of the Z bosons to fermions; the vector component vanishes due to C invariance (Furry's theorem). The diagrams have been omitted in the previous section since they can be computed with their full top mass dependence and, therefore, need no large-mass expansion [52,53].
In brevity we repeat the results from [53] and give the result in terms of our conventions. Let us denote the amplitude for a Z coupling to two gluons by T µνρ AB . We calculate the triangle shown in figure 7, where all momenta are outgoing q 1 + q 2 + q 3 = 0 and to begin with q 2 1 = 0, q 2 2 = 0. The result for the two triangle diagrams (including the minus sign for 6 Both eq. (5) and eq. (7) of [6] contain typographical errors.
where τ f = ±1/2 and, The most general form of Γ consistent with QCD gauge invariance, can be written as, (3.17) By direct calculation it is found that F 4 = 0.
Contracting with the momentum of the Z boson we find that, 18) The divergence of the axial current is found by direct calculation to be, (q 3 ) ρ Γ µνρ = 4m 2 C 0 (q 1 , q 2 ; m, m, m) + 2 Tr[γ µ γ ν q 1 q 2 γ 5 ] (3.19) showing the contribution of the pseudoscalar current proportional to m 2 and the anomalous piece. Summation over one complete quark doublet (τ f = ±1/2) cancels the anomaly term and solely the piece proportional to the top mass remains. For the particular case at hand we are interested in on-shell Z's and in q 2 2 = ε 2 · q 2 = 0, ε 3 · q 3 = 0, q 3 = −q 1 − q 2 , so we get a contribution only from F 1 . The result for F 1 is F 1 (q 1 , q 2 , m) = 1 2q 1 · q 2 2 + 4m 2 C 0 (q 1 , q 2 ; m, m, m)

JHEP08(2016)011
We further define a subtracted F 1 to take into account the contribution of the top and the bottom quarks, Analogous to eq. (3.4) we define the projected matrix element for the anomaly piece (3.23) The amplitude defined in eq. (3.23) is UV and IR finite and requires no renormalisation. Including the effect of both the b quark (taken to be massless) and the t quark we obtain (No statistical factor for identical Z bosons is included).
where N is given in eq. (2.18). Again we include the factors s 2 to indicate the correct dimensionality of F 1 (p 1 , p 2 , m). For completeness we also give the mass expansion of eq. (3.24) in case only the top quark contribution is considered, i.e. F 1 (p 1 , p 2 , 0) → 0. As expected the expansion starts at 1/m 4 .

Visualisation of large-mass expansion results for gg → ZZ
Let us turn towards the graphical representations of the large-mass expansion results for the SM continuum, eqs. (3.10)-(3.13), and their improvements. We proceed analogously to section 2.2.3 and compute the UV+IR renormalised version of eq. (3.7) and again integrate over the ZZ phase space. The setup from eq. (2.33) is utilised. Since we focus our discussion in this section mainly on the different improvements of the large-mass expansions we, again, do not take into account the full NLO correction. We merely focus on the unknown virtual massive two-loop contribution of the SM continuum interfered with the Higgs-mediated process. That is, we set B (m, µ) , (3.26) which also excludes the anomaly style contribution from eq (3.24) since this part can be computed without the necessity of any approximation.
It is important to notice the following conventions for our approximations using Padé approximants below. As in section 2.2.3 the Padé approximants are computed at amplitude level for each finite remainder F A,B , including the conformal mapping. 7 We know from our previous discussion that the best approximation of the LO as well as the virtual NLO contribution of the Higgs-mediated process is given by F (1,2) A, [5/5] . It is understood that we will always use this approximant in the following considerations. In principle, we can also substitute the approximated Higgs-mediated amplitude F (1) A, [5/5] with its exact LO result. Doing so would remove any uncertainties from the Higgs-mediated contribution. On the other hand the numerical difference between both approaches is negligible as discussed in section 2.2.3.
The vector-vector part of the SM continuum gives only a minor contribution to the total cross section, σ V V /σ AA ∼ 10 −3 . This relies on the fact that the mass expansion of the V V part starts only at 1/m 4 whereas the AA part starts at 1/m 2 and additionally The interference including the exact top mass dependence is only known at leadingorder, which is shown in the left panel of figure 8. Comparing the exact result (black) and its naive large-mass approximation up to 1/m 12 (orange) shows excellent agreement up to s ∼ 3m 2 , with approximately 1% deviation from the exact result. At threshold the deviation rises to 12%. In contrast the Padé approximant F     Ineptly this is the region of interest for our later analysis of the Higgs boson width. Going to large values of s the deviations inevitably become larger, but the contribution to the cross section is small due to the suppression by the flux.
This situation seems to continue in case of the next-to-leading-order large-mass expansion as shown in the right panel of figure 8. Evidently no exact result is available and we have to rely on the approximate results. All Padé approximants [2/2], [2/3], [3/2], [3/3] for F (2) B show a stable trend over the entire s/m 2 range. The deviations between the diagonal and non-diagonal Padé approximants are again indicated by the shaded grey area and the approximant F (2) B, [3/3] is shown in orange. The steeper rise near the top threshold suggests a better description of the actual threshold properties of the NLO result with exact top mass dependence in contrast to the naive large-mass expansion (black). Comparing the trend above threshold with its analogous LO situation we can only guess that we have to expect comparable deviations from our Padé approximations with respect to the unknown exact NLO result.
We can also consider rescaling the NLO large-mass expansion as described in eq. (2.26). The resulting curves are shown in the left panel of figure 9. To guide the eye we also include F (2) B, [3/3] (black). The envelope of the different orders n in the expansion σ NLO imp,n is shown as grey area. For s ≤ 20m 2 the envelope is determined from n = {1, . . . , 6}, whereas for s > 20m 2 we only use n = {1, . . . , 5} due to the instabilities for n = 6 in the high energy regime. The most interesting curves, namely the heavy-quark approximation n = 1 and the highest order in the expansion n = 6, are shown in orange and blue, respectively. Factoring out the exact LO result seems to give a more natural description of the threshold behaviour and peak structure in comparison to the plain use of the Padé approximation.  The origin of the numerical instabilities of the n = 6 expansion is probably due to delicate numerical cancellations in the (s/m 2 ) 6 coefficients. One could try to cure this problem by switching to a higher numerical precision or by a proper economisation [49] of the power series. With the Padé approximation we already have an excellent method at hand and we adopt the idea of factoring out the exact LO interference,  figure 9, right panel. We immediately see the advantages of this approach. Firstly we also get a similar, more natural behaviour at threshold and of the peak structure above threshold. Secondly we get a stable result across the entire range of s/m 2 . The grey area is again given by the envelopes due to the variation between the (non-)diagonal Padé approximants [2/2], [2/3]

, [3/2] and [3/3](orange).
Ultimately by using the Padé approximants in contrast to eq. (2.26) we could entirely remove the uncertainty of having to use an approximation for the involved Higgs-mediated amplitude and fall back to using the exactly known result for F (1) A . Some concluding remarks. In contrast to the purely Higgs-mediated case, section 2.2.3, it turns out that we require the Padé approximation in the interference case. Using the conformal mapping alone without an additional Padé approximant on top gives no reasonable approximation for the quantities discussed above. On the other hand we have seen that we hugely benefit by using Padé approximations due to their stability and the possibility of removing any uncertainty besides the approximated virtual massive two-loop gg → ZZ amplitude.

Real corrections to SM ZZ production
Representative diagrams for the real radiation contributions to this process are shown in figures 10 and 11. The Higgs-mediated diagrams have previously been computed in [56].   They can easily be adapted to our calculation by combining those results with the decay amplitude given in eq. (2.16) and N from eq. (2.18). This procedure, together with the strategy for handling the amplitudes for diagrams without a Higgs boson, is described in detail in [16]. We adopt this implementation here. Our calculation of the pure-Higgs contribution involves the computation of the square of the diagrams shown in figure 10, together with all crossings of the quarks in figure 10 (right) into the initial state. Similarly, the interference contribution includes all crossings of the diagrams shown in figure 11. In principle another contribution to the interference occurs at this order, between tree-level amplitudes for the process qg → ZZq and the qg-initiated diagrams shown in figure 10 (right) and 11 (bottom-left). However this contribution is subleading [16], particularly for high invariant masses of the ZZ system, so we do not consider it here.
The real radiation diagrams contain infrared singularities, of soft and collinear origin, that must be isolated and combined with the corresponding poles in the two-loop amplitudes. This is handled using the dipole subtraction procedure [57].

Results
The individual components of the calculation that have been extensively discussed above have been included in the parton-level Monte Carlo code MCFM [58][59][60]. The bulk of the calculation is performed in a straightforward manner using the normal operation of MCFM at NLO. The exception is the finite contribution to the two-loop amplitude containing a closed loop of massless quarks. Since these contributions are computationally expensive to evaluate, we choose to include their effects by reweighting an unweighted sample of LO events.
For the two-loop amplitudes containing massive loops of quarks the approximations used are as follows. The Higgs amplitude is evaluated using the [5/5] Padé approximant to the LME after conformal mapping. As demonstrated in section 2, this is virtually identical to the exact result. The massive quark box contributions are computed by factoring out the exact LO amplitude according to eq. (3.27), with the Padé approximant corresponding to n = m = 3 in the definition given in eq. (2.28). The anomalous diagrams of section 3.4.2 are not included in the discussion of the massive quark loops below, but instead are accounted for only when the sum of all loops is considered.
For massless quarks circulating in the loop the calculation is simplified by the fact that the entire amplitude is proportional to the combination of couplings (v 2 f + a 2 f ), i.e. in the decomposition given in eq. (3.6) the quantities B 0 V V and B 0 AA are equal. The calculation requires the one-loop master integrals up to 2 , for which all orders results are given in ref. [61] for bubble integrals and refs. [62][63][64][65][66] for the easy box (two opposite off-shell legs). The necessary results for the three-mass triangle with massless propagators and the hard box (two adjacent off-shell legs) can be taken from refs. [67] and [68] respectively. We use the coproduct formalism [69,70] to analytical continue the results to the physical phase space regions. All master integrals have been numerically cross-checked with SecDec [71]. The two-loop master integrals for gg → ZZ are taken from ref. [17] and GiNaC is used to evaluate the polylogarithms. Our results for this contribution agree with the earlier calculation of ref. [22].
The parameters for the following results have already been specified in section 2.2.3. Here we make only one change: our central scale corresponds to the choice µ r = µ f = M ZZ /2, where M ZZ is the invariant mass of the ZZ pair. As an estimate of the theoretical uncertainty we consider variations by a factor of two about this value. We also introduce an uncertainty that is based on our combination of LME and Padé approximants in the calculation of the massive quark loops, that has already been explored in figure 9 (right). In order to obtain a more conservative error estimate we multiply the deviations of the extremal values in the grey area with respect to σ NLO imp, [3/3] by a factor of two. The impact of this variation on the complete NLO prediction for the massive loop is shown in figure 12. Even for this choice, the impact of the approximation is estimated to be less than 20% throughout the distribution. For the remaining plots in this section we no longer show the impact of this uncertainty, but it will be explicitly included in tables 2 and 3 later on.
Results for both the massless and massive quark contributions to the interference, including the effects of scale variation, are shown in figure 13. The interference is negative -28 -JHEP08(2016)011  for both the massless and massive quark contributions and is shown in figure 13 reversed in sign. In both cases the K-factor decreases as the invariant mass of the Z-boson pair increases. The K-factor at small invariant masses is larger for the massless loops; as the invariant mass increases, the NLO corrections are more important for the massive loop. The NLO corrections are larger for the top quark loops and exhibit a stronger dependence on M ZZ . In both cases the NLO result lies outside the estimated LO uncertainty bands and the scale uncertainty is not significantly reduced at NLO.
The relative importance of the massive and massless loops can be better-assessed from the NLO predictions shown in figure 14. At smaller invariant masses, below the top-pair  threshold, the massless loops are most important. Around the top-pair threshold the two are of a similar size, but at high energies the massless loops are insignificant. In contrast, the top quark loop quickly becomes the dominant contribution beyond this threshold and exhibits a long tail out to invariant masses of around one TeV. The full prediction for the interference that is obtained by summing over both massless and top quark loops, as well as the numerically-small anomalous contribution discussed in section 3.4.2, is shown in figure 15. The relative size of the massless and top quark loops discussed above means that the behaviour of the K-factor for the sum of both contributions interpolates between the massless-loop K-factor for small M ZZ and the massive loop one for high M ZZ . It therefore decreases from around 3 at the peak of the distribution to approximately 1.8 in the tail. This is to be contrasted with the K-factor distribution for the pure Higgs amplitudes alone, shown in the right panel of figure 15. In that case the K-factor decreases slowly from around 2.2 at small invariant masses to around 1.8 in the far tail. We note that the K-factor for the Higgs amplitudes alone, and the one for the interference with the top quark loops, is almost identical. In the high-energy limit this is guaranteed to be the case, due to the cancellation between these two processes. This behaviour is shown explicitly in figure 16.
The integrated cross-sections for the interference contributions and the Higgs amplitude squared are shown in table 2. Note that, in this table, the total interference differs from the sum of the massive and massless loops by a small amount that is due to the anomalous contribution. At this level the differences between the effects of the NLO corrections on the various contributions is quite small, with all corresponding to a NLO enhancement by close to a factor of two. The K-factor for the massless loops is slightly larger, which is also reflected in the result for the total interference. In addition to the scale uncertainty, we have also indicated our estimate of the residual uncertainty related to the LME expansion that is indicated in figure 12. The impact of this uncertainty is relatively small, at the level of around 5%, due to the fact that the integrated cross-section is dominated by the region M 2 ZZ 5m 2 where the LME is expected to work well.     Table 2. Integrated cross-sections at √ S = 13 TeV, using the input parameters of section 2.2.3 and µ = M ZZ /2. Uncertainties correspond to scale variation as described in the text and, for NLO results that include massive quarks, an estimate of the limitations of the LME. The K-factor is computed using only the central result.   Table 3. Cross-sections at √ S = 13 TeV in the region defined by M ZZ > 300 GeV, using the input parameters of section 2.2.3. Uncertainties correspond to scale variation as described in the text and, for NLO results that include massive quarks, an estimate of the limitations of the LME. The K-factor is computed using only the central result.
by a Higgs boson, with H → ZZ, and by continuum contributions in which the Z bosons couple through loops of quarks. We have considered contributions up to the two-loop level, corresponding to NLO corrections, for the Higgs diagrams alone and also for the interference between the two sets of diagrams.
In the continuum contribution the two-loop corrections containing loops of massless quarks are known and we have reproduced results from the literature. Our treatment of the massive quark loops is based on a large-mass expansion up to order 1/m 12 , that is extended to the high-mass region by using a combination of conformal mapping and Padé approximation. This procedure was shown to provide an excellent approximation of the Higgs contribution alone, where the exact result is known. Additionally, applying the largemass expansion in combination with the conformal mapping and the Padé approximation to the gg → ZZ amplitudes is obviously not limited to the interference calculation alone. The same procedure can also be applied to the virtual two-loop gg → ZZ amplitude including its full tensor structure. It might be desirable to apply the presented procedure also to the Higgs boson pair-production process, because the latter offers identical kinematics. Comparing those results with the recently published results including the full top mass effects [72] could lead to interesting insights concerning the error estimate of the used approximation. However, this is kept as future work.
We have used our calculation to provide theoretical predictions for the impact of the interference contribution on the invariant mass distribution of Z-boson pairs at the 13 TeV LHC. In the high-mass region we have shown that the impact of the NLO corrections to the interference are practically identical to those for Higgs production alone. This explicit calculation justifies using a procedure for estimating the number of off-shell events due to the interference by rescaling the LO prediction by the on-shell K-factor.

JHEP08(2016)011
Comparing each order in α S yields the system of differential equations ⇒ α S 4π F (1) (m, µ) − F (1) (m, µ) = 0 (B.12) (B.13) Solving the homogeneous differential equations for the leading-and next-to-leading-order finite remainder results in  Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.