ZZ production at high transverse momenta beyond NLO QCD

We study the production of the four-lepton final state $l^+ l^- l^+ l^-$, predominantly produced by a pair of electroweak Z bosons, ZZ. Using the LoopSim method, we merge NLO QCD results for ZZ and ZZ+jet and obtain approximate NNLO predictions for ZZ production. The exact gluon-fusion loop-squared contribution to the ZZ process is also included. On top of that, we add to our merged sample the gluon-fusion ZZ+jet contributions from the gluon-gluon channel, which is formally of N^3LO and provides approximate results at NLO for the gluon-fusion mechanism. The predictions are obtained with the VBFNLO package and include the leptonic decays of the Z bosons with all off-shell and spin-correlation effects, as well as virtual photon contributions. We compare our predictions with existing results for the total inclusive cross section at NNLO and find a very good agreement. Then, we present results for differential distributions for two experimental setups, one used in searches for anomalous triple gauge boson couplings, the other in Higgs analyses in the four charged-lepton final state channel. We find that the approximate NNLO corrections are large, reaching up to 20% at high transverse momentum of the Z boson or the leading lepton, and are not covered by the NLO scale uncertainties. Distributions of the four-lepton invariant mass are, however, stable with respect to QCD corrections at this order.


Introduction
The production of a pair of electroweak vector bosons constitutes an excellent avenue to test the electroweak sector of the Standard Model (SM) at the LHC. This class of processes provides for example information on the non-abelian structure of the Lagrangian. Of particular relevance is the production of a pair of Z bosons, where the underlying gauge structure of SU (2) L ⊗ U (1) Y predicts that tri-linear couplings are actually absent at tree-level in the SM. Furthermore, there is a contribution from an s-channel Higgs boson resonance, produced in gluon-fusion (GF) via a heavy-quark mediated effective Higgs-gluon coupling. This part allows to perform interference studies and off-shell Higgs width measurements [1,2,3,4,5,6].
By the experiments, this process is measured indirectly via Z bosons decaying to a pair of charged leptons or neutrinos each. In this article, we will focus on the four charged-lepton final state pp → + 1 − 1 In the past years, the ATLAS and CMS collaborations have provided a rich collection of measurements with increasing accuracy of the four lepton signal, and limits on the trilinear anomalous ZZγ and ZZZ vertices have been deduced [7,8]. The next-to-leading order (NLO) QCD corrections to ZZ production were first computed in Refs. [9,10] for on-shell production, and including the leptonic decays and spin correlations in Refs. [11,12]. They turn out to be sizable, of around 50%, and exhibit a relevant phasespace dependence. In differential distributions, much larger corrections up to the order of 10 can appear. This makes the approximation of correcting the LO differential distribution by the total K-factor (the ratio of the NLO over the LO total cross section), an unreliable estimate of the true NLO differential cross section and can severely underestimate its size. The origin of the large magnitude of the corrections is twofold. At NLO, new partonic sub-processes appear, including those with enhanced gluonic parton distribution functions (PDFs). This explains in part the size of the total K-factor. On the other hand, some new topologies appear for the first time only at NLO. Among them is a topology with a soft or collinear boson emission from a quark or anti-quark, which results in an α s α EW ln 2 (p T,j /m V ) enhancement for a number of observables [13,14,15].
The one-loop gluon-induced corrections gg → + 1 − 1 + 2 − 2 + X are currently known only at LO. They were first reported for on-shell production in Refs. [16,17]. Results including the leptonic decays [18,19] are also available and studies in the framework of Higgs measurements have also been carried out (e.g. Refs. [1,5,20,21]). Formally, they contribute to ZZ production only at next-to-next-to-leading order (NNLO) QCD, but due to the large gluon PDFs at the LHC, their numerical impact is larger than this naive counting of coupling constants suggests. Depending on the selected cuts, their contribution ranges from a few percent up to ten percent. This prediction suffers from large scale uncertainties. However, results at NLO QCD for the gluonic contributions are expected to be available soon -the real corrections, ZZj production, and the virtual two-loop corrections are already known [22,23].
The NLO electroweak corrections for on-shell ZZ production have been computed in Refs. [24,25]. They yield only a modest contribution, ranging from the few percent level for integrated cross sections up to 10-20% for high-p T observables.
At NNLO QCD, further new partonic channels and new topologies contribute. Thus, to match with the expected experimental precision, it is mandatory to assess the size of these NNLO corrections, not only at the total cross section level, but also for the differential distributions. ZZj production at NLO QCD provides the mixed real-virtual and the double real O(α 2 s ) contributions to the NNLO results. They were first computed in Ref. [26] for on-shell production and account for the new sub-processes and the new topologies appearing for the first time at NNLO. Thus, they are expected to provide the dominant contribution to the total NNLO prediction for selected observables.
The NLO QCD corrections to the double real-emission process, pp → + 1 − 1 + 2 − 2 jj + X have been reported recently [27] with corrections around 10%. The size of the corrections is relatively mild, if adequate central scales are chosen, due to the absence of new channels and phase-space regions opening up at this order, although the uncertainty from varying the factorization and renormalization scale gets greatly reduced.
The two-loop virtual corrections for off-shell ZZ production have been presented in several publications [28,29,30,31] and the NNLO QCD result for the inclusive total cross section has been reported recently [32]. The size of the NNLO corrections compared to the NLO result is about 10%. Up to date, no fully differential NNLO predictions are available in the literature.
In this article, we employ the LoopSim method [13,33] together with the NLO predictions for ZZ and ZZj production and the LO ones for GF-ZZ and GF-ZZj, calculated by the Monte Carlo program VBFNLO [34,35,36]. From this procedure, we obtain merged samples that, in certain regions of phase space, are expected to account for the dominant part of the NNLO QCD corrections to the ZZ production process. The abovementioned procedure has been used recently for other diboson production processes [14,15], as well as V+jet processes [37], resulting in corrections ranging from 30-100% for selected distributions. Such sizable corrections should be taken into account in experimental analyses.
The article is organized as follows: In Section 2, the details of the theoretical framework of our calculation are given. In section 3, first, we compare our predictions with the existing ones presented in Ref. [32] for on-shell Z pair production at the integrated cross section level at NNLO. Afterwards, differential distributions are also presented for two set of cuts -in Section 3.2.1, the setup of the ATLAS and CMS experimental analyses on ZZ production has been closely followed, while in Section 3.2.2, Higgs search cuts are imposed, following the CMS analysis of Ref. [38]. Finally, in Section 4, we present our summary and conclusions.

Theoretical Framework
Production of the four-lepton final state happens at leading order mainly via the quark-antiquark t/u-channel diagram (see Fig. 1 (a) for a representative Feynman diagram). The bulk of the contribution comes from on-shell ZZ production, as in that case the electroweak coupling in the Z decay gets effectively replaced by the corresponding branching ratio, which is about 10% for a Z boson decaying into a pair of any charged leptons. Instead of Z bosons, also q q g g q Z/γ * Figure 1: Representative Feynman diagrams for the different production processes contributing Top: qq processes appearing at tree-level. Bottom: GF-initiated processes appearing at the one-loop-squared level. diagrams appear where these are replaced by virtual photons, denoted as γ * in the following. These typically yield lepton pairs with invariant masses much smaller than the Z mass. Their overall contribution strongly depends on the lepton cuts imposed on the final state. Typical experimental invariant mass windows for Z bosons have a lower bound of 66 GeV, which reduces the γ * contribution to a negligible level.
Another possibility to produce this four-lepton final state is production of a single vector boson V ∈ (Z, γ * ), which in turn undergoes a four-body decay into the final-state leptons. An example Feynman diagram is depicted again in Fig. 1 (b). This s-channel contribution is also sub-dominant in the SM, since there are no tree-level tri-linear gauge couplings and selection cuts on the final-state leptons suppress these contributions due to the limited phase space to simultaneously produce the two intermediate vector bosons close to their mass shell.
Finally, there are one-loop-squared GF diagrams that can also generate the same four-lepton final-state. These can either proceed via an s-channel Higgs resonance, which subsequently decays into + The other possibility is a continuum production of two, potentially off-shell, Z/γ * bosons through a quark-loop box diagram. For typical inclusive cuts, the bulk of the GF contributions originates from on-shell ZZ production, similar to the t/uchannel diagrams. Production via an s-channel vector boson resonance is forbidden by gauge invariance and the Landau-Yang theorem [39,40].
For simplicity, if not stated otherwise, we will refer to the whole process by ZZ production, although we will consider all off-shell effects, non-resonant diagrams and spin correlations of the four-lepton final state. Also, adding the contributions from virtual photons γ * is implicitly understood if not mentioned otherwise.

ZZ and ZZj production in VBFNLO
Our calculation relies on the following ingredients. The NLO QCD corrections for the ZZ and the ZZj production processes, as well as the LO one-loop gluon-fusion-induced gg → 2 j contributions are obtained from the VBFNLO package. The NLO QCD corrections to ZZj were included first for this study, and, in the following, we give some details of the methodology used to compute them to make this work self-contained. We follow closely the strategies used for ν γγj production [41].
We use the spinor-helicity amplitude method and the effective current approach [42,43] to factorize the electroweak part of the system, containing the leptonic tensor, from the QCD amplitude. We first generate the generic qq → V 1 V 2 j amplitudes, where, here and in the following, V denotes either a, possibly virtual, Z ( * ) boson or a virtual photon γ * . Additionally, we also calculate the contribution qq → V j. Then the leptonic decays V i → + At NLO, we need to compute the virtual and the real corrections. To regularize the ultraviolet (UV) and infrared (IR) divergences, we use dimensional regularization [44] and the anti-commuting prescription of γ 5 [45]. We employ the Catani-Seymour algorithm [46] to explicitly cancel the IR divergences prior to the phase-space integration.
To evaluate the scalar integrals, we follow the prescription of Refs. [47,48,49] and for the one-loop tensor coefficients, we employ the Passarino-Veltman reduction formalism [50] up to the box level, but avoiding the explicit appearance of Gram determinants [43,51], and the formalism of Refs. [52,53], with the notation laid out in Ref. [43] for pentagon integrals.
Several checks have been applied to the virtual amplitudes, computed with the package described in Ref. [43], among them, the factorization of the poles, gauge invariance, and parametrization invariance [43]. For the virtual contributions with closed fermion loops, we set the mass of all the fermions to zero and check that the poles add to zero. Additionally, we have cross-checked both the continuum and Higgs resonance graphs with a second implementation based on FeynArts and FormCalc [54,55,56,57]. On the amplitude level, there are 10-14 digits of agreement for each case. The integrated leading order and real emission contributions have been checked against Sherpa [58] and agreement at the per mille level has been found. Additionally, we have cross-checked that our predictions agree at the per mille level with the ones provided in MCFM [59].
A similar strategy is used for the calculation of ZZ at NLO QCD and the LO gluon fusion We include the leptonic decays via effective currents using the spin-helicity formalism and all off-shell effects including Higgs graphs and photons are taken into account. More details on the treatment of the challenging numerical instabilities appearing in gg → + 1 − 1 + 2 − 2 j can be found in Ref. [22]. For the final-state leptons, we consider Z decays into the first two generations, i.e. electrons and muons. We neglect contributions from Pauli interference due to identical particles in the final-state, which is below the per mille level. Therefore, all four possible decay combinations (V 1 , V 2 decaying into e or µ each) give the same result and need to be produced only once at the generation level. As we will see later, the cuts are different for same-flavor and differentflavor states, and in the Higgs setup case also for electrons and muons, hence after imposing analysis cuts, the contributions of the four combinations will differ. For the Z resonance, we employ a modified version of the complex-mass scheme [60] where the weak mixing angle is kept real. We work in the five-flavor scheme and use the MS renormalization of the strong coupling constant, with the top quark decoupled from the running of α s . Diagrams with a final state top quark pair, which would appear in the real emission part of ZZj, are considered a separate, experimentally distinguishable process and are therefore discarded. Virtual top-loop contributions in contrast are included in our calculation. We consider a massless bottom-quark, m b = 0, except for the closed bottom-quark loop amplitudes, where it is set to its pole mass value. The latter is important to correctly account for the negative interference between topand bottom-quark loops in the effective Higgs-gluon coupling, while for the continuum diagrams the difference between choosing a massless or massive bottom quark is numerically small.

Computing the dominant part of NNLO with LoopSim
For the calculation of approximate NNLO results for pp → + 1 − 1 + 2 − 2 production, we use the LoopSim method [13,33]. This approach allows us to merge, in a consistent manner, the ZZ@NLO and ZZj@NLO samples, provided by VBFNLO, yielding a result which is simultaneously NLO accurate both for the ZZ and ZZj production processes. This merged sample is expected to provide the dominant part of the NNLO QCD correction to ZZ production in phase space regions where additional jet radiation becomes important.
The NLO ZZj predictions provide the double-real and the mixed real-virtual contributions at NNLO of the ZZ result. They are divergent upon integration over the phase space of the real partons and need to be supplemented with double-virtual contributions to yield the finite NNLO result for ZZ production. LoopSim constructs approximate versions of such doublevirtual terms by utilizing the fact that their structure of divergences has to match exactly that of the higher multiplicity contributions. This procedure guarantees finiteness of the combined result, while providing a dominant part of the NNLO result for a number of relevant observables. In the process of determining the exact divergent terms of the two-loop corrections, some finite pieces are generated, which are not guaranteed to match the exact constant part of the twoloop diagram. They are, however, proportional to the LO born kinematics and are therefore negligible for distributions which receive sizable corrections at NLO.
To link the two programs, we have made use of an interface [14] which consists, on one hand, of an extension of the VBFNLO program that writes down events in the Les Houches event (LHE) [61] format at NLO and, on the other hand, of a class in the LoopSim library that reads and processes the LHE events.
The LoopSim method proceeds in the following steps. First, an underlying structure for each NLO ZZj event is determined with the help of the Cambridge/Aachen (C/A) [62,63] algorithm, as implemented in FastJet [64,65], with a certain radius R LS . This step establishes the sequence of emissions in the input event. For simplicity, we combine each pair of oppositelycharged leptons to a virtual vector boson V and LoopSim processes diagrams at this level. Thereby, we use information from the event generation step to identify the leptons connected by a continuous fermion line. For s-channel-type events, which are dominated by contributions where only a single boson is attached to the QCD part of the amplitude, one might consider to also reflect this in the LoopSim part and combine the four leptons into a single particle. The question would then be how to define the transition between the two regions. As this contribution is strongly suppressed by the cuts applied later, we do not pursue this further, but instead always combine the four leptons into two Z bosons. In the next step, the underlying hard structure of the event is determined by working through the ij → k recombinations in order of decreasing hardness, defined by the k t algorithm measure [66,67]. The first n b particles associated with the hardest merging are marked as "Born". The number of Born particles is fixed by the number of outgoing particles in the LO event and is equal to two for the case of ZZ production. At NNLO, the Born particles can be either both vector bosons, a boson and a parton, or two partons. The remaining particles, which are not marked as "Born", are then "looped" by finding all possible ways of recombining them with the emitters. This step generates approximate one-and two-loop diagrams.
In the next step, a double counting between the approximate one-loop events generated by LoopSim and the exact one-loop events coming from the NLO sample with lower multiplicity is removed. This is done by generating the one-loop diagrams from the tree level events first, and then using them to generate all possible one and two-loop events. This set is then subtracted from the full result, which amounts to removing the two-loop diagrams with both loops simulated by LoopSim and leaving only those with one exact and one simulated loop. To satisfy unitarity, these diagrams are assigned a weight equal to the original event times a prefactor (−1) number of loops . This guarantees that the sum of the weights of the LoopSim-generated events is zero [13]. For a fully inclusive observable, the integrated cross section of our approximated NNLO result would be equal to the NLO one (hence the pure α 4 EW α 2 s contributions would vanish). However, in realistic situations, with finite fiducial volumes or in the case of differential distributions, some of the events generated by LoopSim are removed by cuts or reshuffled within histograms, which results in non-vanishing α 4 EW α 2 s correction and that leads to a genuine, approximate NNLO correction.
The jet C/A and k t algorithms mentioned above depend on the radius R LS , which is a parameter of the LoopSim method. The smaller the value of R LS , the more likely the particles are recombined with the beam, the larger R LS , the more likely they are recombined together. The value of R LS is irrelevant for collinear and soft radiation. It affects only the wide angle (or hard) emissions where the mergings between particles i and j compete with mergings with the beam. In our study, we shall use R LS = 1, and we shall vary it by ±0.5. The R LS uncertainty will therefore account for the part of the LoopSim method which is related to attributing the emission sequence and the underlying hard structure of the events.
In order to distinguish our predictions with simulated loops from those with exact loop diagrams, we denote the approximate loops byn, as opposed to N used for the exact ones. With that notation, for processes whose contributions start at tree level like qq → ZZ,nLO denotes the correction with simulated one-loop diagrams, andnNLO is a result with exact one-loop and simulated two-loop contributions. However, for processes that start contributing only at one-loop, like gg → ZZ,nLO denotes the correction with respect to that first, nontrivial result. Hence, it is formally an N 3 LO contribution with respect to the full process of ZZ production.
The GF contribution formally first contributes at NNLO, and consequently we also include it in our mergednNLO sample generated with LoopSim. Due to the large gluonic PDFs, this process can contribute relevantly despite the α 2 s suppression. Hence, and since it is gauge invariant on its own, by now a common approach in the literature is to add this contribution already to the NLO results. We follow this convention, but make the addition explicit by using the label "NLO+LO-GF" in this case. Additionally, as mentioned before, we also merge the real radiation process GF-ZZj computed at LO to the GF-ZZ result, yielding a contribution appearing only at N 3 LO. Our results will in general also include this contribution, where we label the full results as "nNLO +nLO-GF". In summary, the GF contribution is always implicitely understood to be included at the corresponding order given by the power counting of coupling constants, in particular LO-GF in thenNLO result. If additional GF contributions are added, this is made explicit in the label.

Comparison with inclusive NNLO calculation
We start by comparing the results obtained with LoopSim+VBFNLO with the calculation of the inclusive NNLO cross section of Ref. [32]. Here, in contrast to the rest of the paper, we use the settings as those of Ref. [32], i.e. the Z bosons are on-shell and do not decay, hence no cut is placed on the final state. Also, all numerical values of masses and couplings are taken from there, namely m W = 80.399 GeV, m Z = 91.1876 GeV, G F = 1.16639 × 10 −5 GeV −2 , m t = 173.2 GeV, m H = 125 GeV and µ R = µ F = m Z . As PDFs, the MSTW 2008 set [68] is chosen, evaluated at each corresponding order. Since our matrix elements include all off-shell effects and spin correlations in our prediction of + 1 − 1 + 2 − 2 production, for this comparison, we have to modify our code. All the s-channel contributions, like the example diagram in Fig. 1 (b), or the t/u-channel contributions with γ * in the intermediate state are set to zero. Hence, the only contribution comes from t/u-channel diagrams of Fig. 1 with on-shell Z bosons. For the GF part, both Higgs and continuum diagrams contribute, but for the latter we also have to remove all diagrams with virtual photons. In the phase-space generator, the (Breit-Wigner) distributions for the invariant mass of the Z bosons are replaced by δ-distributions at the Z pole mass. As leptonic decays of the Z bosons are still simulated internally, we finally need to account for this by dividing the result by BR(Z → + i − i ) 2 . The resulting cross sections at 8 TeV are shown in Table 1.
At LO and NLO, agreement at the per mille level is found. The NLO integrated K-factor, defined as the ratio of NLO/LO predictions, is 1.46. For comparison, we also show the result of adding the GF contribution, formally NNLO, to the NLO result, evaluating both with NNLO PDFs. This is the best currently available estimate without requiring the evaluation of twoloop diagrams or merging different jet multiplicities. We see that these give an additional 7.7% contribution to the NLO cross section, or a total K factor of 1.57 comparing NLO+GF to the LO result. Additionally, they increase the scale variation uncertainty significantly. The latter comes from the fact that the GF contribution is the lowest order accuracy result for the gg channel.
With respect to NLO, the total NNLO correction, computed in Ref. [32] and quoted in the fourth row of Table 1, is about 12%. Hence, the GF contribution provides the leading part of those, namely 60%. Our approximatenNLO andnNLO +nLO-GF results are shown in the last two rows of Table 1. Their overall agreement with the full NNLO result is good with the difference at the level of 2% only. This is not a priori guaranteed by our method and is consistent with the assumption of a LO × α 2 s effect coming from the genuine finite pieces of the exact twoloop virtual amplitudes. These terms are not properly determined by our method, but they are covered by the remaining scale uncertainty. The scale uncertainty of ournNLO result is similar to NLO+GF, and not reduced like for the NNLO result, since the LoopSim method does not attempt to reconstruct higher-order terms proportional to the scale dependence in order not to underestimate the variation, although, technically, this would be possible.

Differential distributions
In the following sections, results at the LHC at √ s = 8 TeV will be given for two different sets of cuts. In Section 3.2.1, we closely follow the ATLAS and CMS experimental analyses on ZZ production, while in Section 3.2.2, we impose Higgs search cuts, following the CMS analysis of Ref. [38]. Below, we describe the common settings.
The mass of the top and bottom quarks, which run in the closed fermion loops, are set to All other quarks, including external bottom quarks, are taken as massless.
The jets are defined with the anti-k t algorithm [69], as implemented in FastJet [64,70], with the radius R = 0.4. Independently of the order of a prediction, we use the NNLO MSTW2008 [68] PDF set, provided by the LHAPDF [71] implementation with α s (m Z ) = 0.11707.
At fixed order in perturbation theory, the cross section depends on the renormalization and factorization scale. As central values for both of those scales, we choose the scalar sum of the transverse energy of the system where p T,V 1,2 and m V 1,2 are the transverse momenta and invariant masses of the recombined, opposite-signed charged lepton pairs, respectively. The scale uncertainty is obtained by varying simultaneously the factorization and renormalization scale by a factor two around the central scale. Additionally, to assess the uncertainties associated with the recombination method used by LoopSim, we show the uncertainty bands associated with variations of ±0.5 around R LS = 1.

ZZ analysis
In the analysis of ZZ production, we use the following cuts, inspired largely by the ATLAS paper [7]. The settings of the corresponding CMS analysis [8] are comparable. The transverse momenta and pseudorapidities of leptons, as well as those of jets (for observables exclusive with respect to jet activity), and the distance between leptons and leptons and jets are required to stay in the following fiducial volume:  using the cuts of Eq. (5), separated into ZZ and ZZ * event categories defined in Eq. (6). The errors in brackets are the statistical error from Monte Carlo integration, while the percentages give the scale variation error, obtained from varying µ = µ F = µ R ∈ [ 1 2 µ 0 ; 2µ 0 ], with µ 0 given by Eq. 4. For thenNLO results, we additionally give the error due to a variation of R LS between 0.5 and 1.5.
To reconstruct the Z bosons from the leptons, we employ the following algorithm. First, all invariant-mass pairs of same-flavor and opposite-sign lepton pairs are formed. If all leptons are of the same generation, there are in total four possibilities, while for different generations only two exist. The invariant-mass pair closest to the physical Z boson mass is labeled Z 1 and it is required to satisfy the cut 66 GeV < m inv,Z 1 < 116 GeV, otherwise the event is discarded. If the second pair of leptons, which we denote as Z 2 , falls into the same mass window, the event is labelled as ZZ, if not, it is called a ZZ event, provided that m Z 2 > 20 GeV. If the latter is not satisfied, the event is rejected. Hence, our two selection types can be summarized as where m Z,max is the maximal mass that can be obtained for a given energy of the system of the incoming partons. Note that the terminology adopted for our study differs slightly from that of Ref. [7], where ZZ was used for the union of both categories defined in Eq. (6).
In Table 2, we present the inclusive cross section at different levels of accuracy. As one can see, the overall behaviour is similar to what we have already observed for total on-shell ZZ production in the previous subsection. The GF contribution gives a significant correction to the NLO result of about +7.2% for the ZZ case, while for ZZ * it is only +3.9%. In both cases the scale dependence is strongly increased. The additional integratednNLO corrections are modest with 1.7% and 3.8% for the ZZ and ZZ * cases, respectively. The dependence on the LoopSim-Parameter R LS is clearly smaller than the remaining scale variation error.
A more important aspect for our method are, however, differential distributions, where the effects can be much larger. To this, we will turn next. Fig. 2 shows distributions of the effective mass observable, H T , defined as a scalar sum of the transverse momenta of leptons and jets, and missing transverse energy The left panel corresponds to the ZZ and the right panel to the ZZ * types of cuts. In the former case, the K factor is very large, both at NLO and atnNLO. In the case of ZZ * selection, the K factor is visibly smaller. The leading correction to this observable at NLO comes from  configurations, shown in the middle diagram of Fig. 3, with one of the bosons emitted collinearly and the other with large transverse momentum recoiling against a hard jet p T, jet p T,Z [72], which results in a dependence given by A similar enhancement occurs at NNLO, where both Z bosons are allowed to be soft or collinear and the result is dominated by the dijet type configurations shown in Fig. 3 (right). This explains both why the K factors grow with transverse momentum and shows that the rate of this growth depends on the selection of the vector boson mass. In Fig. 4, we split the result of Fig. 2 (right) into the two separate mass regions for the offshell bosons: 20 < m Z 2 < 66 GeV Fig. 4 (left) and m Z 2 > 116 GeV Fig. 4 (right). We see that the NLO and thenNLO K factors are much larger in the former case, as putting a smaller mass in Eq. (8) leads to a stronger logarithmic enhancement. We also see that, at large H T , it is the right plot of Fig. 4 that contributes more to the sum shown in Fig. 2 (right), in terms of absolute values. This comes from the fact that the born cross section is already significantly larger in this case. Therefore, Fig. 4 (right), with larger m Z values, which based on Eq. (8) results in a lower K factor, dominates the total result of Fig. 2 (right). That is why the enhancement is smaller there, as compared to Fig. 2 (left), where only the mass region around the Z peak is included.
Let us now turn to Fig. 5, where we show distributions of the transverse momentum of the leading Z boson and that of the leading lepton for the ZZ selection. We see that, in both cases, thenNLO correction is significant, reaching up to 20% with respect to the NLO result at high p T . This magnitude ofnNLO correction at high p T is similar to the one already found in W Z [14] and W W [15] production. We note that, above 200 GeV, the R LS uncertainty is much smaller than the uncertainty coming from variations of the factorization and renormalization scales. Hence, the error related to assigning emission sequence by LoopSim is negligible. The overall theoretical uncertainty decreases by 20-30% as we go from NLO tonNLO for these distributions. We consider thenNLO predictions shown in Fig. 5 as one of the highlights of our study. We emphasize that the inclusion of QCD corrections of that size should be mandatory in all related diboson analyses at high transverse momenta. Such large corrections are of huge importance and should be accounted for both in computations of SM backgrounds and in searches for anomalous couplings.
We notice that the kink around 900 GeV in the p lead T,Z distribution is due to the ∆R , cut of Eq. (5). For a high-p T boson, the decay products are highly collimated with the Z transverse momentum and they share approximately equal amounts of its p T . Hence, the mass of the dilepton system can be approximated by m 2 1 4 p 2 T,Z ∆R 2 , . The imposed phase-space cut ∆R , > 0.2 then leads to the condition p T,Z 10 m . Since the m , distribution is strongly peaked at the Z boson mass, we obtain that the typical separation between the leptons drops below 0.2 at p T,Z 900 GeV. This leads to the kink observed in the LO distribution of Fig. 5. Also, because the two Z bosons are back-to-back in LO configurations and hence have the same p T , the same effect happens simultaneously for both Z bosons. With additional parton radiation, this effect is smoothed out. The sub-leading Z boson will in general have a smaller p T , as some transverse momentum is carried by the parton. The kink at about 800 GeV in the p max T, distribution, shown in the right panel of Fig. 5, is also due to the ∆R , cut. We have checked explicitly that changing the ∆R , cut to 0.4 moves the kink to 450 and 400 GeV for the p lead T,Z and p max T, distribution, respectively, consistent with the above discussion. Another interesting distribution is the mass of the ZZ system or, equivalently, the mass of the system of four leptons, m 4 , which is shown in Fig. 6. This observable is particularly important from the point of view of anomalous triple gauge coupling (aTGC) searches, as the effects of anomalous couplings are enhanced in events where large momentum is transferred through the triple-boson vertex [7,8]. Such events result in a large invariant mass of the diboson system. It is interesting to find that the m 4 distribution receives only modestnNLO   corrections from QCD, which stay always below 5%. This is to be compared with typical sizes of electro-weak corrections, which are not accounted for, and can be of the order of 10% in the tail of the distributions, and with the PDF uncertainties, estimated at 5-10%. Hence, compared to the above sources of uncertainties, thenNLO QCD corrections are small and we conclude that the m 4 distribution becomes stable at this order and can be safely used for setting aTGC limits [7,8]. Fig. 6 (right) shows the distribution of the transverse momentum of the four-lepton system. At LO, this distribution is just a δ-function at p T,4 = 0, since there is no jet the four-lepton system could recoil against. We do not show this first bin in the figure. The first non-trivial order for this observable is then NLO. As we go tonNLO, the p T,4 distribution receives significant corrections of the order of 60-80% for the range presented in the plot, consistent with the NLO predictions shown for ZZ+jet production in Ref. [26].
Finally, in Fig. 7, we show the distribution of the x jet variable defined as which has been introduced in Ref. [73] in the context of a dynamical jet veto. Large values of this variable correspond to configurations where most of the energy of the final state is carried by jets. For the case of ZZ production, x jet = 0 at LO and the distribution starts to be nontrivial at NLO. Here we now also show the first bin, which almost exclusively consists of the x jet = 0 contributions, and which we have scaled down by a factor 30 to fit the plot range. At LO, including LO-GF, the result is simply the inclusive cross section of the process. At higher orders, this bin corresponds effectively to the contribution after placing a jet veto on the process. As we see in Fig. 7, at NLO, where one jet is allowed, the x jet distribution is peaked around 0.1-0.15, hence, for most events, the jet carries 10-15% fraction of the final state energy. We also see, however, that there is a non-negligible tail reaching out to x jet = 0.5. When we move tonNLO, the yield of x jet events increases over the entire range of the distribution at the expense of the energy carried by the Z bosons. In particular, the tail receives corrections with K factors of the order of 5 around x jet = 0.5. ThenNLO result suggests that the x jet cut, used in dynamical jet veto analyses [73], should be set around x jet = 0.2, somewhat lower than what could be inferred from the NLO distribution. The R LS uncertainty is negligible for this observable and the renormalization and factorization scale uncertainty is comparable at NLO andnNLO.

Higgs analysis
We turn now to the discussion of the four-lepton production in the context of Higgs analyses. As explained in Section 2 and shown in Fig. 1, the same ZZ → 4 signature can be produced with and without the intermediate Higgs boson. The latter constitutes an irreducible background from the point of view of studies of Higgs production in gluon-gluon fusion, therefore, its precise determination is of utmost importance. The analyses optimized for Higgs studies use slightly different set of cuts with respect to those employed for aTGC searches.
Following the CMS analyses of Ref.
We start by showing in Fig. 8 the distributions of transverse momenta of the leading Z and the leading lepton, focusing only on the GF part. The LO result corresponds exactly to the diagrams of Fig. 1 (c) and (d). ThenLO correction is computed with LoopSim using the gg → 4 + j result of Ref. [22] provided by VBFNLO. We see that the correction to p lead T,Z and p lead T, is large practically over the entire range shown in Fig. 8. The K factor reaches up to 50% for central values and the scale variation does not decrease as we go from LO tonLO (bottom panel). It is important to notice that the R LS uncertainty is negligible for those distributions, hence our prediction is very reliable within the uncertainty coming from the factorization and renormalization scale.
Let us now turn to the distribution of the invariant mass of the ZZ pair. As argued in Ref. [1], even though the Higgs width is extremely small in the Standard Model, interference effects between continuum and Higgs-mediated contributions (diagrams (c) and (d) in Fig. 1) lead to enhanced four-lepton mass spectra in the region m 4 > 2m Z . Hence, the off-shell effects cannot be neglected, and, as shown recently [2,3], by comparing the yield at the Higgs mass peak with that off the peak, they can be used for setting bounds on the Higgs decay width.
Our predictions for the four-lepton mass spectra are shown in Fig. 9. Similarly to the ZZ selection discussed in the previous section, also here,nNLO QCD effects bring a very small correction to this observable. Hence, the m 4 distribution is stable at this order. This is important from the point of view of the off-shell effect studies in GF fusion, since a precise determination of the expected theoretical yield has impact on setting the Higgs width limits.

Summary
We have studied the ZZ production process at the LHC beyond NLO in QCD by merging the ZZ and ZZ+jet NLO samples with help of the LoopSim method. We have included the exact, loop-induced GF predictions, which are part of NNLO, as well as the GF ZZ+jet contribution, which is formally of N 3 LO. The NLO samples were obtained with the VBFNLO package. The leptonic decays of the Z bosons, including all off-shell and spin correlation effects, have been fully taken into account.
We have compared the fully inclusive cross section predictions from our framework with the exact results computed in Ref. [32] and found a very good agreement within 2%, covered by the remaining scale uncertainties, despite the fact that the LoopSim method is missing some finite parts originating from the two-loop virtual contributions. Following closely two experimental setups, one used in the SM ZZ production and aTGC searches, and the other used in Higgs analyses, we obtained results for a selection of differential distributions.
For the observables sensitive to QCD radiation, the corrections exceed the errors bars of the NLO predictions combined with the LO-GF results and range from 20%, in the case of the transverse momentum of the Z boson or the leading lepton, up to 100% for the effective mass variable H T . A study of the x jet observable, defined as ratio of the transverse energy sum of the jets over the sum of transverse energies of jets and Z bosons, shows that thenNLO corrections start to increase quickly when this variable exceeds 0.2. Therefore, when this variable is used to impose a dynamical jet veto, a cut should be placed at around this value. For Higgs-type selection, we investigated also the ZZ production via GF atnLO. The size of the corrections are around 50% for the transverse momenta of the leading Z and 20% for those of the leading lepton.
For observables which favor the LO kinematics, like m 4l , the approximated NNLO QCD corrections are small, of the order of 5%, and comparable with the size of the remaining scale or PDF uncertainties.
Modifications to the VBFNLO program used in this article are available on request and will be part of a future release. The LoopSim library, together with the Les Houches Event interface, is publicly available at https://loopsim.hepforge.org.