Loop-induced $ZZ$ production at the LHC: an improved description by matrix-element matching

Loop-induced $ZZ$ production can be enhanced by the large gluon flux at the LHC, and thus should be taken into account in relevant experimental analyses. We present for the first time the results of a fully exclusive simulation based on the matrix elements for loop-induced $ZZ + 0,1,2$-parton processes, merged and matched to parton showers. The new description is studied and validated by comparing it with well-established simulation with jets from parton showers. We find that the matched simulation provides a state-of-the-art description of the final state jets. We also briefly discuss the physics impact on vector boson scattering (VBS) measurements at the LHC, where event yields are found to be smaller by about 40\% in a VBS $ZZjj$ baseline search region, compared to previous simulation. We hence advocate relevant analyses to employ a more accurate jet description for the modeling of the loop-induced process.


I. INTRODUCTION
Pair production of Z bosons is an important background for Higgs boson production or new physics searches at the CERN LHC. The loop-induced gluon-fusion process gg → ZZ [1] contributes formally only at the next-to-next-to-leading order in perturbative quantum chromodynamics (QCD). Nevertheless, it can get enhanced by the large gluon flux at the LHC, and thus should be taken into account in relevant experimental analyses, including, for example, Higgs-boson related measurements using the ZZ decay channel, both in the onshell [2, 3] and off-shell regions [4,5], tests of the standard model through diboson inclusive production [6,7] and vector-boson scattering (VBS) [8,9], as well as searches for new physics in various forms of heavy resonances [10,11].
The next-to-leading order (NLO) QCD calculation to gg → ZZ has been evaluated in [12][13][14][15], and can increase the loop-induced Born-level result by an amount ranging from 50% to 100%, depending on the renormalization and factorization scale choices. The loop induced gg → ZZg was evaluated in [16] and can contribute by more than 10% to the next-to-leading order QCD cross section.
On the other hand, experimental analyses employ the full kinematical properties of the events, thus focusing not just on the ZZ-related quantities but also on additional jets, both in terms of production rate in fiducial regions and of their phase-space variables. It is therefore crucial to get as precise predictions as possible for such exclusive observables, which include the dijet invariant mass m jj and the absolute dijet pseudorapidity separation |∆η jj |. It is possible to get such a full exclusive control at hadron level on the complex event topology at the LHC, while still reaching approximately next-to-leading logarithmic accuracy, with the help of recent sophisticated matching methods between matrix elements and parton showers [17,18].
We present here the results of a fully exclusive simulation of gluon-mediated Z pair production based on the matrix elements for loop-induced ZZ + 0, 1, 2 parton(s), merged and matched to parton showers. We examine and validate this new description by comparing it with established simulations where jets are described from parton showers. We find that the matched simulation provide a state-of-the-art accurate description of the final state jets, and is not in agreement with previous simulation for the jet kinematics. We finally focus on the phase space region with two on-shell Z bosons and discuss the impact on the VBS ZZ measurement. We find a large event yield discrepancy (up to 43%) using the matched ZZ simulation.
The paper is organized as follows. We begin by describing our methodology, describing the steps of the matrix element generation in Sec. II and the matching to parton showers in Sec. III. Then we provide the computational details in Sec. IV, and the validation of our results in Sec. V. We show that the matching procedure provides reliable results at the LHC and that the effects are significant. Finally, we discuss briefly the physics impact in Sec. VI and conclude in the last section.

II. MATRIX-ELEMENT SIMULATION
We simulate the loop-induced ZZ + 0, 1, 2 parton(s) processes with LHC settings, using MadGraph5 aMC@NLO [19] version 2.6.5. The event sample is generated at leading order (LO) using a specialized mode to simulate loop-induced processes [20], using the commands below.
generate g g > z z [noborn=QCD] add process p p > z z j [noborn=QCD] add process p p > z z j j [noborn=QCD] It is known that genuine loop-induced diagrams cannot be automatically selected out of all one-loop diagrams in the MadGraph loop-induced mode. In our case, one-loop diagrams also consist of those serving as one-loop corrections to the tree diagrams, which should be excluded in this simulation, as they do not pertain to the gg initial state. A "diagram filter" is especially designed following a suggestion from MadGraph authors [21] to discard those diagrams, using the criteria below.
• The loop in the diagram should not contain any gluon line, so that all vertex-and box-correction diagrams are discarded.
• The loop in the diagram should be connected to at least one Z, W boson or photon, to avoid diagrams representing gluon self-energy corrections through quark lines, and diagrams mediated by a Higgs boson.

3
Once the filter is applied, only genuine gg → ZZ loop-induced diagrams remain. The Higgs-mediated gg → H → ZZ process is also excluded due to the second condition 1 .
For the generation commands, it is worth noting that we specify "pp" instead of "gg" as the initial state to include initial-state radiation (ISR), where an initial-state quark can transform to a gluon through radiation, which then takes part to the hard process. The use of "pp" initial state brings significantly more Feynman diagrams. Only for the 0-jet process, using "pp" is equivalent to "gg" since it does not introduce extra loop-induced diagrams. 1 The main reason is just computing time saving. Although it is known that the Higgs-mediated ZZ production and its interference with the main ZZ production process plays an important role, the main impact is on the total cross section, but negligible on the jet kinematics which is the main topic of this paper. Besides, the benchmark study in Sec. VI requires Z bosons to be on-shell, further diminishing the Higgs contribution. Therefore, the spin correlations of the outcoming leptons are not simulated, and both Z bosons are generated at their pole mass.

III. MATCHING WITH PARTON SHOWERS
The "Les Houches" events (LHE) produced by MadGraph are interfaced to pythia for parton showering. The MLM matching procedure with the shower-k T scheme [18,22] is applied on the 0, 1, 2-parton sub-processes to avoid overlapping of the jet phase space as modeled from the matrix elements and the parton showers. The method introduces a cutoff scale Q ME min in the matrix-element level to remove events with soft partons, and applies another scale Q jet min in the parton showering step, more specifically onto the n + 1 → n differential jet rates (DJR) of the different n-parton sub-process. Under such selection, the final event sample is a combined subset of events coming from each n-parton sub-process.
The jets in each event thus include both harder jets stemming from the matrix-element calculation and softer ones from the parton showers.
The optimal Q ME min value depends on the specific process, while Q jet min has the default value of max Q ME min + 10 GeV, 1.2 Q ME min . To validate the matching procedure in the loop-induced gg → ZZ + 0, 1, 2 jet(s) process, the effect of varying the matching cutoff parameters Q ME min and Q jet min on several distributions, including the DJR, have been extensively studied. We observe for the first time that, in a loop-induced process, the optimal matrix-element level scale Q ME min is smaller than conventionally expected. Fig. 2 (a) shows the DJR distribution for the first and second jet with Q ME min set to 5 GeV. The smooth transition between curves for the different sub-processes signals a good matching result under such parameters. As a comparison, Fig. 2 (b) shows the less smooth DJR distributions with Q ME min set to 10 GeV, i.e., a commonly used cutoff scale. In each case, Q jet min is set to the default value, namely 15 and 20 GeV. For a cross-validation, Fig. 3 further shows the DJR distribution for the first jet with Q ME min = 5 and Q jet min = 15 GeV, under a similar gg → ZZ + 0, 1-parton process. The DJR distributions again illustrate a good matching result. We also check that, by choosing the optimal cutoff value rather than the conventional one, the matched cross section for both the 0, 1, 2-jet and 0, 1-jet processes agrees better with the inclusive result. Besides, we find that moving from the default MadGraph dynamical scale choice to µ = m ZZ does not influence the goodness of the matching. This enables the application of an NLO/LO k-factor correction, which is evaluated using the latter scale choice [12].

IV. COMPUTING PERFORMANCE
The matrix-element level simulation of a loop-induced process with up to two extra partons is implemented for the first time. Because of the complexity of the loop calculation and the large number of diagrams compared to a tree-level process, the event generation turns out to be very time-consuming. We use the MadGraph "gridpack" mode to produce this sample, as it separates the phase-space integration and event Monte Carlo simulation into two steps, making the complicated process easier to handle. On a local cluster with 184 CPU cores, the MadGraph gridpack production takes 24 hours to generate the matrix element code for a loop-induced ZZ + 0, 1, 2 parton(s) process, which contains 24,066 loop diagrams in total. This is followed by a 62-hour concurrent run for phase-space integration.
The events are thereafter generated with the use of the gridpack. Because of the complex phase-space topology, event generation is even more expensive in time: the raw LHE event production rate is 8 min/event, which, when considering a MLM matching rate of about 8%, reaches a net production rate of 100 min/event.

V. VALIDATION
The MLM matched simulation has an intrinsic advantage in describing the jet phase space. As can be seen from the Feynman diagrams in Fig. 1 (b, d,  We compare those simulations for various generator-level jet observables, where jets are reconstructed from final-state particles, adopting an anti-k T algorithm with a distance parameter of 0.4. Fig. 4 shows the spectrum of the generator-level leading and sub-leading jet transverse momenta (p T,j 1 , p T,j 2 ), the dijet transverse momentum p T,jj , and the dijet invariant mass m jj . We first notice that the MadGraph and mcfm simulation with pure parton-shower jets agree well in shape. It is of interest to see that, starting from the Mad-Graph 1-jet simulation, the jet p T and mass distributions gradually turn softer, which could be considered as effects brought by the matrix-element modeled jets. Meanwhile, comparing the MadGraph 0, 1-jet matched simulation with the 0, 1, 2-jet one, we observe similar behavior in the first jet kinematics, but slight discrepancies in the second jet. This is in agreement with our expectations since the in 0, 1-jet simulation the leading jet is modeled by the matrix-element, similarly to the 0, 1, 2-jet simulation. It turns out that the 0, 1, 2-jet simulation gives the softest jets, while it should be the most realistic description amongst all the methods. The shaded region in Fig. 4 represents the combined uncertainties from renormalization and factorization scale (dominant source) and from the parton distribution function, which do not allow to account for the shape differences. The distributions illustrate the sizable discrepancy in dijet phase space modeling between an 0, 1, 2-jet or 0, 1-jet MLM matched description and a full parton-shower description.
It is important to note that the majority of the gg → ZZ (or more generally, gg → V V where V = Z, W ) loop-induced simulations as implemented in various experimental LHC analyses, have the jets modeled in a full parton-shower approach [3,7,8,10,11] or with some approximation using the similarity of the H → ZZ process [2, 5]. Some analyses [4,6,9] apply an alternate approach, i.e., to use Sherpa for the matrix-element simulation of ZZ with zero or one extra jet and merge them with the Sherpa parton showers, but still has not reached an accuracy of dijet matrix-element simulation. Thus, the discrepancy in Fig. 4 should be considered carefully in the relevant analyses.

VI. PHYSICS IMPACT
The innovation from the earlier, less accurate, modeling of the dijet phase space to the new 0, 1, 2-jet matched description may bring potential impacts in relevant analyses. We take a typical VBS ZZ measurement on LHC at an integrated luminosity of 150 fb −1 as an example and discuss the impact based on generator-level simulation. Referring to the object reconstruction and event selection strategy in jj channel in the ATLAS and CMS study [8,9], we design an algorithm to select four generator-level leptons and reconstruct the ZZ pair. A ZZjj baseline selection [8] is designed to select the signal, imposing a jet requirement on the leading and sub-leading jets, namely p T,j 1,2 > 30 GeV and m jj >  Table I shows the yields with statistical uncertainties for the four simulated samples after the ZZjj baseline and VBS-enrich selection respectively. We see that the event yields after the ZZjj baseline selection decrease by 43% for this 0, 1, 2-jet description and, less significantly, by a factor of 34% for the 0, 1-jet description and 9% for the 1-jet description. The event reduction reaches up to 56% when moving to the tighter  result in Table I. We summarize as follows: • As a consequence of validation results from Sec. V, the softness of jets modeled in the MLM matched simulation may cause lower baseline selection efficiency, since a typical VBS region favors high-p T jets. Hence the yields are generally smaller in each distribution of Fig. 5 and in Table I. The compared 1-jet and 0, 1-jet matched simulations also show a decrease in event yields, but not as significant as for the 0, 1, 2-jet case.
• As shown in the |∆η jj | distribution, the 0, 1, 2-jet simulation with dijet simulated from matrix elements induces larger separation of jets, compared to the 0, 1-jet simulation where one jet is produced from matrix-element and another from parton showers, however the discrepancy in |∆η jj | vanishes after the VBS-enriched selection.
• Since the ZZ pair recoils against the emitted jets, the softness of jets may in turn cause a larger m ZZ . This explains the increasing ratio of yields in higher bins of the m ZZ distribution.
It is evident that the improvement in gg → ZZ simulation may impact the total background estimation, and thereafter influence the fit results of VBS signal searches. We note that a similar behaviour may appear in other analyses with the change of jet description using improved loop-induced simulation.

VII. CONCLUSIONS
In this study, we present for the first time the results of a fully exclusive simulation based on the matrix elements for loop-induced ZZ + 0, 1, 2 parton processes, merged and matched to parton showers. We find the optimal MLM matching cutoff scale to be smaller in this simulation compared to commonly used values. We examine and validate this new description by comparing it with various loop-induced ZZ simulations, including two mcfm and MadGraph ZZ simulations with jets from parton showers, a MadGraph simulation with ZZ + 1-jet simulation, and an analogous MadGraph matched simulation for 0, 1partons. We find that the 0, 1, 2-parton matched simulation provides the most state-ofthe-art exclusive description of the final state jets, despite its high complexity in event generation. Jets modeled from the 0, 1, 2-parton matched description are found to exhibit a generally softer transverse momentum spectrum compared to pure parton-shower jets.
We also briefly discuss the physics impact on vector boson scattering (VBS) ZZ measurements at the LHC. By replacing the earlier loop-induced gg → ZZ simulation with the new 0, 1, 2-jet matched simulation, we observe a decrease of 43% in event yields for a typical VBS ZZjj baseline region, and of 56% for a tighter VBS-enriched region. We also observe significant discrepancies in the generator-level jet or reconstructed Z boson kinematics among the different modeling approaches. We hence suggest the implementation of a more accurate description of the emitted jets in this process.