Next-to-next-to-leading order event generation for Z boson pair production matched to parton shower

We present a novel next-to-next-to-leading order (NNLO) QCD calculation matched to parton shower for the production of a pair of Z bosons decaying to four massless leptons, pp→ `+`−`′+`′−+X, at the LHC. Spin correlations, interferences and off-shell effects are included throughout. Our result is based on the resummed beam-thrust spectrum, which we evaluate at next-to-next-to-leading-logarithmic (NNLLT0) accuracy for the first time for this process, and makes use of the Geneva Monte Carlo framework for the matching to Pythia8 shower and hadronisation models. We compare our predictions with data from the ATLAS and CMS experiments at 13 TeV, finding a good agreement.


Introduction
Diboson production at the LHC is of paramount importance to the precision study of electroweak (EW) physics in the Standard Model (SM) and beyond, as it directly probes the non-Abelian EW couplings. The four-lepton channel that proceeds predominantly through intermediate Z-boson pairs is of particular importance due to its very clean signature. Consequently, several cross-section measurements and studies of anomalous couplings were carried out by both the ATLAS and CMS collaborations at 7 TeV [1][2][3][4], 8 TeV [5-9] and 13 TeV [10][11][12][13][14][15][16]. Moreover, analyses of the four-lepton channel have helped to constrain the width and couplings of the Higgs boson [17][18][19][20][21][22][23][24][25][26][27][28]. In addition, given that Z resonances are common features in models of new physics, measurements of Z pair production can help to set limits on Beyond the SM scenarios.
In order to take advantage of the increasing precision of experimental data, it is useful to provide predictions in the form of fully exclusive Monte Carlo (MC) event generators. These allow hadron-level events to be produced that can be directly interfaced to detector simulations and experimental analyses.
For the qq → ZZ → 4 channel, the state-of-the-art in this regard is still matching NLO QCD correction to parton shower (NLOPS), as implemented in Powheg [51,52], also including NLO EW effects [53].
Currently, the only available NNLOPS calculation featuring two massive bosons in the final state is for the process pp → W + W − → ν ν [68] via the MiNLO method [56], which made use of a differential reweighting to the Matrix predictions [69] to achieve NNLO accuracy. Naïvely, the complexity of the final state would require a nine-dimensional reweighting, something unfeasible in practical terms. In Ref. [68], the authors were able to circumvent this limitation by rewriting the differential cross section in terms of angular coefficients, which they used to reweight each angular contribution separately. The dimensionality of the reweighting procedure was thus reduced to just three. This relied, however, on the approximation that the vector bosons are close to being on-shell, and so cannot be easily applied to the case of pp → Z/γ * Z/γ * → 4 , given the non-negligible contribution from photon exchange.
Since the Geneva method does not need a multidifferential reweighting to reach NNLO accuracy, we are able to include the full off-shell effects and deliver the first NNLOPS for ZZ production. 1 The process pp → + − + − +X features contributions from channels with very different resonance structures. In order to increase efficiency in event generation in such a situation, it is necessary to make use of a phase space generator which samples these channels separately. To this end, we have built an interface between Geneva and the multi-channel integrator Munich [70] which is completely general and allows for the integration for any SM process.

Process definition
In this Letter, we consider the process pp → + − + − + X, including off-shell effects, Z/γ * interference and spin correlations. The derivation of the Geneva formulae has been presented in several papers, see e.g. [59,65]. Here we content ourselves with presenting the final results for the differential weights of the 0−,1− and 2− jet partonic cross sections.
These are defined in such a way that they correspond to physical and IR-finite events at a given perturbative accuracy, with the condition that IR singularities cancel on an event-by-event basis. The Geneva method achieves this by mapping IR-divergent final states with M partons into IR-finite final states with N jets, with M ≥ N . Events are classified according to the value of N -jet resolution variables T N which partition the phase space into different regions according to the number of resolved emissions. In particular, the Geneva Monte Carlo cross section dσ mc N receives contributions from both N -parton events and Mparton events where the additional emission(s) are below the resolution cut T cut N used to separate resolved and unresolved emissions. The unphysical dependence on the boundaries of this partitioning procedure is removed by requiring that the resolution parameters are resummed at high enough accuracy. For the production of a pair of Z bosons at NNLO accuracy, we need to introduce the 0-jettiness resolution variables T 0 and one-jettiness T 1 defined as with N = 0 or 1, q a , q b representing the beam directions and q k any final state massless direction that minimise T N . These separate the 0-and 1-jet exclusive cross sections from the 2-jet inclusive one.
We find that dσ mc where the B j , V j and W j are the 0-, 1-and 2-loop matrix elements for j QCD partons in the final state.
We have introduced the shorthand notation to indicate that the integration over a region of the Mbody phase space is performed while keeping the N -body phase space and the value of some specific observable O fixed, with N ≤ M . The Θ O (Φ N ) term in the previous equation limits the integration to the phase space points included in the singular contribution for the given observable O. For example, when generating 1-body events we use where the map used by the 1 → 2 splitting has been constructed to preserve T 0 , i.e.
and Θ T (Φ 2 ) defines the projectable region of Φ 2 which can be reached starting from a point in Φ 1 with a specific value of T 0 . The use of a T 0 -preserving mapping is necessary to ensure that the pointwise singular T 0 dependence is alike among all terms in eqs. (3) and (5) and that the cancellation of said singular terms is guaranteed on an event-by-event basis. The expressions in eqs. (4) and (6) encode the nonsingular contributions to the 1-and 2-jet rates which arise from non-projectable configurations below the corresponding cut. This is highlighted by the appearance of the complementary Θ functions, Θ O , which account for any configuration that is not projectable either because it would result in an invalid underlying-Born flavour structure or because it does not satisfy the T 0 -preserving mapping (see also [65]).
The term V C 1 denotes the soft-virtual contribution of a standard NLO local subtraction (in our implementation, we follow the FKS subtraction as detailed in [71]). We have that with C 2 a singular approximation of B 2 : in practice we use the subtraction counterterms which we integrate over the radiation variables dΦ 2 /dΦ C 1 using the singular limit C of the phase space mapping. U 1 is an NLL Sudakov which resums large logarithms of T 1 and U 1 its derivative with respect to T 1 .
The term P(Φ N +1 ) represents a normalised splitting probability which serves to extend the differential dependence of the resummed terms from the N -jet to the (N+1)jet phase space. For example, in eq. (3), the term P(Φ 1 ) makes the resummed spectrum in the first term (which is naturally differential in the Φ 0 variables and T 0 ) differential also in the additional two variables needed to cover the full Φ 1 phase space. These splitting probabilities are normalised, i.e. they satisfy The two extra variables are chosen to be an energy ratio z and an azimuthal angle φ. The functional forms of the P(Φ N +1 ) are based on the Altarelli-Parisi splitting kernels, weighted by PDFs where appropriate.
For the specific details of the implementation of the above formulae we refer the reader to Ref. [59].
The resummed contribution in the previous formulae are obtained from Soft-Collinear Effective Theory (SCET), where a factorisation formula for the production of a colour singlet can be written as The sum in the equation above runs over all possible qq pairs ij = {uū,ūu, dd,dd, . . .}. It also depends on the hard H ij , soft S and beam B i,j functions which describe the square of the hard interaction Wilson coefficients, the soft emissions between external partons and the hard emissions collinear to the beams respectively. In SCET, the resummation of large logarithms is achieved by means of renormalisation group (RG) evolution between the typical energy scale of each component (µ H , µ B and µ S ) and a common scale µ. This proceeds via convolutions of the single scale factors with the evolution functions U i (µ i , µ). The resulting resummed formula for the T 0 spectrum is then given by where the convolutions between the different functions are written in a schematic form. In order to reach NNLL accuracy, we need to know the boundary conditions of the evolution, namely the hard, beam and soft functions up to NNLO accuracy, and the cusp (non-cusp) anomalous dimensions up to three-(two-)loop order. The beam and soft functions at NNLO accuracy are available in the literature [72][73][74] as well as and the cusp (non-cusp) anomalous dimensions up to three-(two-)loop order [75][76][77][78][79].
The two-loop hard function H ij was instead computed starting from the form factors calculated in Refs. [80,81] and using the numerical implementation in the public code VVamp [82]. Finally, we use OpenLoops 2 [83][84][85]. for the calculation of all the remaining matrix elements.
Starting from these accurate parton-level predictions we can interface to the Pythia8 parton shower to produce the high multiplicity final states that can in turn be compared to experimental data. The shower adds extra radiation to the exclusive 0-and 1-jet cross sections and extends the inclusive 2-jet cross section by including higher jet multiplicities. This means that we can expect the shower to modify the distributions of exclusive observables sensitive to the radiation, preserving the leading logarithmic accuracy for observables other than T 0 . At the same time, we have verified that the shower does not modify any distribution which is inclusive over the radiation, as was the case for the Geneva implementation of similar colour-singlet production processes. As a consequence we maintain NNLO accuracy for inclusive observables.
All the details about the interface between Geneva and Pythia8 can be found in section 3 of [59].

Results and comparison to LHC data
In the following, we focus on the process at a hadronic centre-of-mass energy of 13 TeV and we require that the masses of both lepton-antilepton pairs are between 50 and 150 GeV. We use the PDF set NNPDF31_nnlo_as_0118 [86] from LHAPDF6 [87] and set both the renormalisation and factorisation scales to the mass M 4 of the four-lepton system. We choose the resolution cutoffs to be T cut 0 = 1 GeV and T cut 1 = 1 GeV. For the validation, we focus only on the quark-antiquark channels and neglect the loop-induced gluon fusion channel, which only starts appearing in the calculation at NNLO and can therefore be added as a nonsingular fixedorder contribution. At this energy this amounts to ∼ 6% of the total cross section and as such its inclusion will be important when comparing to data.
In Fig. 1 we validate the NNLO accuracy of the Geneva results by comparing with those obtained via an independent NNLO calculation implemented in Matrix. In particular, we show the distributions of the rapidity y 4 , mass M 4 of the four leptons and the transverse momentum p − T of the electron. We observe a good agreement, with the only differences appearing in the shape of the p e − T distribution. This is likely to be a consequence of the additional higher-order effects provided by Geneva, as observed in previous Geneva predictions for other colour-singlet production processes. We observe, however, that in almost all the bins the theoretical uncertainty bands computed by Geneva and Matrix still overlap.
Next we turn on the shower effects by interfacing to Pythia8. In order to maintain a simple analysis routine and focus only on QCD corrections, we do not include QED effects and multiple particle interactions (MPI) in the shower.
In Fig. 2 we present our predictions for the NNLL T0 +NNLO beam-thrust spectrum (partonic result) and study the effect of the shower and hadronisation on the T 0 distribution. In order to highlight the peak and transition regions where the resummation effects are more important, we show the results on a semi-logarithmic scale, which is linear up to the value of 30 GeV and logarithmic beyond. In the first ratio plot, we compare the T 0 distribution before and after the shower, observing that the NNLL T0 +NNLO accuracy reached at parton level is numerically very well preserved. The largest difference is, as expected, in the first bin where the shower generates all the events with T 0 < T cut 0 which were previously integrated out in the inclusive 0-jet cross section. In the second ratio plot we instead show the effect of hadronisation on the showered distribution. Owing to its nonperturbative origin, hadronisation affects mostly the region at small T 0 and becomes less and less important in the tail of the distribution.
Finally, in Fig. 3 we compare our predictions with data obtained at the LHC at 13 TeV both from the ATLAS [11] and CMS [16] experiments. For this comparison we now include the loop-induced gluon fusion channel and MPI effects. The binning of observables and the acceptance cuts applied to the events in each analysis have been modelled according the experimental requirements. They are briefly listed in the figure, but we refer the reader to the original papers for all the details.
In the first row we show the comparison with ATLAS data, obtained with an integrated luminosity of 36.1 fb −1 . We consider the distributions for the transverse momentum of the four leptons p T,4 , and two more inclusive distributions, the absolute value |y 4 | of the rapidity of the four leptons and the absolute value |∆y Z1Z2 | of the difference in rapidity between the two bosons. The agreement with data is reasonably good but the reduced experimental statistics make it difficult to draw more precise conclusions. We only notice a possible tension in the distribution of |y 4 |, where for small values of the rapidity the data seem to be systematically above the Geneva predictions. This behaviour has, however, been observed previously with several other Monte Carlo event generators in the original ATLAS publication [11].
In the second row, we show the comparison with CMS data, obtained with an integrated luminosity of 137 fb −1 . We focus on the normalised distributions for p T,4 , the transverse momentum p T, + − of the vector bosons and the transverse momentum p T, of the leptons. The latter two distributions are obtained averaging over all vector bosons and leptons present in an event. Due to the increased luminosity and the fact that we are comparing normalised distributions, here we observe smaller experimental uncertainties and a better agreement between the data and the Geneva predictions. We highlight the fact that according to the original definition of bins in Ref.
[16], the last bin of each distribution we show must be interpreted as an overflow bin containing the contributions from the lower border up to the maximum possible value for that observable. This justifies for example the peculiar behaviour of the last bin of the p T, distribution.
We conclude by noticing that the transverse momentum p T, + − of the vector bosons shows a sizable difference for values larger than 150 GeV, where EW effects are known to be important [46] and should thus be included to improve the agreement with data.

Conclusions
In this Letter, we have presented the first event generator for the production of a Z boson pair decaying into four leptons at NNLO accuracy matched to the Pythia 8 parton shower.
This was obtained using the resummation of the 0-jettiness T 0 resolution variable at NNLL accuracy, matched to a fixed-order calculation at NNLO precision. This was performed within the Geneva Monte Carlo framework, which allowed us to interface with the Pythia8 parton shower and hadronisation models.
After successfully validating the NNLO accuracy of the results against Matrix, we studied the effect of the shower on the differential distributions, observing that they are numerically small for inclusive quantities, as expected. Consequently, their NNLO accuracy is correctly maintained. We also verified that the hadronisation has a significant impact only for exclusive observables in the region of small T 0 , where the nonperturbative effects are not negligible.
Finally, we compared to LHC data, both from the AT-LAS and CMS experiments. We found good agreement, except for a tension observed in the region of small absolute rapidity of the four lepton system. A similar behaviour had already been observed for other Monte Carlo generators. We also observed a general overestimation of the production rate for vector bosons at large transverse momentum, motivating the need for the inclusion of EW corrections to improve the agreement with data in this region.
Possible future directions for improvement for this calculation would be the inclusion of the NLO QCD corrections to the gluon fusion channel and of the aforementioned NLO EW corrections.
The code used for the simulations presented in this work is available upon request from the authors and will be made public in a future release of Geneva. [5] Measurement of the total ZZ production cross section in proton-proton collisions at √ s = 8 in 20 fb −1 with the ATLAS detectorarXiv:ATLAS-CONF-2013-020,ATLAS-COM-CONF-2013-020.