Production of two Z-bosons in gluon fusion in the heavy top quark approximation

We compute QCD radiative corrections to the continuum production of a pair of Z-bosons in the annihilation of two gluons. We only consider the contribution of the top quark loops and we treat them assuming that $m_t$ is much larger than any other kinematic invariant in the problem. We estimate the QCD corrections to $pp \to ZZ$ using the first non-trivial term in the expansion in the inverse top quark mass and we compare them to QCD corrections of the signal process, $pp \to H \to ZZ$.


I. INTRODUCTION
Production of pairs of vector bosons in proton collisions is one of the most interesting processes studied at the LHC at the end of the Run I [1][2][3]. Indeed, pp → ZZ, pp → W + W − , and pp → γγ, play an important role in Higgs boson physics, provide stringent tests of the Standard Model and give constraints on anomalous electroweak triple gauge boson couplings. In the case of Higgs physics, such processes are essential for understanding backgrounds to Higgs boson signals, for constraining anomalous Higgs boson couplings, for measuring the quantum numbers of the Higgs boson and for studying the Higgs boson width, see e.g. Refs. [4][5][6][7].
Production of electroweak gauge boson pairs occurs mainly due to quark-antiquark annihilation. This contribution is known through next-to-next-to-leading order (NNLO) in perturbative QCD [8,9]. However, as was pointed out in Refs. [10][11][12], there is a sizeable contribution of the gluon annihilation channel gg → V 1 V 2 whose significance depends on the selection cuts. For example [13], agressive cuts applied to pp → W + W − to separate the Higgs boson signal from the continuum background can increase the fraction of gluon fusion events in the background sample to O(30) percent. Since gg → V 1 V 2 is the one-loop process and since production of electroweak boson pairs at leading order occurs only in the qq channel, the gluon fusion contribution to pp → V 1 V 2 through NNLO only needs to be known at the leading, one-loop, approximation. Thus, all existing numerical estimates of the significance of the gluon fusion mechanism in weak boson pair production ignore radiative corrections to gg → ZZ that are, potentially, quite large [14]. The need to have a more accurate estimate of QCD corrections to gluon fusion processes was strongly emphasized in Ref. [7], in the context of the Higgs width and generic off-shell measurements [15][16][17][18].
The largest contribution to gg → V 1 V 2 comes from quarks of the first two generations that can be taken to be massless (for a recent discussion, see Ref. [17]). The contribution of the third generation is, in general smaller. For example, in the case of W + W − production it is known that the third generation changes the gg → V 1 V 2 production cross-section by O(10) percent at the 13 TeV LHC [17]. Since gluon fusion contributes O(5) percent to the pp → W + W − cross-section, the impact of top quark loops on the cross-section is marginal. On the other hand, studies of off-shell Higgs boson production may be senstitive to the third generation of quarks and, especially, to massive top quark loops. Of particular concern in this context is the interference of gg → V 1 V 2 and gg → H * → V 1 V 2 amplitudes, as discussed recently in Refs. [17,18].
The recent progress in calculating two-loop integrals with two massless and two massive external lines [19][20][21][22] enables computation of scattering amplitudes and, eventually, QCD corrections to the production of pairs of vector bosons in gluon fusion through loops of massless quarks. A similar progress towards computing gg → V 1 V 2 contributions mediated by loops of massive quarks is very desirable but, probably, it will not be immediate since two-loop computations of four-point functions with internal massive lines are beyond existing technical capabilities.
In this situation, it is useful to think about alternative approaches that will allow an estimate of QCD radiative corrections to gluon fusion processes mediated by heavy quark loops. A practical opportunity is provided by the expansion of amplitudes in the inverse quark mass. Indeed, this approach reduces the calculation of the oneloop gg → ZZ amplitude with massive internal particles to the calculation of tadpole diagrams which makes generalization to higher-order corrections relatively straightforward. 1 While the expansion of cross-sections in 1/m t cannot be fully justified, particularly for large invariant masses of Z-pairs, we have significant evidence that such computations do provide a reasonable estimate of the size of QCD corrections. Indeed, this is an approach that is taken in calculations of single- [23] and double-Higgs [24] production at the LHC where exact one-loop computations supplemented with QCD corrections calculated in the m t → ∞ approximation are believed to provide reasonably accurate descriptions of these processes for realistic values of top quark and Higgs boson masses. It is clear that a similar approach should be applicable to the production of pairs of Z-bosons in gluon fusion through the top quark loop. In fact, the m t → ∞ approximation for gg → ZZ should work better than for the case of Higgs pair production since 2m Z is smaller than 2m H . The goal of this paper is to make the first step towards estimating the NLO QCD correction to the production of Z-boson pairs in gluon fusion. To this end, we take continuum production of two on-shell Z-bosons in gluon fusion through the top quark loop and compute the NLO QCD corrections to it in the heavy top approximation. This allows us to compare, for the first time, the QCD corrections to the "background" gg → ZZ and the signal gg → H → ZZ processes. We find that the corrections to the two processes are indeed similar, in accord with the arguments in Ref. [14].
It should be clear from the previous discussion that computation of QCD corrections to the total cross-section is just one of many interesting physics questions, including interference with the Higgs signal on and off the mass shell, combination of light and heavy quark contirbutions to the gg → ZZ amplitude, estimates of 1/m t corrections to cross-sections etc., that can be discussed in the context of vector boson pair production in gluon fusion, once the two-loop amplitude for gg → ZZ becomes available. We plan to address these questions in the near future.
The paper is organized as follows. In Section II, we describe the general set up of the computation and present the analytic result for the two-loop amplitude gg → ZZ in the large-m t approximation. In Section III, we derive the analytic formulas for gg → ZZ partonic crosssections. In Section IV we discuss numerical results. We present our conclusions in Section V.

II. THE SET UP OF THE COMPUTATION
We consider the process g(p 1 ) + g(p 2 ) → Z(p 3 ) + Z(p 4 ) in a theory where Z-bosons only couple to top quarks. 2 (top) and massless (bottom) quarks. 2 Such a theory is anomalous and, in principle, one should carefully consider diagrams where a Z-boson couples to two gluons. Given the order in the 1/mt expansion that we work to in this paper, Feynman diagrams where each Z independently couples to gluon pairs do not contribute. Figure 1: Representative two-loop diagrams that describe production of Z-boson pairs in gluon fusion.
Contributions of massless quarks are not included in the computation except in the running of the coupling constant where the complete β-function is employed. The coupling of top quarks to Z-bosons is given by a linear combination of vector and axial couplings The scattering amplitude for gg → ZZ can be written as a sum of axial, vector and mixed terms where a 1,2 are the color indices of the colliding gluons, and α s (µ) is the MS QCD coupling constant in the theory with five active flavors. 3 We note that thanks to charge parity conservation, the axial-vector term vanishes, i.e. A av = 0. The remaining two terms -axial-axial A aa and vector-vector A vv -do not vanish but they behave differently under the 1/m t expansion.
Indeed, consider a vector-current interaction of Z-bosons with top quarks. The gg → ZZ amplitude behaves as A vv ∼ s 2 /m 4 t . This is a direct consequence of the vector current conservation which requires that, in the expression for the amplitude, each polarization vector for either a gluon or an electroweak gauge boson is accompanied by its momentum. A similar suppression in the QED case is familiar in the context of Euler-Heisenberg Lagrangian.
However, if the interaction of Z-bosons with top quarks is mediated by the axial current, the situation is different since the axial current is not conserved. As a consequence, the scattering amplitude can involve two polarization vectors of the Z-bosons without corresponding momenta while the gluon polarization vectors should still be accompanied by their momenta to satisfy the vector current conservation constraint. Hence, we expect that the axial amplitude A aa behaves as A aa ∼ s/m 2 t and, therefore, exhibits weaker suppression in the m t → ∞ limit compared to A vv . Since our goal in this paper is to study the leading term of the A gg→ZZ amplitude in the m t → ∞ expansion, we conclude that we only need to study terms induced by the axial coupling of the Zbosons to top quarks.
The production of Z-boson pairs in gluon fusion is a loopinduced process. There are eight one-loop and nintythree two-loop diagrams that contribute to gg → ZZ. Some examples are shown in Fig. 1. We compute these diagrams using asymptotic expansions in the inverse top quark mass [25]. The essense of this procedure is that the loop momenta in each of the Feynman diagrams are separated into soft l ∼ p 1,..4 and hard l ∼ m t . All possible assignments must be considered. The integrand of a Feynman diagram is then Taylor expanded in all quantities that are considered small. Upon such an expansion, computation of Feynman diagrams significantly simplifies. Consider one-loop diagrams as an example. In this case the momentum can only be hard, l ∼ m t , and so integrands for all diagrams are expanded in Taylor series in their external momenta. 4 All one-loop integrals then become vacuum tadpole integrals and it is straightforward to evaluate them.
The situation with two-loop integrals is similar although somewhat more involved. Indeed, in this case two momentum configurations are possible: either both loop momenta are hard or one of the loop momenta is hard and the other one is soft. If both loop momenta are hard, the calculation is reduced to the calculation of two-loop vacuum tadpole diagrams. If one of the loop momenta is soft and the other one is hard, the diagram factorizes into a product of one-loop integrals, the most complicated of which is a three-point function with all internal and two external lines massless.
We will now present our results for the gg → ZZ amplitude. We write it as an expansion in the strong coupling constant To emphasize constraints on the amplitude that follow from gauge invariance, we introduce the Fourier transform of the field-strength tensor for each of the gluons The one-loop amplitude reads 4 If the loop momentum is assumed to be soft, each propagator is expanded in l/mt generating scaleless integrals. where and ǫ 3,4 are the polarization vectors of the two Z-bosons. We emphasize that the dependence of the amplitude on the dimensional regularization parameter ǫ in Eq.(5) is exact.
In addition to virtual corrections, we require an amplitude for the real emission process, g(p 1 )g(p 2 ) → Z(p 3 ) + Z(p 4 ) + g(p 5 ). To order 1/m 2 t the corresponding amplitude can, in principle, be obtained from the amplitude of the one-loop scattering process in Eq.(5), if the latter is written as a term in an effective Lagrangian, and then used to generate amplitudes with additional gluons in the final state. However, it is also convenient to apply the asymptotic expansion procedure to the computation of the relevant diagrams since this approach can be used to obtain the amplitude for gg → ZZ + g beyond the leading order in 1/m t .
We use the second approach to compute the gg → ZZ +g amplitude. There are fifty diagrams that contribute to this process and we compute the relevant diagrams using the 1/m t expansion. The technical details of the calculation are identical to the calculation of the scattering amplitude for the gg → ZZ process and we do not repeat it here. Unfortunately, the expresssion for the amplitude appears to be too complex to be presented here.
To calculate the production cross-section, we square the elastic and the inelastic scattering amplitudes and integrate them over the corresponding phase-spaces. In order to make the cross-section finite, we need to remove collinear singularities by performing renormalization of parton distribution functions. All of these steps are relatively standard and well-known; for this reason we refrain from describing them in detail.

III. PRODUCTION CROSS-SECTION
We are now in position to present results for the gluon fusion contribution to the production cross-section pp → ZZ. As explained previously, we only consider loops of top quarks and we work to leading order in the 1/m t expansion. We take the invaraint mass of the Z-boson pair to be q 2 and write the differential cross-section as a convolution of the partonic production cross-section and the parton distribution functions In Eq. (8), we used the following notation: f g (x 1,2 ) are the gluon parton distribution functions, τ = q 2 /S hadr and S hadr is the hadronic center-of-mass energy squared. We note that dependencies on the renormalization and factorization scales in Eq.(8) are suppressed. In what follows, we take the factorization and the renormalization scales to be equal.
It is conventional to parametrize the partonic crosssection as where and G(z, q 2 ) can be written as series in the strong coupling constant. To present it, we introduce a parameter r defined as r = q 2 /(4m 2 Z ). We find where D i (z) = ln(1 − z) i /(1 − z) + are the different plus-distribution functions and We note that ∆ 0 has a strong dependence on q 2 . The leading growth caused by the O(r 2 ) ∼ q 4 /m 4 Z term in Eq. (12) is the consequence of the fact that pairs of longitudinal bosons can be produced. It is this growth that should, eventually, get tamed by the destructive interference of gg → ZZ and gg → H * → ZZ amplitudes.
The contributions of hard emissions, not proportional to the leading order cross-section read where ω(z) = 1 − z + z 2 and κ(z) = 1 + z 2 .

IV. NUMERICAL RESULTS
We have implemented the above formulas in a numerical Fortran program that allows us to compute QCD corrections to the top quark loop contribution to the gluon fusion process pp → ZZ as a function of the invariant mass of the Z-bosons, q 2 . We employ NNPDF3.0 parton distribution functions [26] and use leading order parton distributions to compute the production cross-section at leading (one-loop) approximation and next-to-leading order parton distirbutions to calculate it in the two-loop Figure 3: LO pp → ZZ production cross-section ( gluon fusion through top loop only) dσ/dq 2 in fb/GeV 2 , as a function of the invariant mass squared of the Z-boson pair, q 2 , in GeV 2 with a top mass of 173 GeV. We compare our cross-section, which is valid in the mt → ∞ limit, with the one implemented in MCFM, which has exact mt dependence. The dots and squares correspond to the results from MCFM and this paper respectively. We set the renormalization scale and the factorization scale to 200 GeV. The difference between the values ranges from ∼ 20% to ∼ 220% for the values of q 2 considered. The inset shows the MCFM cross-section computed for a larger range of the invariant masses.
approximation. We set the renormalization and factorization scales equal to each other.
To assess the magnitude of QCD corrections, in the main plot of Fig. 2 we show the K-factor defined as the ratio of NLO and LO cross-sections, depending on the invariant mass of the Z-boson pair. We find that the K-factor is a slowly rising function of q 2 and that K ∼ 1.5 − 1.8 for µ = q/2 for the invariant masses considered. The NLO QCD corrections to the gg → ZZ process are therefore similar to what has been observed for other processes where gluons annihilate into colorless final states. As an illustration, we compare the above results with the NLO K-factors for Higgs boson production pp → H, shown in the inset of Fig. 2. We take the Higgs bosons mass to be equal to the invariant mass of the Z-boson pair. The K-factors disagree by about 10 − 15 percent at low values of q 2 and agree almost perfectly at high(er) values of q 2 . This is in accord with the suggestion of Ref. [14] where it was proposed to employ the signal K-factor for the description of the complete process gg → ZZ including the continuum contribution.
Finally, it is interesting to assess how the cross-section expanded in 1/m t compares with the exact result. To this end, we show in Fig. 3 the leading order contribution to gg → ZZ with exact dependence on m t and the m t → ∞ limit. The exact result is obtained using the program MCFM [27]. We see that near threshold values of q 2 , the two cross-sections are similar, within ∼ 20%. At larger values of q 2 the predictions start to diverge as 1/m t supressed terms become more important. As a check of our LO cross-section we performed a similar comparison with MCFM using a top mass of m t = 400 GeV to simulate the m t → ∞ limit. In this case, our calculation is within ∼ 5% of the MCFM predictions for values of q 2 between (200 GeV) 2 and (400 GeV) 2 . The inset in Fig. 3 shows the top quark loop contribution to the pp → ZZ cross-section with full mass dependence, as obtained with MCFM. The cross-section peaks slightly above 400 GeV and then starts to decrease. Our K-factor calculation is valid to the left of the peak, where the cross-section exhibits rapid growth. However, it can be extended beyond that by reweighting the exact gg → ZZ leading order partonic cross-section with K-factors computed in m t → ∞ limit.

V. CONCLUSIONS
In this paper we studied QCD corrections to the production of a pair of Z-bosons in gluon fusion through loops of massive top quarks. This process occurs at one-loop and belongs to the class of processes where two gluons annihilate into a colorless final state. Similar to other processes of this type, such as gg → H and gg → HH, we find large, O(50 − 100) percent radiative corrections. Radiative corrections of this magnitude suggest that the significance of gluon fusion is, perhaps, underestimated by existing NNLO QCD computations of vector boson pair production in proton collisions.
There are several avenues that are interesting to explore as the direct continuation of this work. First, a straightforward extension of this calculation should allow us to compute QCD radiative corrections to the interference of gg → ZZ and gg → H → ZZ amplitudes, both on and off the mass shell of the Higgs boson. Although such a computation will, at the moment, be restricted to top quark contributions to gg → ZZ, it will already give us important information on whether or not the radiative corrections to gg → H → ZZ, gg → ZZ and the interference are related.
Second, it will be interesting to extend our calculation to include higher powers in the expansion of s/m 2 t , to estimate the impact of mass suppressed effects on QCD radiative corrections. In addition, as we explained at the beginning of the paper, the effects of the vector coupling of Z-bosons to top quarks do not appear at leading order in the 1/m t expansion which means that it is important to go one order higher in 1/m t to fully incorporate physics of Ztt interactions into the description of the process. Of course, it is to be expected that since the vector coupling of Z-bosons to top quarks is almost three times smaller than the axial coupling, the inclusion of the vector coupling should only lead to small changes in the cross-section.
Third, it is interesting to incorporate decays of Z-bosons, off-shell effects and realistic selection criteria into our calculation. This should, in principle, be straightforward since the primary objects that we compute are the scattering amplitudes for both gg → ZZ * and gg → ZZ * + g processes.
Finally, it is important to combine contributions of massless quarks and the top quark to gg → ZZ amplitudes in higher orders of QCD. Since the two-loop amplitudes for gg → ZZ with massless intermediate quarks are within reach [8,28], it should be relatively straightforward to incorporate both massless and massive quark loops into the description of gluon fusion contribuitions to Z-boson pair production.
To conclude, we described the calculation of NLO QCD corrections to continuum production of Z-boson pairs in gluon fusion. The results of this computation present first direct evidenece that the gluon fusion production cross-section of Z-bosons receives large O(100%) QCD radiative corrections. Our result should encourage further studies of QCD radiative corrections to weak boson pair production in gluon fusion processes, mediated by massless and massive quark loops.