NNLO predictions for Z-boson pair production at the LHC

We present a calculation of the NNLO QCD corrections to Z-boson pair production at hadron colliders, based on the N-jettiness method for the real radiation parts. We discuss the size and shape of the perturbative corrections along with their associated scale uncertainties and compare our results to recent LHC data at $\sqrt{s}=13$ TeV.


Introduction
The pair production of Z-bosons at the LHC is an important process to test the electroweak sector of the Standard Model (SM). It is sensitive to anomalous gauge boson couplings and constitutes an irreducible background to the production of a Higgs boson decaying into vector bosons and to New Physics searches.
calculation of the NNLO corrections to on-shell Z-boson pairs using a different method, based on N -jettiness subtraction [44,45]. The effect of massive quark loops has been estimated to be at the level of a permille contribution to the total cross section in Ref. [31]. However, calculations performed in an s/m 2 t expansion framework [46,47] indicate that the contributions may be larger, and they certainly will be important in the region of large values of the 4-lepton invariant mass m 4l , which is sensitive to the coupling of the longitudinal Z-boson components to the top quarks loops.

Details of the calculation
The NNLO computation requires the evaluation of the tree-level scattering amplitudes with two additional partons (double-real (RR) contribution), of the one-loop amplitudes with one additional parton (real-virtual (RV) contribution) and the two-loop corrections to the Born process (double-virtual (VV) contribution). In this way we systematically combine all the amplitudes containing two additional powers in the strong coupling constant with respect to the Born process such that the final result is NNLO accurate in perturbation theory. In Table 1 we list the matrix elements for ZZ production at NNLO.
Although the sum of virtual and real corrections yields a finite result, the individual contributions contain singularities of infrared (IR) and ultraviolet nature, such that a direct numerical evaluation is not possible. Virtual and real corrections come from phase space integrals of different multiplicity; therefore a framework to combine them must be such that the divergent regions in the real-radiation contribution (corresponding to soft and collinear emissions which map to configurations with one or two particles less, and therefore are degenerate with the virtual contribution) can be extracted and cancelled with the singularities of the virtual matrix elements.
In this work we employ the N -jettiness subtraction scheme [44,45,48,49] to perform the evaluation of the NNLO cross section. We begin by reviewing the definition of the N -jettiness variable introduced in Refs. [50,51], 1) where N denotes the number of jets desired in the final state and the sum runs over all QCD radiated particles. In Eq. (2.1) the q a , q b and q 1 , . . . , q N are a fixed set of massless reference momenta for the two beam jets and the N observed jets, the p k are the parton momenta, and the dimensionful parameter Q 2 is the hard interaction scale. For the specific case of a colourless diboson system in the final state, Eq. (2.1) reduces to the 0-jettiness or beam thrust which in the leptonic frame reads [50,52], where n a = (1, 0, 0, 1) and n b = (1, 0, 0, −1) define the beam axis and the p k are defined in the hadronic centre-of-mass frame. In the context of N -jettiness subtractions, taking into account the boost with rapidity Y ZZ of the Born system ensures that the power corrections are independent of Y ZZ [52].
Looking at the definition of 0-jettiness in Eq. (2.2) we can observe that T 0 → 0 in the limit where the QCD emission p k is soft or collinear to the initial state. For this reason values of T 0 close to zero indicate a final state containing the ZZ pair and only IR (soft and collinear) emissions. In this way the N -jettiness variable can be used as a slicing parameter in any real-radiation phase space integral to separate infrared singular regions from hard and resolved configurations. In that sense the approach extends the slicing methods developed in the early 90's to compute higher-order corrections at NLO [53] to NNLO.
To proceed we employ a T cut 0 in the real-radiation NNLO phase space and split the cross section into regions above and below T cut 0 [44,45,49], In Eq. (2.3) we have abbreviated θ < 0 = θ(T cut 0 − T 0 ) and θ > N = θ(T 0 − T cut 0 ), and have suppressed any (infrared-safe) measurement function under the phase space integral. The first three terms in this expression all have T 0 < T cut 0 , and are collectively denoted as σ N N LO (T 0 < T cut 0 ), while the remaining two terms, with T 0 > T cut 0 , are collectively denoted as σ N N LO (T 0 > T cut 0 ). Contributions with Born-level kinematics necessarily have T 0 = 0. 1 Contributions with T 0 > T cut 0 necessarily contain one or more well separated hadronic energy deposits and thus reproduce the ZZ+jet cross section at NLO. The contributions with T 0 < T cut 0 correspond to the limit of the ZZ+jet NLO cross section where the jet is unresolved. The key advantage that allows the computation of the cross section at NNLO below T cut 0 is the fact that in the limit where all QCD emission is soft or collinear, the cross section can be approximately computed using the machinery of Soft-Collinear Effective Theory (SCET) [56]. In particular, the existence of a factorization theorem that gives an all-orders description of N -jettiness for small T N less than some value T cut N allows the cross section to be written in the schematic form, where H describes the effect of hard radiation from the purely virtual corrections to the process, B encodes the effect of radiation collinear to one of the two initial beam directions, S describes soft radiation and J n contains the radiation collinear to hard final-state jets. The ellipsis denote power-suppressed terms which become negligible for T N Q. We have expanded the formula in Eq. (2.4) to second order in the strong coupling constant to obtain the σ N N LO (T 0 < T cut 0 ) cross section for ZZ production at hadron colliders. In particular this includes contributions from the universal quark beam function at two loops [57] and the 0-jettiness soft function at two-loops [58,59]. The process dependent hard function has been extracted from the two-loop amplitude computed in Ref. [28] via an interface to the program qqvvamp. We do not include massive top-quark loops in the qqZZ two-loop amplitude. Using N f = 5 therefore introduces the chiral anomaly stemming from subdiagrams where one Z-boson and two gluons are attached to a b-quark triangle. However, we neglect this anomalous contribution in our calculation as the anomaly must cancel once the top quark loops are included, following the same strategy as advocated in Ref. [33].
In SCET renormalised form, the IR finite one-and two-loop amplitudes can be written at renormalisation scale µ 2 as [60,61] Ω (1),finite N = Ω (1) where the N -jettiness subtraction operators are defined by 2 The subtraction operators are identical to those given for diphoton production in Appendix A, Eq. (A.17) of Ref. [62], but these formulae appear to contain two typographical errors. Specifically, we find that Γ 0 should be multiplied to β0 in the second term of Eq. (2.6) and that the last term has a factor 16 in the denominator. and constants In Ref. [28] the finite remainder of the one-and two-loop form factors for vector boson pair production are provided in the q T - [63] and  subtraction schemes. By comparing the definition of the subtraction schemes, the form factors in the N -jettiness scheme can be derived from those of the q T -scheme according to, with the scheme conversion coefficients given by where, for brevity, we have set µ 2 = s. Finally we have obtained the σ N N LO (T 0 > T cut 0 ) contribution of the ZZ NNLO cross section using the tree level matrix elements from VBFNLO [65,66] for the double-real emission phase space integral cross-checked with MadGraph5 [67], while the one-loop amplitudes for the real-virtual phase space were generated with GoSam [68,69] and cross-checked with OpenLoops [70]. GoSam uses QGRAF [71], FORM [72] and Spinney [73] for the generation of the Feynman diagrams, and offers a choice from Samurai [74], golem95C [75][76][77] and Ninja [78,79] for the reduction. At run time the amplitudes were computed using Ninja, which calls OneLOop [80] for the master integrals, and rescued using an implementation of Ninja in quadruple precision for unstable phase space points. We also include the loop induced one-loop squared corrections in the gg → ZZ channel, which are formally of NNLO accuracy, keeping full dependence on the top quark mass and on the Higgs mediated contributions using GoSam.

Discussion of the IR subtraction procedure
Before we present our numerical results for ZZ production at hadron colliders we would like to make a few remarks concerning the IR subtraction scheme employed for this calculation. As mentioned in the previous section, the N -jettiness subtraction scheme is a non-local subtraction scheme. In local subtraction schemes the IR divergent phase space integrals are regulated by the introduction of suitable IR real-radiation counterterms that satisfy two basic properties, • reproduce locally, for each phase space point in a singular region, the physical IR divergent soft and collinear limits of the matrix elements of the process under consideration; • be simple enough to allow their analytic integration and generate a local and pointwise analytic pole cancellation between the explicit 1/ -poles of the virtual corrections and the 1/ -poles of the integrated real-radiation counterterms.
Some flexibility however exists on how the singular limits are locally subtracted and how the real-radiation phase space is parametrised. For antenna subtraction [81][82][83] physical matrix elements with three partons [84] at tree-level and one-loop suffice to reproduce single unresolved limits in QCD amplitudes, while four parton antennae [85,86] subtract doubly unresolved configurations at NNLO [87][88][89][90]. Examples of local IR subtraction schemes which employ a structured decomposition of the real-radiation phase space based on singular IR limits of QCD amplitudes [91][92][93][94] have also been developed and applied to specific NNLO calculations. For the specific case of colourless systems in the initial state local subtractions have been developed in Refs. [95,96]. The extension of the N -jettiness method towards local subtractions has also been conceptually discussed in Ref. [45].
On the other hand, non-local IR subtraction schemes regulate the singularities of the real-radiation phase space integrals by the introduction of a suitable variable (N -jettiness or the transverse momentum q T of a colourless system for q T -subtraction [43] 3 ) that regulates the phase space integration by separating IR divergent regions from hard and resolved ones according to Eq. (2.3). In this way contributions to the cross section for T 0 above and below T cut 0 are separately logarithmically divergent. However, in the sum all the logarithmic dependence on T cut 0 should cancel, provided the value of T cut 0 employed in the phase space integration is small enough such that the SCET approximation to the cross section is valid. In particular, for each 1/ -IR pole in dimensional regularisation there is a corresponding logarithmically divergent coefficient predicted from SCET in the T 0 → 0 limit, according to 1 n ∼ log n T 0 µ . (2.10) The infrared pole cancellation in this case is observed through the cancellation between the universal and analytically known terms predicted by SCET, integrated over the Born phase space, and the Monte Carlo integration over the real-radiation phase space of the real-emisson matrix elements for small T 0 . The method of N -jettiness meanwhile has been applied successfully to various processes calculated at NNLO [44,48,49,62,[98][99][100][101][102][103][104].
Due to the non-local IR subtraction method employed in our calculation of the NNLO corrections we found it necessary to do the following optimisations at the Monte Carlo integration level in order to observe the independence of our results on the choice of the slicing parameter value T cut 0 : • introduce a phase space generator where the 0-jettiness variable T 0 is directly sampled by VEGAS. This can be achieved by applying the transformation p ±,i = E i ± p z,i , dE i dp z,i = 1 2 dp +,i dp −,i (2.11) to the integration over the real radiation phase space. 4 With the momenta defined in the center-of-mass system of the Z-bosons, the 0-jettiness as defined in Eq. (2.2) then reads T 0 = min(p +,1 , p −,1 ) + min(p +,2 , p −,2 ). (2.12) With further transformations, where the regions with p +,i < p −,i and p +,i > p −,i have to be distinguished, it is then possible to directly sample T 0 , followed by the sampling of min(p +,1 , p −,1 ), which also fixes the value of min(p +,2 , p −,2 ). Afterwards, the two missing values of p ±,i can be sampled.
• have a fast implementation of the double-real and real-virtual matrix elements for ZZ-production which is stable in the multiple soft and collinear limits.
The first optimisation ensures that the real-radiation phase space generator properly samples the phase space boundaries determined by the choice of slicing parameter T cut 0 and that the phase space integral converges. The second optimisation ensures that the matrix elements are fast and stable enough to be integrated near the singular IR limits which give the bulk of the cross section for the (T 0 > T cut 0 ) phase space integrals when T cut 0 is small. This is in contrast with a local subtraction scheme, where the real radiation subtraction terms ensure that the integrand vanishes as we approach a singular region.

Results
Our numerical studies for proton-proton collisions at centre-of-mass energy √ s =13 TeV are for on-shell Z-boson pair production. We use the MSTW2008 [105] and NNPDF-3.0 [106] sets of parton distribution functions via the LHAPDF [107] interface, with densities and α s evaluated at each corresponding order (i.e. we use (n+1)-loop α s at N n LO, with n = 0, 1, 2) and we consider N f = 5 massless quark flavours. The default renormalisation (µ R ) and factorisation (µ F ) scales are set to µ R = µ F = m Z . We use the G µ EW scheme where the EW input parameters have been set to G F = 1. We show in Fig. 1 the NLO and NNLO coefficients of the ZZ cross section as a function of T cut 0 . On the left-hand side we compare the NLO results obtained using either antenna subtraction or N -jettiness subtraction and observe full agreement in the evaluation of the NLO corrections ∆σ NLO . We present separate results for the phase space integration with real emission kinematics (2 → 3) and Born-like kinematics (2 → 2). Obviously, using a local subtraction scheme (antenna subtraction), all phase space integrals contributing to ∆σ NLO are independent of the choice of T cut 0 . Moreover, the bulk of the NLO coefficient comes from the 2 → 2 phase space integral which determines where more statistical precision is needed to obtain an accurate result. On the other hand, using N -jettiness we observe that both phase space integrals are separately double-logarithmically divergent and therefore a very good numerical precision is needed for both contributions to improve the accuracy of the final result. In Fig. 1b we present the NNLO coefficient of the ZZ cross section as a function of T cut 0 . In this case we observe that the phase space integrals for the contributions σ N N LO (T 0 < T cut 0 ) and σ N N LO (T 0 > T cut 0 ) are logarithmically divergent to the fourth power in log T cut 0 , and for typical values of T cut 0 in the range 10 −2 − 10 −3 GeV need to be known with better than permille level accuracy to achieve an accurate determination of the NNLO coefficient.
In order to study in more detail the independence of the NNLO coefficient on the choice of slicing parameter T cut 0 we present in Fig. 2 on a smaller scale the NNLO coefficient after combining the contributions for σ N N LO (T 0 < T cut 0 ) and σ N N LO (T 0 > T cut 0 ) as black data

ATLAS [7]
17.3 ± 0.6(stat.) ± 0.5(syst.) ± 0.6(lumi.) CMS [8] 17.2 ± 0.5(stat.) ± 0.7(syst.) ± 0.4(theo.) ± 0.4(lumi.) Table 2: Inclusive cross section for ZZ production at the LHC run II √ s =13 TeV at LO, NLO and NNLO with µ R = µ F = m Z , together with the measurements from ATLAS [7] and CMS [8]. Uncertainties in the theory calculation at each order are obtained by varying the renormalisation and factorisation scales in the range 0.5m Z < µ R , µ F < 2m Z with the constraint 0.5 < µ F /µ R < 2. Uncertainties in the experimental measurements denote absolute statistical, systematic and luminosity uncertainties. with respect to other processes [44,49,102] seems to indicate that for ZZ production their contribution is small. Nonetheless the leading power correction can be modeled after integration over the final-state phase space as [52,108,109] where Q is an appropriate hard scale of the process and c 2 , c 3 are unknown constants. We have performed a fit of the results of our Monte-Carlo runs to this functional form of the N -jettiness NNLO coefficient for ZZ and show the resulting fit as a black dotted line in Fig. 2. The fit allows us to numerically extract the value of the NNLO coefficient in the limit where T 0 → 0. This value can be compared to the reconstructed NNLO coefficient for ZZ production obtained in Ref. [31] 5 , which is shown as a red line. We use the extrapolated value for our result for the ZZ cross section at NNLO shown in Table 2, which is in excellent agreement with the result σ N N LO = 16.91 pb obtained in Ref. [31]. As a consistency check we have also fitted a constant to the plateau region (T cut 0 < 10 −2 GeV or T cut 0 < 10 −1 GeV) and these fits yield compatible results for ∆σ N N LO . Further, we have also fitted the leading power corrections using (3.1) including only results for T cut 0 < 1 GeV. When fitting the leading power corrections with T cut 0 < 1 GeV there is a strong correlation between c 3 and Q as well as c 2 and Q; fixing Q to values in the range 50 − 5000 GeV we obtain compatible results for ∆σ N N LO . Including in the fit results up to T cut 0 < 10 2 GeV, as shown in Fig. 2 , we obtain a stable fit also when Q is treated as a free parameter. 5 The NNLO coefficient was reconstructed by subtracting from the total NNLO ZZ cross section quoted in Table 1 Table 2. In the same Table we present an updated value for the NNLO cross section computed as described in the previous section using the more recently determined NNPDF-3.0 [106] PDF sets and an updated value for the W -boson mass of M W =80.385 GeV; these settings are also used for our phenomenological results in the following. We observe a significant improvement in the agreement with the data after the inclusion of the NNLO corrections.
In order to study in more detail the scale uncertainty of the cross section we present in Fig. 3 the renormalisation and factorisation scale dependence of the ZZ cross section at LO, NLO and NNLO. The figure shows largely non-overlapping scale uncertainty bands which demonstrate that for this process, the scale variations are insufficient to estimate missing higher order terms in the perturbative expansion. This however is not unexpected since ZZ production at the LHC is an electroweak process which exhibits no renormalisation scale dependence at LO. For this reason we obtain large NLO QCD corrections to the cross section which are outside the LO scale band. Moreover, when going from NLO to NNLO, the loop-induced gluon fusion channel gg → ZZ opens up, and due to the large gluon flux it represents a numerically significant contribution. Since this new channel contributes for the first time at NNLO its contribution cannot be captured by the scale variation of the NLO cross section. Therefore, when increasing the perturbative order, we can observe a systematic reduction of the factorisation scale dependence of the cross section (indicated by the thickness of the scale uncertainty band), while there is no significant reduction of the renormalisation scale dependence. To show that this effect can be attributed to the gluon fusion channel opening up at NNLO, we also show the NNLO result excluding this channel, leading to an improved convergence of the perturbative expansion.
The appearance of new channels that open up at NNLO and their importance in the various kinematic regions can be studied by considering differential results. Due to the observed mild power corrections in this process we chose to fix the value of the 0jettiness slicing parameter to T cut 0 = 10 −2 GeV for all our histograms. In Fig. 4 we present the invariant mass of the ZZ system and the average transverse momentum distribution p T,Z of any Z-boson, defined as p T,Z = (|p Z 1 T | + |p Z 2 T |)/2. We also present results for the loop-induced gg → ZZ channel.
In Fig. 4a we show our results for the ZZ invariant mass. In the first and second sub-panels we show the effect of the NLO and NNLO corrections, respectively. We observe in the first sub-panel large NLO QCD corrections which vary between 40% at low m ZZ and 60% at high m ZZ , and change both the shape and normalisation of the predicted cross section with respect to the LO result. Going to NNLO we observe an approximately flat increase of the cross section of about 18% with respect to the NLO result, where approximately 60% of this effect comes from the loop-induced gg → ZZ channel, which is outside the scale uncertainty band of the NLO prediction. Similarly, in the transverse momentum distribution (Fig. 4b), we observe large NLO corrections of approximately 30% at low p T,Z , which can reach almost 100% at high p T,Z . The shape of the NNLO  corrections in the second sub-panel largely follows the contribution of the loop-induced gg → ZZ channel at low p T,Z , and we observe a 30% effect at low p T,Z which decreases to about 18% at high p T,Z . For both distributions we observe good convergence of the perturbative expansion, however the scale uncertainty bands do not overlap between the orders in the perturbative expansion that we have computed. These results show that the inclusion of NNLO effects in ZZ production at the LHC is essential to obtain a reliable theoretical description of this process.

Conclusions
In this work we have calculated the NNLO QCD corrections to on-shell Z-boson pair production. Our calculation of the real emission contributions uses N -jettiness to isolate the infrared divergent contributions. We discussed our setup in some detail, also showing a comparison between results based on antenna subtraction and results based on N -jettiness for the NLO corrections.
After the inclusion of the NNLO correction in the theory prediction we found good agreement with the results of the recent ATLAS and CMS measurements. Due to the fact that the numerically sizeable loop-induced gg → ZZ contribution appears for the first time at NNLO, the scale uncertainties in the NNLO prediction are not reduced with respect to NLO. The NNLO corrections increase the NLO result by about 18%, where almost 60% of this increase stems from the loop-induced gg → ZZ channel. In view of the numerical importance of this channel, it is desirable to add the two-loop diagrams, including massive top quark loops, to this channel, which will be left for a subsequent publication.