Four-lepton production in gluon fusion at NLO matched to parton showers

We present a calculation of the next-to-leading order (NLO) QCD corrections to gluon-induced electroweak gauge boson pair production, $gg \to ZZ$ and $gg \to W^+W^-$, matched to the PYTHIA8 parton shower in the POWHEG approach. The calculation consistently incorporates the continuum background, the Higgs-mediated $gg\to H^* \to VV$ process, and their interference. We consider leptonic decay modes of the massive vector bosons and retain offshell and non-resonant contributions. The processes considered are loop-induced at leading order and thus contain two-loop virtual contributions as well as loop-squared real contributions. Parton-shower effects are found to be marginal in inclusive observables and quite sizeable in observables that are exclusive in additional jet radiation. The Monte Carlo generator presented here allows for realistic experimental effects to be incorporated in state-of-the-art precision analyses of diboson production and of the Higgs boson in the offshell regime.


Introduction
One of the main objectives of Run 3 of the Large Hadron Collider (LHC) will be the further investigation of the Higgs sector. Most studies directly targeting the Higgs boson will focus on its onshell production and subsequent decay. Indeed, one might naively expect that the cross section to produce an offshell Higgs boson is negligible, due to the extremely narrow width of the Higgs a e-mail: simone.alioli@unimib.it b e-mail: silvia.ferrarioravasio@physics.ox.ac.uk c e-mail: j.lindert@sussex.ac.uk d e-mail: raoul.rontsch@cern.ch boson of about 4 MeV in the Standard Model (SM). However, contrary to this expectation, it is known that approximately 10% of gg → H * → V V events are produced with an invariant mass m V V above the 2m V production threshold [1]. The importance of the offshell region for Higgs phenomenology was further highlighted in Ref. [2], which showed that a comparison of onshell and offshell data can provide stringent constraints on the width of the Higgs boson (see also Refs. [3,4]). While later work indicated that such constraints are not model-independent, they also revealed the potential of using offshell data to probe the couplings of the Higgs boson [5][6][7][8][9][10][11]. Offshell analyses have been performed by both ATLAS [12,13] and CMS [14][15][16][17], and have succeeded in constraining the Higgs boson width to O(10 MeV). This is several orders of magnitude smaller than a direct constraint, which is limited by the detector resolution. Nevertheless, offshell analyses are currently still limited by the available statistics. Further studies of offshell Higgs boson production will therefore be a key component of the investigations of the Higgs sector during both Run 3 and in the high luminosity phase of the LHC.
In this paper, we will focus on the production of an offshell Higgs boson through gluon fusion and its subsequent decay into a pair of electroweak gauge bosons. To this end we consider the signal Higgs production process gg → H * → V V together with the corresponding continuum background process gg → V V and their interference. We study the two diboson modes V V = {ZZ, W + W − } and we assume leptonic decays of the diboson pair. In the following, for brevity, we often denote the processes according to the intermediate diboson res-onances (ZZ, W + W − ). However by this we always refer to the full four-lepton offshell processes, including the interference between Z and offshell photon production.
The signal process proceeds predominantly through a top-quark loop. For onshell Higgs production, the top-quark mass is the largest scale in the process and can be approximated as infinitely heavy, allowing this loop-induced process to be reduced to a tree-level one. Using this approximation, the next-to-next-to-next-toleading order (N3LO) corrections to Higgs production have been computed [18][19][20]. However, this approximation is not valid for offshell Higgs production, since the virtuality of the Higgs may be comparable to (or even larger than) the top-quark mass. This means that a leading-order (LO) prediction for offshell Higgs production requires the computation of a one-loop amplitude with the full top mass dependence, while the next-to-leading order (NLO) correction requires a twoloop amplitude. By itself, this would not be so onerous, but there is a second reason why predictions for offshell Higgs production are more demanding than for onshell Higgs production. It is well-known that the interference effects between the signal gg → H * → V V and the background process gg → V V can be sizeable and thus must be taken into account [1]. Moreover, as we discussed above, the impact of top quarks in the loops cannot be neglected, and this means that in the computation of the background amplitudes gg → V V the contribution from both massless and massive quarks circulating in the loops should be considered.
Results for offshell Higgs production including the mass dependence of quarks in the loop and interference effects are known at LO [1,3,4,21]. Results in the presence of an additional radiated jet have also been presented [22,23]. At NLO, the two-loop gg → V V amplitudes for massless quarks circulating in the loop have been known for several years [24,25]. However, the corresponding amplitudes for massive quark loops have only recently become available [26][27][28]. This means that a fully consistent NLO prediction with the exact dependence on the top-quark mass for offshell Higgs production is in sight but still not available.
However, NLO calculations including interference effects have been presented based on an expansion in 1/m t [29][30][31]. This expansion is not valid for high energies, but has been shown to work well below the toppair production threshold 2m t . In fact, Ref. [30] uses a conformal mapping and Padé approximants to extend the results beyond the top-pair threshold. More recently, it has been demonstrated that using an expansion in 1/m t together with a threshold expansion as inputs for Padé approximants can lead to improved estimates for both gg → HH and gg → V V amplitudes [32,33]. In Ref. [34] the massive two-loop amplitude for gg → ZZ has been computed in the highenergy expansion s, |t| m 2 t , which opens the door for a NLO description of this process in the phase space m V V > 2m t . However, even disregarding these methods, there is a significant region of the offshell phase space with m V V < 2 m t in which the 1/m t expansion is expected to be reliable, and hence where a good approximation to the NLO corrections can be obtained. We base the Monte Carlo generator presented here on such an approximation, following the calculation of Ref. [31]. In the future, the generator can easily be extended to also cover the region m V V > 2 m t by replacing the massive two-loop amplitudes.
Reliable NLO corrections to the continuum background gg → V V alone can be obtained ignoring heavy quark contributions (or these can be incorporated via a reweighting of the massless two-loop amplitude with the LO mass dependence). They are available in the literature both for gg → ZZ [35,36] and gg → W + W − [37,38]. 1 Formally these are of O(α 3 S ) with respect to the LO pp → V V process, i.e. they contribute beyond the order of the known NNLO corrections to the quarkinduced channels [39][40][41][42][43]] -yet they yield phenomenologically relevant contributions.
The NLO results of Refs. [30,31,[35][36][37][38] are at fixedorder parton level, meaning that they do not account for radiation beyond one additional jet. This, together with the fact that unweighted events are not available, hinders the use of these calculations in experimental analyses. In this paper, we report on NLO calculations for offshell Higgs production, including interference effects, matched to parton showers using the POWHEG method [44][45][46][47]. The implementation extends earlier work by two of us [48] that considered the background process gg → ZZ → 4 only. Furthermore, in contrast to Ref. [48], here we also include the contribution from qg-and qqinitiated channels. This implementation allows the generation of unweighted events with additional radiation included through the parton shower, and should facilitate the use of the NLO calculations in experimental analyses. The corresponding POWHEG-BOX-RES generator gg4l will be made publicly available in due time.
The paper is organized as follows. In Sec. 2, we briefly discuss the technical details involved in the parton-level calculation as well as in the matching procedure. In Sec. 3, we summarize the numerical inputs that we use. In Sec. 4, we present fixed-order results validating our calculation and investigate the applied approximations. Finally in Sec. 5 we present numerical results for ZZ and W W production matched to parton showers. We conclude in Sec. 6.

Computational setup
In this section, we describe the matching of the NLO calculation of gluon-induced four-lepton production to parton showers through the POWHEG method implemented in POWHEG-BOX-RES. We first describe the structure of the fixed-order NLO computation and then discuss several details relevant for the matching to PYTHIA 8.

Structure of the NLO computation
We begin by summarizing the salient features of the NLO calculation, and refer the reader to Ref. [31] for additional discussion. As mentioned in the previous section, we need to consider both Higgs-mediated amplitudes gg → H * → V V as well as continuum production gg → V V amplitudes. We therefore write the full amplitude for gluon-induced V V production as where A signal refers to Higgs-mediated amplitudes, while A bkgd refers to amplitudes without any Higgs propagators. Squaring this equation gives (2) Upon integrating over the phase space for the final state particles, the first two terms on the right-hand side give the signal and background results, respectively, while the third term gives the interference contribution dσ full = dσ signal + dσ bkgd + dσ intf . ( In Secs. 4 and 5 we will present results for these contributions separately, as well as for their sum dσ full . As mentioned in the previous section, the LO amplitudes for both A signal and A bkgd are well known [1,3,4,21,[49][50][51]. At NLO, we have to compute the real and virtual corrections to A signal and A bkgd . The corrections to A signal have been known for some time [52][53][54][55]. On the other hand, the NLO corrections to the background amplitude A bkgd are more involved, and deserve a separate discussion. We begin by examining the virtual corrections to the gg → ZZ process. In this case, one can clearly separate massless loops of the first five flavours, and massive topquark loops. The virtual (two-loop) amplitudes for the former are known [24,25], and we construct these using the ggVVamp library [25]. Results for two-loop amplitudes with massive quarks were presented very recently [26,28]. However, here we follow the approach of Fig. 1 Real corrections to gluon-induced V V production, with different partonic channels.
Refs. [29,31] and use an expansion in 1/m t for the massive amplitudes. This implies that our NLO results for the ZZ production process are only valid below the top pair production threshold m ZZ < 2m t . Finally, we need to include double-triangle amplitudes, where each triangle can have either massless or massive quarks in the loop. We employ analytic results for these amplitudes taken from Refs. [56,57]. We now discuss the case of W W production. Since top and bottom quarks mix in the loop, there is no longer a clear division into massive and massless loops. For this reason, Ref. [31] only considered four massless quark flavours in the loop for W W , neglecting the third generation entirely. Here we take a slightly different approach. We compute the two-loop amplitudes A 2loop bkgd using ggVVamp assuming massless quarks for four active flavours. We then reweight the amplitudes as follows where A 1loop bkgd (u, d, s, c, t, b) is the one-loop amplitude with full mass dependence for the third-generation quarks and A 1loop bkgd (u, d, s, c) is the one-loop amplitude with four active flavours. 2 We will comment on the accuracy of this approach in Sec. 4.2.
The real corrections to gg → V V include both the purely gluonic channel gg → V V +g as well as channels with initial state quarks qg → V V +q and qq → V V +g (see Fig. 1) and their crossings. At O(α 3 s ), the former can be unambiguously identified as corrections to the loop-induced process gg → V V . The qg and qq channels are more intricate. These channels also appear in the O(α 3 s ) corrections to the qq → V V process, and it is not possible to parametrically distinguish these corrections from the corrections to the loop-induced process that we are interested in. For this reason, these channels were not included in Ref. [31]. On the other hand, there is no obstacle to computing these corrections, as they form a gauge invariant subset and their only infrared singularities are removed through the collinear renormalization of the parton distribution functions. Indeed, results for these channels were included in Refs. [36,38]. In this paper, we choose to include these channels in our nominal predictions at NLO and also in the results after matching to the parton shower. At the same time we will also investigate the impact of these channels in Sec. 5.4, so that both their magnitude and their impact on the scale variation can be properly assessed. Note that we define these contributions to include any amplitudes with at least one vector boson attached to a closed fermion loop. In particular, amplitudes such as those represented in Fig. 1(d) (contributing to gg → ZZ) are included. 3 We compute all the real correction loop-squared amplitudes using OpenLoops 2 [58,59] including massless and massive quark contributions in the loop, allowing us to retain the full dependence on the top quark mass (and bottom quark mass where applicable). This is in contrast to the approach of Ref. [31] where the real amplitudes for gg → ZZ + g involving a top-quark loop were computed using an expansion in 1/m t . Finally, we note that our calculation includes single-resonant amplitudes in all partonic channels.
In order to ensure numerical stability across the whole phase-space, including in the IR regions close to the soft or collinear limits of the real radiation, we rely on the OpenLoops stability system, which automatically reevaluates all phase-space points with the two reduction methods implemented in Collier [60]. For unstable points matrix elements are set to zero. We verified that varying the corresponding threshold stability kill2 by a factor of 10 around a central value of 10 −2 leaves all results unchanged (see Ref. [59] for documentation of this threshold parameter).
In order to optimize the treatment of all the colorless resonances, we take advantage of the POWHEG-BOX-RES framework [47], which, despite being specifically designed to handle the subtractions when intermediate colored resonances are present, can at the same time improve the phase-space sampling of any resonance. This is achieved by first manually specifying the resonance histories. 4 The POWHEG-BOX-RES then decom- 3 In our code, the user can choose to switch off amplitudes with one vector boson attached to an external quark line using the ol noexternalvqq flag. The user can also turn off the qg and qq channels altogether with the select real flag. 4 The automatic resonance-finding algorithm in POWHEG-BOX- RES is not yet able to handle resonances in loop-induced contributions.
poses the cross section into contributions associated to a well-defined resonance structure, which are enhanced on that particular cascade chain. Each contribution is separately integrated at this point with a dedicated resonance-aware phase space sampling which makes use of a resonance-aware subtraction procedure. The resonance-aware subtraction makes use of a mapping from a real to the underlying Born configuration that preserves the virtuality of intermediate resonances.
Due to the absence of QCD divergences in the resonances, the resonance-aware subtraction is strictly speaking not necessary for the processes considered here. However, we choose to adopt it because its usage improves the statistical errors for observables directly probing the resonance structure. The last essential feature of the POWHEG-BOX-RES implementation is the ability to generate remnants and regular events even when the corresponding cross section is negative, which was not possible in previous versions of the POWHEG-BOX that were instead assuming them to be positive. Despite usually being squares of matrix elements, in this process remnants and regulars contributions might indeed assume negative values in the calculation the interference terms. Technical details about necessary modifications in POWHEG-BOX-RES to deal with the processes at hand are given in Appendix A.

Matching to PYTHIA 8
We next discuss the matching of the NLO calculation of gg → V V to the PYTHIA 8 parton shower in the framework of POWHEG-BOX-RES.
The resonance structure that we construct at the partonic level is further preserved by the parton shower by specifying the input resonance cascade chain at the Les Houches event level (LHE) and making sure that the shower does not distort it through recoil effects. This is achieved by using the PowhegHooks class in PYTHIA 8. However, the PowhegHooks class needs to know the number of final state particles involved in the LO process once the resonance decays are stripped. Since in the POWHEG-BOX-RES this number is not fixed, we modified the PowhegHooks class accordingly, following the recipe adopted in Ref. [61].
The PYTHIA 8 parton shower implements two recoil schemes for initial-state radiation (ISR): in both cases the recoil is always applied to all final-state particles to absorb the transverse momentum imbalance due to ISR off an initial-initial dipole. In the default scheme [62] the same is also done for initial-final dipoles. There is, however, also the option to use a fully local scheme [63], in which when an initial-state emission takes place from an initial-final dipole, the final-state spectator absorbs the transverse-momentum recoil and the other particles in the event are left unchanged.
The default recoil scheme is the recommended option to handle the s-channel production of colour singlets, while the alternative one was originally designed to handle deep inelastic scattering and vector boson fusion events. Since at LO the process considered in this paper describes the production of a colour singlet, we maintain as our baseline the default recoil scheme of PYTHIA 8. However, since in principle we can also generate events with a hard final state jet, in our numerical results we compare the two recoil prescriptions for exclusive observables that might be sensitive to this choice.
In all our showered predictions we include underlying event simulation and hadronization effects. However, in order to simplify the identification of the leptons, we turn off QED radiation and the decay of unstable hadrons.

Numerical setup
In this section we present the numerical inputs for the results presented in the following sections.
Coupling and mass input parameters are fixed to the following values: where G F denotes the Fermi constant and with the real-valued weak mixing angle In general, the gg4l generator allows for finite bottomquark masses. In this work, we mostly use the N F = 5 flavour scheme and treat the bottom quark as massless m b = 0 GeV. Only in Sec. 4.1, where we validate against the results of Ref. [31], do we use a nonzero bottom-quark mass, and there we choose m b = 4.5 GeV.
We use the partonic luminosities and strong coupling from the NNPDF30 lo as 0130 and the NNPDF30 nlo as 0118 sets [64] for the validation against the results of Ref. [31] that we present in Sec. 4.1. For all other results, we use the NNPDF31 nlo as 0118 set [65].
We consider center-of-mass energies of 13 TeV, and set as renormalization and factorization scales for all modes where We obtain scale uncertainty bands by independently varying the renormalization and factorization scales by a factor of two and omitting antipodal variations. At the generator level the following kinematic cuts are applied in the ZZ channel, 70 GeV < m 4 < 340 GeV .
We need to impose such an upper cut on m 4 because, as discussed in the previous sections, the virtual corrections are computed using a 1/m t expansion which is no longer valid for large values of m 4 [31]. For W W production we only require to ensure the renormalization and factorization scales remain inside the perturbative domain. We do not impose any transverse momentum or rapidity requirements on the final-state leptons.
In order to avoid numerical instabilities of the loopinduced amplitudes we need to impose additional mild technical cuts at the generation level. For the Born kinematics, we discard configurations where the transverse momentum of the vector boson is smaller than 100 MeV. For the real corrections, we neglect configurations where the transverse momentum of the radiated parton is smaller than 100 MeV, as also done in Ref. [48]. Indeed this region only gives rise to powersuppressed contributions that do not significantly change the total cross-section. We verified that our results are independent of these technical cuts varying them by a factor of 5 from 0.1 GeV to 0.5 GeV.
Finally, we reconstruct jets with the anti-k T algorithm [66] as implemented in the Fastjet package [67,68], with jet radius R = 0.4 and p T,j > 20 GeV.

Fixed-order NLO results
In the following we present selected fixed-order results that we used to validate our implementation and to investigate the accuracy of the applied approximation for the treatment of mass effects in the virtual corrections. ZZ: gg → e + e − µ + µ −

Validation
As a validation of our implementation, we compare the fixed-order LO and NLO cross sections for gg → ZZ → e + e − µ + µ − and gg → W + W − → e + ν e µ −ν µ against the results of Ref. [31] in Tab. 1, with selection cuts as specified in Ref. [31]. The signal, background and interference contributions are shown separately. Following the approach of Ref. [31], a finite bottom-quark mass m b is used everywhere in the signal (|A signal | 2 ). For the ZZ channel, the background (|A bkgd | 2 ) is computed with m b = 0 and the interference (2Re(A signal A * bkgd )) is computed with a finite m b besides from the virtual contribution which is computed analytically in a mixedmass scheme, where the bottom mass is neglected in the background amplitude A bkgd . 5 For the W W channel, A bkgd is evaluated with n f = 4, for both the background and the interference channels.
For this validation we do not consider quark-initiated channels in the real contribution, consistent with Ref. [31]. Within numerical accuracy we find convincing agreement between the two calculations at LO and NLO.

Mass effects
As discussed in Sec. 2.1, the massive contributions to the two-loop virtual amplitudes are incorporated via approximations in our calculation. All other ingredients including the real amplitudes retain full mass dependence. For the ZZ process the approximation for the 5 Conversely to the calculation in Ref. [31], we cannot easily use a different bottom mass value for A bkgd and A signal when evaluating the interference using the Born and real matrix elements provided by OpenLoops.   massive two-loop amplitudes is based on an expansion in 1/m t . As discussed in detail in Ref. [31] the resulting accuracy is estimated to be at the percent level for m ZZ < 2m t and quickly deteriorates beyond this. For the W W process we reweight the massless two-loop amplitudes as detailed in Eq. 4. In order to gauge the impact of this reweighting and thus the accuracy of the applied approximation, we compare against results obtained using only two massless generations for twoloop A bkgd amplitudes, while all other contributions are computed using full mass dependence as usual. We plot these results for the transverse mass m T of the four lepton system in W + W − production in Fig. 2. This observable is defined as where the missing transverse energy E T,miss is given by the neutrino momenta at truth level and φ miss, + − is the angle between the sum of the neutrino momenta and the sum of the lepton momenta. For the background contribution the effect of the reweighting is at the few percent level for the bulk of the cross section and increases up to about 15%-20% at large transverse masses. For the interference the impact is at the 15% level inclusively and mildly increases in the tail of the transverse mass distribution. The sum of all contributions -which also includes the signal where no approximations are needed -receives inclusive variations due to the reweighting procedure of 7%. In the tail this increases to 10%-15%.

NLO results matched to parton showers
In this section we present our numerical results matched to the PYTHIA 8 parton shower. We consider the differentflavour decay modes gg → e + e − µ + µ − and gg → e + ν e µ −ν µ and for simplicity denote them ZZ and W + W − production respectively. The same-flavour leptonic decay modes will also be made available in the gg4l generator, but are not the focus of this study.

ZZ production
In Figs. 3-6 we present numerical results at NLO, LHE level and NLO matched to PYTHIA 8 (NLO+PS) for gluon-induced ZZ production, showing the full result as well as the signal, background and interference contributions separately. In Fig. 3 the invariant mass of the four-lepton system is shown. The Higgs-mediated signal shows the resonance peak at the Higgs boson mass together with the well-known significant offshell tail. The background clearly exhibits a single-resonant peak at m 4 = m Z 6 and increases significantly for m 4 > 2m Z , where both intermediate Z bosons can become onshell. In this region the interference also starts to become relevant. As a consequence of the very inclusive phase-space cuts employed in our numerical analysis, both the signal and the interference reach about 10% of the full result at large m 4 ≈ 2m t , with the interference being destructive. It is well known that the interference provides an even larger destructive contribution at higher values of m 4 , which are however beyond the validity of the 1/mt expansion used in our calculation. The m 4 observable is inclusive in QCD radiation and consequently partonshower corrections are marginal for all contributions (individually and in their sum). In fact, for all production modes the fixed-order NLO prediction agrees at the percent level with both the LHE level prediction and the fully showered prediction. Scale uncertainties at the fully showered level are approximately 20%. At small invariant masses (m 4 < 150 GeV) the interference becomes very small and consequently Monte Carlo statistics deteriorate quickly in this regime. Fig. 4 shows the distribution in where the sum over the transverse momenta considers all leptons and reconstructed jets. In this distribution the signal peaks at H T = m H , while the background peaks at H T = 2m Z . For small H T parton-shower corrections are mostly driven by the first radiation already present at the LHE level. For the background contribution, these corrections are small, but for the signal contribution they lead to a negative correction of about 50%. A possible explanation is that the signal distribution is strongly peaked around m H and therefore very sensitive to additional radiation that moves events away from the peak. For large H T , the parton showers provide substantial positive corrections up to a factor of 2, while the scale uncertainties can be as large as 50%. This effect can be understood as follows. The upper cut on the invariant mass of the four leptons Eq. 8 also restricts H T < 340 GeV at LO. However, the phase space for H T > 340 GeV can be filled via additional QCD radiation. This leads to significant NLO corrections (not shown here), as well as to sizable parton-shower corrections and LO-like scale uncertainties.
Figs. 5 and 6 display the transverse momentum of the four-lepton system and of the hardest jet respectively. For the latter no lower cut on the jet transversemomentum is applied. The two distributions are identical at fixed-order (they only differ in the first bin which for p T,4 includes the Born and virtual contributions proportional to δ(p T,4 )). The fully showered predictions include a Sudakov suppression which can clearly be seen at the lower end of both the p T,4 and the p T,j1 distributions. We also observe that the parton shower changes the sign of the lowest bin in the p T,4 spectrum. This can be understood as follows: the virtual contribution, proportional to δ(p T,4 ), always comes with an opposite sign of the corresponding real contribution.After the shower (and even after the first POWHEG emission) the virtual contribution gets spread out at finite values of p T,4 . This results in a change of sign in the first bin. Turning now to the opposite end of the spectrum, the p T,4 distribution corresponds to the entire QCD recoil of the four-lepton system and for all contributions receives large parton shower corrections in the tail, while LHE level corrections are largest at p T,4 ≈ 40 − 50 GeV and become small in the tail, where the Sudakov suppression fades away. As already discussed in Ref. [48] the large parton-shower corrections can be explained by the fact that, by adding further radiation, the shower increases the transverse momentum of the colour-neutral four lepton system, which has to recoil against the sum of all emitted particles. On the contrary, in the tail of p T,j1 no such enhancement of the corrections due to the parton shower is observed. In fact, by construction the shower emissions are subdominant with respect to the leading jet and on average are separated enough not to be clustered with it. With respect to the LHE level we observe small and negative parton-shower corrections, being compatible within scale uncertainties.

W + W − production
In Figs. 7-10 we present numerical results at NLO, LHE level and NLO matched to PYTHIA 8 for gluon-induced W + W − production, showing again the signal, background, interference, and full results. In contrast to the corresponding results for ZZ production, here we consider the distribution in the transverse mass m T of the four-lepton system, as defined in Eq. 10, instead of the invariant mass of the coloursinglet system. This is shown in Fig. 7. As for the invariant mass in ZZ production, the impact of the parton-shower corrections on the transverse mass in W + W − production is marginal, as expected from its inclusive (with respect to QCD radiation) nature. It is noteworthy that the interference becomes very large at high m T and eventually contributes beyond −50% for m T > 300 GeV. However, also for the interference alone parton-shower corrections are marginal for the entire m T range considered.
A similarly strong enhancement of the interference can also be observed at large H T , as shown in Fig. 8. In the tail of this observable, parton-shower corrections are again sizable. However in contrast to ZZ production, here no upper boundary on the four-lepton invariant mass is applied and the parton-shower corrections level off for large H T at around 50% for the background and 70% for the full. We finally consider the QCD recoil for W + W − production in Fig. 9 and the transverse-momentum distribution of the hardest jet in Fig. 10. We observe similar behaviour as for ZZ production: the anticipated Sudakov suppression at the low end of both spectra, very large parton-shower corrections in the tail of p T,2 2ν , and mild corrections in the entire p T,j1 spectrum.

Shower recoil scheme
As discussed in Sec. 2.2 PYTHIA 8 implements two alternative shower recoil schemes: the default scheme in which the transverse momentum imbalance after an initial-final dipole emission is democratically distributed among all final-state particles, including the four lepton system, and a fully local scheme, in which the recoil is entirely absorbed by the coloured spectator. 7 In Fig. 11 we compare these two schemes considering the transverse momentum of the four lepton system in the background contribution to ZZ production. As already anticipated in Sec. 5.1 the default recoil scheme leads to a very hard spectrum in the tail (with a 50% increase with respect to the LHE distribution around 100 GeV). Conversely the dipole scheme remains close to the LHE level at large p T,4 . For small values of p T,4 , the dominant contribution should arise from several (soft) emissions whose total transverse momentum sum up to zero. However, in the dipole scheme, the transverse momentum recoil for ISR is not always absorbed by the final-state colour singlet. This explains why for very small values of p T,4 the local recoil leads to a significantly smaller cross section compared to the default scheme. Thus, the default scheme yields a better description of the logarithmically enhanced region, while it also overpopulates the hard region of the spectrum.
A detailed discussion of the logarithmic accuracy of the parton shower goes beyond the purposes of this article and the choice of the recoil scheme has important implications at higher logarithmic orders [69][70][71][72][73][74]. However, since the choice of the recoil scheme only affects our predictions beyond the claimed accuracy, a comparison of the two options available in PYTHIA 8 should help assess the size of the total theoretical uncertainty. Differential distribution in transverse momentum of the four charged leptons in gg → e + e − µ + µ − for the background contribution at the LHE level (black), and at the particle level, using the default PYTHIA 8 recoil (blue) and the fully local dipole recoil (violet). In the lower panel the ratio with respect to the LHE level is shown.

Effect of qg and qq channels
In contrast to the calculations in Ref. [31,48], in this study we do include the qg and qq induced channels contributing to the real radiation at NLO. 8 Here we would like to explicitly highlight the impact of these production channels. To this end in Figs. 12-13 we illustrate at the LHE level the impact of the qg and qq channels with respect to only the gg channels for the different production modes, considering the m 4 and H T distributions in ZZ production. We find very similar results also for the W + W − production mode. In the m 4 distribution the impact of the qg/qq channels is rather flat and about 25% for all production modes. For H T it is increasing with increasing H T and reaches up to 50% in the considered range. Clearly, any precision analysis of gg-induced four-lepton production should include these additional partonic channels opening up at NLO.

Conclusions and Outlook
Gluon-induced four-lepton production offers a unique laboratory for the measurements of offshell Higgs bosons. At the same time precision studies of diboson processes 8 These channels were also considered in the fixed-order NLO study of [36,38].  Fig. 12 Differential distribution in invariant mass m 4 of the four-lepton system in gg → e + e − µ + µ − at NLO and LHE level. Colour coding of the different production modes as in Fig. 12. The lower ratio plots show the full LHE contribution including all partonic channels (gg + qg + qq) over only the gg channel contribution for the different production modes.
and corresponding background estimates in new physics searches are becoming sensitive to the accuracy of the modelling of the gluon-induced production modes. Having this in mind, in this paper we presented an implementation of the loop-induced processes gg → ZZ and gg → W + W − including offshell leptonic decays and non-resonant contributions at NLO matched to the PYTHIA 8 parton shower event generator. We consistently include the continuum background contribution, the Higgs-mediated signal, and their interference. All of these are loop-induced processes and therefore their implementation in a fully-exclusive NLO event generator matched to parton showers poses a significant technical challenge.
In inclusive observables, such as the four-lepton invariant-mass distribution in ZZ production, the parton-shower corrections are found to be marginal, while in more exclusive observables like the recoil of the four-lepton system they can become substantial. For the latter we highlighted the importance of the parton-shower recoil scheme. Furthermore we investigated the relevance of the qg/qq induced production  Fig. 13 Differential distribution in H T in gg → e + e − µ + µ − at NLO and LHE level. Predictions, colour coding and bands as in Fig. 12.
channels, which partly overlap with the higher-order corrections to quark-induced diboson production. In our calculation all ingredients have been treated exact at the NLO level apart from the massive amplitudes contributing to the two-loop virtuals, which are incorporated via approximations. Exact results for the latter have become available very recently and could be incorporated in an updated version of the gg4l generator presented here. Moreover, the generator will be made publicly available in the POWHEG-BOX-RES framework.

Acknowledgments
We are grateful to Fabrizio Caola for useful discussions during the preliminary stages of this work. We also thank Paolo Nason for discussing the modifications to POWHEG-BOX-RES necessary for implementing the processes described in this paper. J.L. is supported by the In this section we outline the modifications we implemented to the POWHEG-BOX-RES framework in order to be able to deal with a loop-induced process and with non-positive defined LO processes. These modifications are available in the gg4l process folder and will be incorporated in a future release of the POWHEG-BOX-RES.
As already discussed in Sec. 2, we discard configurations where the one-loop real amplitude becomes unstable, on the ground that this happens only when the p T of the radiated parton is very small and thus, once we include the subtraction terms, we are left only with negligible power corrections. To do so, setting to zero the amplitude computed in setreal is not sufficient and we have modified the subroutine btildereal to ensure that the subtraction terms are included only when a non-zero real amplitude is found.
For the real corrections, we have the possibility to select only the gg → V V g channel. By doing so, we are excluding the collinear divergent term q(q)g → V V q(q). Thus, we have modified the subroutine btildecoll to include the integrated q(q) → gq(q) collinear remnant only when the real q(q)g → V V q(q) contribution is considered.
The NNPDF30 nlo as 0118 PDF set includes nonperturbative corrections which becomes sizeable for scales smaller than m b = 4.5 GeV. This may cause problems in estimating the upper bound for the strong coupling when generating the radiation according to the method detailed in Ref [46]. In essence, the upper bound is computed by using the LO running of α s and a suitable choice of the infrared cutoff Λ rad , which controls the magnitude of the running, so that with b 0 = 33−2×5 12π , and "cmw" denoting the Catani-Marchesini-Webber prescription for the running coupling [75]. We have modified the appropriate subroutine (init rad lambda) in such a way one spans over scales smaller than the bottom threshold to find the appropriate value of Λ rad .
When performing event generation, the subroutines in POWHEG-BOX-RES implicitly assume that both the Born and the real squared amplitudes are positive. This is not the case when we consider the interference contribution alone, which can also be negative. Thus, we have modified the subroutines gen rad isr, pick random and do maxrat to work with absolute values. Furthermore, away from the singular limits, Born and real amplitudes can have opposite signs. When this happens, we always assume that the real contribution is nonsingular and we do not apply any POWHEG Sudakov suppression to it. Therefore we move these nonsingular contributions into the remnants by means of a modified bornzerodamp subroutine.