QCD corrections to ZZ production in gluon fusion at the LHC

We compute the next-to-leading order QCD corrections to the production of two Z-bosons in the annihilation of two gluons at the LHC. Being enhanced by a large gluon flux, these corrections provide distinct and, potentially, the dominant part of the N$^3$LO QCD contributions to Z-pair production in proton collisions. The $gg \to ZZ$ annihilation is a loop-induced process that receives the dominant contribution from loops of five light quarks, that are included in our computation in the massless approximation. We find that QCD corrections increase the $gg \to ZZ$ production cross section by ${\cal O}(50\%-100\%)$ depending on the values of the renormalization and factorization scales used in the leading order computation, and the collider energy. The large corrections to $gg \to ZZ$ channel increase the $pp \to ZZ$ cross section by about six to eight percent, exceeding the estimated theoretical uncertainty of the recent NNLO QCD calculation.


I. INTRODUCTION
Production of pairs of vector bosons in proton collisions is one of the most interesting processes studied by ATLAS and CMS during the LHC Run I [? ? ? ]. Indeed, pp → ZZ, pp → W + W − , and pp → γγ were instrumental for the discovery of the Higgs boson. As the focus of Higgs physics shifts from the discovery to precision studies of the Higgs boson properties, di-boson production processes become essential for constraining anomalous Higgs boson couplings, for measuring the quantum numbers of the Higgs boson and for studying the Higgs boson width, see Refs. [? ? ? ? ]. Additionally, these processes provide important tests of our understanding of the Standard Model and can be used to constrain anomalous electroweak gauge boson couplings.
Production of electroweak gauge boson pairs occurs mainly due to quark-antiquark annihilation qq → V 1 V 2 . This contribution is known through next-to-next-to-leading order (NNLO) in perturbative QCD [? ? ? ? ? ? ]. However, as was pointed out in Refs. [? ? ? ], there is a sizable contribution from the gluon annihilation channel gg → V 1 V 2 , whose significance depends on the selection cuts. For example, aggressive 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 [? ]. Since gg → V 1 V 2 is a one-loop process and since production of electroweak boson pairs at leading order (LO) occurs only in the qq channel, the gluon fusion contribution to pp → V 1 V 2 through NNLO only needs to be known at leading order, i.e. the 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 [? ]. The need to have an accurate estimate of QCD corrections to gluon fusion processes for the Higgs width [? ? ] and generic off-shell measurements [? ? ? ] was strongly emphasized in Ref. [? ].
In this paper, we will focus on the calculation of the next-to-leading order (NLO) QCD corrections to the gluon fusion contribution to pp → ZZ process. The largest contribution to gg → ZZ comes from quarks of the first two generations; these quarks can be taken to be massless. The situation is more complicated for quarks of the third generation. Ideally, we would like to include the (massless) bottom quark contribution and ignore the contribution of the massive top quark since, at leading order, the top-quark contributions change the cross section by only about 1% (cf. Refs. [? ? ]). 1 We can separate bottom and top contributions everywhere except in triangle diagrams that involve anomalous correlators of vector and axial currents. In these triangle diagrams, when bottom and top contributions are combined, the residual contributions are suppressed by the top quark mass, provided that we can assume it to be larger than any other energy scale in the problem. Unfortunately, in these diagrams top and bottom contributions can not be separated because the resulting theory is anomalous. To deal with this issue, we adopt the following strategy: we include quarks of the first two generations and the b-quark in our calculation in the massless approximation and we neglect all triangle diagrams whose contribution is then naturally associated with the quark contributions to gg → ZZ process. We note that the evaluation of the NLO QCD corrections to top quark mediated contribution to gg → ZZ process is not yet possible because the relevant two-loop amplitudes are not available. However, such contributions were recently studied in Ref. [? ] in the approximation of a very large mass of the top quark. In that calculation quite large QCD corrections were found.
Computing NLO QCD corrections to gg → ZZ process is challenging because it is loop- The paper is organized as follows. In Section ?? we present a brief review of the calculation of the two-loop scattering amplitude for gg → ZZ process. In Section ?? we discuss the calculation of the one-loop helicity amplitudes for gg → ZZg and present numerical results for a kinematic point. In Section ?? we present numerical results for gg → ZZ contribution to pp → ZZ process at 8 and 13 TeV LHC at leading and next-to-leading order in perturbative QCD. We conclude in Section ??.

II. THE TWO-LOOP SCATTERING AMPLITUDES FOR gg → ZZ
We start with a brief discussion of the two-loop scattering amplitudes for gg → ZZ process.
Helicity amplitudes for this process were recently computed in Refs. [? ? ]. In these references, each of the two independent helicity amplitudes for the process gg → ZZ → 4l was written as linear combinations of nine form factors that depend on the Mandelstam invariants of the "prompt" process gg → ZZ and the invariant masses of the two Z bosons.
The form factors are expressed in terms of polylogarithmic functions, including both ordinary and Goncharov polylogarithms.
In this paper we use the results of Ref. [? ] which are implemented in a C++ code that can produce numerical results with arbitrary precision. In order to detect possible numerical instabilities, the code compares numerical evaluations obtained with different (double, quadruple and, if required, arbitrary) precision settings. If the results differ beyond a chosen tolerance, the precision is automatically increased. Of course, switching to arbitrary precision increases the evaluation time substantially. Fortunately, we found that for phenomenologically relevant situations, the number of points where the code switches to arbitrary precision is negligible. Such points originate from kinematic regions where the two Z-bosons have either vanishing kinetic energies or vanishing transverse momenta. The amplitude squared is integrable in both of these regions, but, in practice, it can become numerically unstable. Since the contribution of these regions to the gg → ZZ cross section is relatively small, cutting them away, in principle, leads to an opportunity to perform stable numerical 2 For recent reviews see Refs. [? ? ]. integration of the two-loop virtual correction over the four-lepton phase-space, resorting to quadruple precision only. However, we found that the improvement in performance achieved by cutting away the problematic regions is rather limited, so we used the default arbitrary precision implementation of the two-loop amplitude in practice.
Since the gg → ZZ amplitude is one of the most complicated amplitudes that are currently known analytically, it is interesting to point out that the required evaluation times are acceptable for phenomenological needs. Indeed, calculation of all helicity amplitudes requires about two seconds per phase-space point in quadruple precision and, since the phase-space for gg → ZZ is relatively simple, one does not need excessively large number of points to sample it with good precision.
For further reference we provide numerical results for the finite remainder of the one-and two-loop scattering amplitudes defined in q t -subtraction scheme, see Ref. [? ]. The numerical results are presented for the choice of the renormalization scale µ = √ s, where s is the partonic center-of-mass energy squared. The q t -subtraction scheme [? ] is discussed in detail in Ref. [? ]. We consider the kinematical point with (in GeV units) and define a normalized amplitude through the following equation d σ gg→(Z/γ)(Z/γ)→4l = (N 2 c − 1) 512s × 10 −6 × λ 1 ,λ 2 ,λe,λµ A(1 λ 1 g , 2 λ 2 g ; 3 λe e − , 4 −λe e + , 5 (2) Note that in Eq.(??) all the color factors have been factored out and dLIPS 4 is the standard Lorentz-invariant phase-space of the four final leptons. The color-stripped amplitude admits