Next-to-leading-order electroweak corrections to the production of four charged leptons at the LHC

We present a state-of-the-art calculation of the next-to-leading-order electroweak corrections to ZZ production, including the leptonic decays of the Z bosons into μ+μ−e+e− or μ+μ−μ+μ− final states. We use complete leading-order and next-to-leading-order matrix elements for four-lepton production, including contributions of virtual photons and all off-shell effects of Z bosons, where the finite Z-boson width is taken into account using the complex-mass scheme. The matrix elements are implemented into Monte Carlo programs allowing for the evaluation of arbitrary differential distributions. We present integrated and differential cross sections for the LHC at 13 TeV both for an inclusive setup where only lepton identification cuts are applied, and for a setup motivated by Higgs-boson analyses in the four-lepton decay channel. The electroweak corrections are divided into photonic and purely weak contributions. The former show the well-known pronounced tails near kinematical thresholds and resonances; the latter are generically at the level of ∼ −5% and reach several −10% in the high-energy tails of distributions. Comparing the results for μ+μ−e+e− and μ+μ−μ+μ− final states, we find significant differences mainly in distributions that are sensitive to the μ+μ− pairing in the μ+μ−μ+μ− final state. Differences between μ+μ−e+e− and μ+μ−μ+μ− channels due to interferences of equal-flavour leptons in the final state can reach up to 10% in off-shell-sensitive regions. Contributions induced by incoming photons, i.e. photon-photon and quark-photon channels, are included, but turn out to be phenomenologically unimportant.


Introduction
The physics programme of the LHC at Run I was particularly successful in the investigation of electroweak (EW) interactions and culminated in the discovery of a Higgs boson, but no evidence for physics beyond the Standard Model (SM) was found.While the community is looking forward to a major discovery at Run II, an important task is the precise measurement of the properties of the Higgs boson and the other particles of the SM.Small deviations from the predictions of the SM in the observed event rates or distributions might reveal signs of new physics.
One class of processes particularly relevant for tests of the EW sector of the SM is EW gaugeboson pair production.These reactions allow to measure the triple gauge-boson couplings and to study the EW gauge bosons in more detail.Moreover, they constitute a background to Higgs-boson production with subsequent decay into weak gauge-boson pairs and to searches for new physics such as heavy vector bosons.In the Higgs-signal region below the WW and ZZ production thresholds, off-shell effects of the W and Z bosons are of particular importance.In this paper we focus on the production of Z-boson pairs with subsequent decays to four charged leptons, covering all off-shell domains in phase space.While this channel has the smallest cross section among the vector-boson pair production processes, it is the cleanest, as it leads to final states with four charged leptons that can be well studied experimentally.
At Run I both ATLAS and CMS measured the cross section of Z-boson pair production [1][2][3][4] using final states with four charged leptons or two charged leptons and two neutrinos.The results of these measurements are in agreement with the predictions of the SM and permitted to derive improved limits on triple gauge-boson couplings between neutral gauge bosons [5][6][7].Run II allows to improve these measurements, and first analyses have already been published [8,9].
Precise measurements call for precise predictions.The next-to-leading order (NLO) QCD corrections to Z-boson pair production were calculated a long time ago for stable Z bosons [10,11] and including leptonic decays in the narrow-width approximation [12].Once the oneloop helicity amplitudes were available [13], complete calculations including spin correlations and off-shell effects became possible [14,15].Gluon-induced one-loop contributions were evaluated for stable Z bosons [16,17], including off-shell effects [18,19], and studied as a background to Higgs-boson searches [20].NLO QCD corrections were matched to parton showers in various frameworks with [21] and without [22][23][24][25] including leptonic decays.In Ref. [26], a comprehensive NLO-QCD-based prediction was presented for off-shell weak diboson production as a background to Higgs production.Recently, the next-to-next-to-leading order (NNLO) QCD corrections to Z-pair production have been calculated for the total cross section [27] and including leptonic decays [28].Although formally being beyond NNLO in the pp cross section, even the NLO corrections to the loop-induced gluon-fusion channel were calculated [29][30][31] because of their particular relevance in Higgs-boson analyses.
Besides QCD corrections also EW NLO corrections are important for precise predictions of vector-boson pair production at the LHC.EW corrections typically increase with energy owing to the presence of large Sudakov and other subleading EW logarithms [32][33][34][35][36][37] and reach several 10% in the high-energy tails of distributions.In addition, photonic corrections lead to pronounced radiative tails near resonances or kinematical thresholds.Logarithmic EW corrections to gauge-boson pair production at the LHC were studied in Ref. [38] and found to reach 30% for Z-pair production for ZZ invariant masses in the TeV range.Later, the complete NLO EW corrections were calculated for stable vector bosons and all pair production processes including photon-induced contributions [39,40].The size and in particular the non-uniform effect on the shapes of distributions were confirmed.Leptonic vector-boson decays were first included in NLO EW calculations in the form of a consistent expansion about the resonances for W-pair production [41], and in an approximate variant via the Herwig++ Monte Carlo generator for WW, WZ, and ZZ production [42].Most recently, NLO EW calculations based on full 2 → 4 particle amplitudes, including all off-shell effects, have been presented for W-pair [43] and Z-pair production [44] for four-lepton final states of different fermion generations (i.e.without identical particle effects or WW/ZZ interferences).For Z-pair production, the off-shell effects include also the contributions of virtual photons that cannot be separated from the Z-pair signal, but only suppressed by using appropriate invariant-mass cuts.Note that these full off-shell calculations are essential to safely assess the EW corrections below the WW and ZZ thresholds, i.e. in the kinematical region where WW * /ZZ * production appears as background to Higgs-boson analyses.Moreover, a detailed comparison of the full four-lepton calculation [43] to the double-pole approximation for W-boson pairs [41] revealed limitations of the latter approach for transversemomentum distributions of the leptons in the high-energy domain where new-physics signals are searched for.In Ref. [44] we have presented some selected results for the NLO EW corrections to off-shell ZZ production in a scenario relevant for Higgs-boson studies.In this paper we provide more detailed phenomenological studies in various phase-space regions relevant for LHC analyses for pp → µ + µ − e + e − + X and completely new results on pp → µ + µ − µ + µ − + X, including interference effects from identical final-state leptons.We follow the same concepts and strategies as in Refs.[43,44], i.e. finite-width effects of the Z bosons are consistently included using the complex-mass scheme [45][46][47], so that we obtain NLO EW precision everywhere in phase space.We also include photon-induced partonic processes originating from γγ or qγ/qγ initial states.
The paper is organized as follows: Some details on the calculational methods are presented in Sec. 2. Phenomenological results for two different experimental setups are discussed in Sec. 3. Our conclusions are given in Sec. 4.

Details of the calculation 2.1 Partonic channels
The leading-order (LO) cross sections of the two processes pp → µ + µ − e + e − + X and pp → µ + µ − µ + µ − + X receive contributions from the quark-antiquark annihilation channels qq/q q → µ + µ − e + e − , µ + µ − µ + µ − , (2.1) with q ∈ {u, d, c, s, b}.Sample diagrams for these channels, which are generically called qq channels in the following, are shown in Figs.1(a) and 1(b).Note that all LO diagrams involve Z-boson and photon exchange only.There are LO channels with two photons in the initial state as well, γγ → µ + µ − e + e − , µ Z, γ Z, γ Z, γ with corresponding diagrams shown in Fig. 1(c).Due to their small numerical impact, verified in Sec. 3, we consider their contribution as a correction δ γγ to the qq-induced LO cross section and do not include higher-order corrections to these processes.The NLO EW corrections comprise virtual and real contributions of the qq channels, qq/q q → µ + µ − e + e − (+γ), and the real photon-induced contributions with one (anti)quark and one photon in the initial state, generically referred to as qγ channels in the following.

Virtual corrections
The one-loop virtual corrections to the qq channels are computed including the full set of Feynman diagrams.We employ the complex-mass scheme for the proper handling of unstable internal particles [45][46][47].This approach allows the simultaneous treatment of phase-space regions above, near, and below the Z resonances within a single framework, leading to NLO accuracy both in resonant and non-resonant regions.Sample diagrams for the virtual EW corrections are shown in Fig. 2. A first set of diagrams is obtained by exchanging Z bosons or photons in all possible ways between the fermion lines of the tree-level diagrams in Fig. 1: Diagram types (a) and (b) of Fig. 2 would also be present in narrow-width or pole approximations for the Z bosons and contain separate corrections to the production and the decay of the Z boson.Diagrams (c) and (d) feature correlations between the initial and final states or between different Z-boson decays and are only present in a full off-shell calculation.The sample diagrams (e)-(i) cannot be obtained by naive vector-boson insertions between fermion lines.They involve, for example, closed fermion loops (e) or the exchange of W or Higgs bosons.
In our calculation, we perform a gauge-invariant decomposition of the full NLO EW corrections ∆σ NLO into a purely weak part ∆σ weak NLO and a photonic part ∆σ phot NLO .The virtual photonic part is defined as the set of all diagrams with at least one photon in the loop coupling to the fermion lines.The weak contribution is then the set of all remaining one-loop diagrams, including also self-energy insertions and vertex corrections induced by closed fermion loops.The contributions to the renormalization constants have to be split accordingly.This splitting is possible, because the LO process with four charged leptons in the final state does not involve charged-current interactions, i.e. there is no W-boson exchange at tree level.The vector-boson insertions between fermion lines exemplarily shown in the diagrams (a)-(d) of Fig. 2 thus exhaust all generic possibilities how a photon appears in a loop propagator and can systematically be used to construct the virtual photonic contribution.Figure 3 shows the decomposition of the eight diagrams represented by Fig. 2(d) into the purely weak part with only Z bosons in the loop coupling to fermion lines (upper row) and the photonic part with one or two photons in the loop (lower row).Note that the criterion for the splitting considers only the vector bosons in the loop, while it does not refer to the tree-level part of the diagram.The contributions to the field renormalization constants of the fermions are decomposed in an analogous manner.Since only loops with internal photons lead to soft and collinear divergences, the purely weak contribution is infrared (IR) finite.The full and finite photonic corrections to the qq channels, on the other hand, comprise the virtual photonic corrections plus the real photon emission described in the next section.In general, the decomposition of an amplitude into gauge-invariant parts requires great caution.The gauge-invariant isolation of photonic corrections in processes that proceed in LO via neutral-current interactions is only possible because the LO amplitude can be interpreted as a result obtained in a U(1) γ × U(1) Z gauge theory with the same fermion content and the same couplings as the SM.Here, U(1) γ refers to the electromagnetic gauge group with the photon as massless gauge boson and U(1) Z to an Abelian gauge theory with the Z boson as gauge boson.The corrections within this theory, which is a consistent and renormalizable field theory on its own (see e.g.Ref. [48] or Ref. [49], Chapter 12.9), form a gauge-invariant subset of the full SM corrections to our neutral-current process.Since the electric charges of the fermions are free, independent parameters in the U(1) γ × U(1) Z theory, the photonic contributions form a gauge-invariant subset of the corrections (which could even be further decomposed into subsets defined by the global charge factors of the fermions linked by the photon).Note that the subset of closed fermion loops could be isolated as another gauge-invariant part of the EW correction, since each fermion generation (in the absence of generation mixing) delivers a gauge-invariant subset of diagrams.For the process at hand, we, however, prefer to keep the closed fermion loops in the weak corrections.

Real corrections
The real corrections to the qq channels include all possible ways of photon radiation off the initial-state quarks or off the final-state leptons, schematically depicted in Fig. 4. The phasespace integrations over the squared real-emission amplitudes diverge if the radiated photon becomes soft or collinear to one of the external fermions.For IR-safe (i.e.soft-and collinearsafe) observables, however, the collinear final-state singularities and the soft singularities cancel exactly the corresponding soft and collinear divergences from the virtual corrections after integration over the phase space.The collinear initial-state singularities do not fully cancel between real and virtual corrections; the remnants are absorbed into the parton distribution functions (PDFs) via factorization, in complete analogy to the usual procedure applied in QCD.
We employ the dipole subtraction formalism for the numerical integration of the real corrections.In detail, we apply two different variants based on dimensional regularization [50] and mass regularization [51], respectively.The results obtained with the two approaches are in perfect agreement.The underlying idea is to add and subtract auxiliary terms |M sub | 2 at the integrand level which pointwise mimic the universal singularity structure of the squared real matrix elements |M real | 2 and which, on the other hand, can be integrated analytically in a process-independent way.In the dipole subtraction approach, the subtraction function is constructed from "emitter-spectator pairs", where the "emitter" is the particle whose collinear splitting leads to an IR singularity and the "spectator" is the particle balancing momentum and charge conservation in the emission process.Schematically, the NLO correction ∆σ NLO then reads where 2Re [M * Born M virt ] denotes the interference of the Born amplitude M Born with the oneloop amplitude M virt .The last term represents the (IR-divergent) factorization contribution from the PDF redefinition, which takes the form of a convolution over the momentum fraction x quantifying the momentum loss via collinear parton/photon emission.In order to achieve the form of Eq. (2.5), the phase-space measure dΦ of the (n + 1)-particle real phase space, which includes the bremsstrahlung photon, is decomposed into a "reduced n-particle phase space" d Φ without photon emission and the remaining one-particle phase space [dk] of the bremsstrahlung photon according to dΦ = d Φ × [dk].The integration over [dk] involves an integration over the variable x which controls the momentum loss from initial-state radiation.Splitting off the soft singularities developing in this integration of the subtraction terms over x, decomposes the integral 1 [dk]|M sub | 2 into a part proportional to δ(1 − x) and a continuum part.The former contribution can be analytically combined with the virtual corrections, the latter with the factorization contribution, to produce individually IR-finite terms which may be integrated separately in a fully numerical way: Here, |M sub (x)| 2 fin and 2Re [M * Born M virt ] fin are the finite parts resulting from this splitting of the x integration into continuum and endpoint parts.The momenta pi of the reduced n-particle phase space in the first integral are related to the momenta p i of the (n + 1)-particle phase space such that pk → p k + p γ for the emitter k in the collinear limit p k • p γ → 0. The integration over x is a remainder from the factorized one-particle phase space and is only present for radiation off initial-state particles.The Θ(p 1 , . . . ) functions represent the phase-space cuts applied to the particle momenta {p 1 , . . .}, possibly after applying a recombination procedure which is discussed in the next paragraphs.
Quark-induced channels-the collinear-safe setup: Observables that are collinear safe with respect to final-state radiation-as described in Ref. [51]-are constructed by applying an appropriate procedure for recombining radiated photons with nearly collinear final-state leptons.In collinear regions, a photon-lepton pair with photon and lepton momenta p γ and p ℓ is treated as one quasi-particle with momentum p γℓ = p γ + p ℓ if the two-particle separation ∆R γ,ℓ , as defined in Eq. (3.8) below, is beneath a given threshold.Any phase-space cut or any evaluation of an observable is performed for the recombined momenta, while the matrix elements themselves are evaluated with the original kinematics.Both the local and the integrated subtraction terms and the virtual corrections are cut in terms of the n-particle kinematics.Note that the difference |M real | 2 Θ(p 1 , . . ., p n , p γ ) − |M sub | 2 Θ(p 1 , . . ., pn ) is only integrable after photon recombination, because this procedure ensures that the two Θ functions become equal in the soft/collinear regions (up to edge-of-phase-space effects which do not spoil integrability).Since collinear lepton-photon configurations are treated fully inclusively within some collinear cone defined by the photon recombination, the conditions for the KLN theorem are fulfilled, guaranteeing the cancellation of the collinear mass singularity.The formation of such a quasi-particle is close to the experimental concept of "dressed leptons", as e.g.described by the ATLAS collaboration in Ref. [52].
Quark-induced channels-the collinear-unsafe setup: It is not a priori necessary that observables sensitive to photon radiation off a final-state lepton are defined in a collinear-safe way.The reason is that photons and charged leptons may be detected in geometrically separated places, i.e. the photons in the electromagnetic calorimeter and the muons in the muon chamber.This allows the measurement of an arbitrarily collinear photon emission off a muon.In the absence of photon recombination, the lepton masses serve as a physical cutoff for collinear singularities.On the computational side, this simply forbids the recombination of a muon with momentum p µ and a photon with momentum p γ to a quasi-particle of momentum p µγ = p µ + p γ in the collinear regions as one would do in a collinear-safe setup.In the case of photon emission off electrons, the detection of the two particles takes place in the electromagnetic calorimeter.The finite resolution of the detector then defines a natural "cone size" for the recombination of the lepton-photon pair to a single quasi-particle.In our collinear-unsafe setup, we exclude the muons from recombination, while the electrons are recombined with photons like in the collinear-safe case.
In Ref. [53], the dipole subtraction formalism was extended to collinear-unsafe observables.As in the collinear-safe case, the starting point of the formalism is Eq.(2.5), the fundamental difference being that without recombination of a lepton-photon pair, some observables may now be sensitive to the individual lepton and photon momenta within the collinear region.While this is obvious for the real-emission matrix element, this is also required from the subtraction terms in order to guarantee the local subtraction of the singularities.To this end, the reduced n-point kinematics of the local subtraction terms (which are integrated over the (n + 1)-particle phase space) is a posteriori extended to an effective (n+1)-particle configuration with a resolved collinear lepton-photon pair with momenta p1 , . . ., p i = z pi , . . ., pn , p γ = (1 − z)p i . ( Here pi denotes the momentum of a particular final-state emitter in the reduced kinematics, and p i its momentum after collinear γ emission.The energy fraction of the emitter in the leptonphoton pair is denoted by z; it is constructed from kinematical invariants of the (n + 1)-particle phase space.The local subtraction terms are evaluated in the same way as in the collinear-safe case.However, any collinear-unsafe contribution is now cut with respect to the unrecombined (n + 1)-particle phase space: The schematic shorthand notation applies separately to every dipole with final-state emitter momentum p i = z pi after γ emission.Note that the one-particle phase-space integral [dk] in the integrated subtraction terms in Eq. (2.5) is modified, because the z dependence is required in the re-added subtraction contribution in order to allow for cuts on the bare lepton momentum also in the collinear region.This is in contrast to the collinear-safe case where z could be integrated in a process-independent way [51].Splitting off the soft-and collinear-singular contributions in the z-integration, the subtraction terms can be separated into an inclusive part for the collinear-safe case plus extra terms for the collinear-unsafe case.The detailed form of |M sub (x, z)| 2 fin can be found in Ref. [53] and is not repeated here.
For collinear-unsafe observables, the integration over z is not inclusive, so that the conditions for the KLN theorem are not fulfilled.Hence, the collinear singularities from the virtual Figure 5: Photon-induced contributions with the two initial-state splittings γ → q + q * and q → q + γ ⋆ .The star indicates the particle belonging to the initial state of the reduced Born matrix element.corrections do not entirely cancel against those from the real corrections.Using the muon mass m µ as a physical regulator, terms of order α ln(m µ ) remain in the integral and modify the cross section, which often leads to significant shape distortions of differential distributions.Since partonic scattering energies at the LHC are much larger than the muon mass, all terms suppressed by factors of m µ can be safely neglected.From a practical point of view this means that all kinematics is evaluated with exactly massless muons and the relicts from collinearly-sensitive observables remain in the finite but possibly large contributions of order α ln(m µ ).

Quark-photon-induced channels:
The qγ-induced real contributions include all channels with one photon and one (anti)quark in the initial state.Since an external soft quark does not lead to a singularity and since there is no collinear divergence when the final-state quark becomes collinear to one of the final-state leptons, the matrix elements exhibit only collinear initial-state singularities.As illustrated in Fig. 5, they can be grouped into two classes: First, the incoming photon splits into a quark-antiquark pair with the final-state (anti)quark becoming collinear to the incoming photon, or second, the incoming (anti)quark splits into a photon-(anti)quark pair with the final-and initial-state (anti)quark becoming collinear.With the dipole subtraction method each of the two collinear singularities may be locally subtracted with a single dipole whose functional form is given in Ref. [53].The singularities in the two corresponding integrated subtraction terms cancel against the collinear counterterm from the PDFs, which can, e.g. , be found in Ref. [54].

Numerical implementation and independent checks of the calculation
We have performed a complete calculation of all contributions using the publicly available matrix element generator Recola [55] for the evaluation of the virtual corrections and for all tree-level amplitudes at Born level and at the level of the real corrections.The phase-space integration has been carried out with a multi-channel Monte Carlo integrator with an implementation of the dipole-subtraction formalism [50,51,53,56] for collinear-safe and collinear-unsafe observables.The calculation has been cross-checked both at the level of phase-space points and differential cross sections with two other independent implementations.The one-loop matrix elements of the equal-flavour case have been checked against amplitudes from the Mathematica package Pole [57], which employs FeynArts [58,59] and FormCalc [60].The one-loop matrix elements of the mixed-flavour case have been checked against a calculation based on diagrammatic methods like those developed for four-fermion production in electron-positron collisions [46,61], starting from a generation of amplitudes with FeynArts [58,59] and further algebraic processing with in-house Mathematica routines.In all three calculational approaches, the one-loop integrals are evaluated with the tensor-integral library Collier [62] containing two independent implementations of the tensor and scalar integrals.Collier employs the numerical reduction schemes of Refs.[63,64] for one-loop tensor integrals and the explicit results of one-loop scalar integrals of Refs.[65][66][67] for complex masses.The phase-space integration in all three approaches is carried out with independent multi-channel Monte Carlo integrators which are further developments of the ones described in Refs.[68][69][70].For all differential and total cross sections obtained with the different implementations we find agreement within the statistical uncertainty of the Monte Carlo integration.

Input parameters
In the numerical analysis presented below, we consider the LHC at a centre-of-mass (CM) energy of 13 TeV and choose the following input parameters.For the values of the on-shell masses and widths of the gauge bosons we use In the complex-mass scheme, the on-shell masses and widths need to be converted to pole quantities according to the relations [71] The complex weak mixing angle used in the complex-mass scheme is derived from the ratio For details on the complex renormalization of the EW parameters, we refer to Ref. [46].Since the Higgs boson and the top quark do not appear as internal resonances in our calculation, their widths are set equal to zero.For the corresponding masses we choose the values All the charged leptons ℓ = e, µ, τ and the quarks q = u, d, c, s, b are considered as light particles with zero mass throughout the calculation.In the computation of collinear-unsafe observables, the physical muon mass appears as a regulator with numerical value m µ = 105.6583715MeV.(3.4)Note that this non-zero value for the muon mass is only kept in otherwise divergent logarithms from photon radiation off muons, while everywhere else the muons are strictly treated as massless particles.
We work in the G µ scheme where the electromagnetic coupling α is derived from the Fermi constant according to This choice absorbs the effect of the running of α from zero-momentum transfer to the EW scale into the LO cross section and thus avoids mass singularities in the charge renormalization.
Moreover, α Gµ partially accounts for the leading universal renormalization effects originating from the ρ-parameter.The fine-structure constant, is only used as coupling parameter in the relative photonic corrections, i.e. the NLO contribution ∆σ phot NLO to the cross section scales as α 4 Gµ α(0).This choice is motivated by the dominance of corrections from the emission of real photons coupling with α(0).The relative genuine weak corrections, on the other hand, are parametrized with α Gµ , i.e. ∆σ weak NLO scales as α 5 Gµ , since the dominating weak high-energy corrections involve this coupling factor.As we do not consider QCD corrections in this paper, our cross-section predictions do not depend on the renormalization scale µ ren .The dependence of the relative NLO EW corrections on the factorization scale µ fact is only marginal, so that we simply set it equal to the Z-boson mass, µ fact = M Z , without the need to investigate residual scale dependences or alternative scale choices.As PDFs we use the NNPDF23 nlo as 0118 qed set [72]. 1 Throughout our calculation of EW corrections, we employ the deep-inelastic-scattering (DIS) factorization scheme, following the arguments given in Ref. [74].The corresponding finite terms for the EW corrections to be included in the subtraction formalism can be found in Ref. [54].

Definition of observables and acceptance cuts
In the following we define two different event selections: an "inclusive" and a "Higgs-specific" setup.The former uses typical lepton-identification cuts without any further selection criteria; the latter is motivated by specific criteria designed for Higgs-boson analyses by ATLAS [75] and CMS [76].
Inclusive setup: In the collinear-safe case, photons emerging in the real-emission contributions are recombined with the closest charged lepton (cf.Sec.2.3) if their separation in the rapidity-azimuthal-angle plane obeys Here, y j denotes the rapidity of a final state particle j and ∆φ ℓ i γ the azimuthal-angle difference between a charged lepton ℓ i and the photon γ.Note that we only take into account photons with |y γ | < 5 in the recombination procedure, while we consider photons with larger rapidities as lost in the beam pipe.In the collinear-unsafe case, photons are recombined only with electrons/positrons, while no recombination with muons/antimuons is performed.For observables of the equal-flavour-lepton final state, the leading lepton pair (ℓ + 1 , ℓ − 1 ) is defined as the one whose invariant mass is closest to the nominal Z-boson mass; the subleading lepton pair (ℓ + 2 , ℓ − 2 ) is then formed by the remaining two leptons.Leading leptons are thus labelled as ℓ ± 1 , and subleading leptons as ℓ ± 2 .As default setup, we consider a minimal set of selection cuts inspired by the ATLAS analysis [1].For each charged lepton ℓ i , we restrict transverse momentum p T,ℓ i and rapidity y ℓ i according to Any pair of charged leptons (ℓ i , ℓ j ) is required to be well separated in the rapidity-azimuthalangle plane, ∆R ℓ i ,ℓ j > 0.2. (3.10) Higgs-specific setup: For the Higgs-specific setup, motivated by the ATLAS and CMS analyses [75,76], we replace the cuts of Eq. (3.9) by the less restrictive criteria retain the cut (3.10), and complement them with additional invariant-mass cuts on the charged leptons.For the mixed-flavour final state, we require for the two same-flavour lepton pairs with referring to the invariant masses of the lepton pair that is closer to or further away from the nominal Z-boson mass, respectively.For the same-flavour final state, we apply the cuts of Eq. (3.12) after selecting the leading and subleading lepton pairs (ℓ + 1 , ℓ − 1 ) and (ℓ + 2 , ℓ − 2 ) in the same way as described above.The invariant mass M 4ℓ of the four-lepton system is subjected to the cut which is independent of the flavour of the final-state leptons.
In both setups we treat the additional (anti)quark in the final state of the photon-induced contributions in a fully inclusive way, i.e. we do not apply any jet veto.Finally, we note that the employed single-particle lepton identification cuts are chosen to be equal for all charged leptons.Since the lepton pairing in the equal-flavour final state is flavour independent, the results for the e + e − e + e − final state are equal to the results of the µ + µ − µ + µ − final state within the collinear-safe setup.In the following, [4µ] denotes the equal-flavour final state, while [2µ2e] denotes the mixed-flavour final state.
Whenever possible, we have compared the results of our full off-shell calculation with the available results for on-shell Z-boson pair production [39,40].Since our calculation sets phasespace cuts on the charged leptons, a direct comparison is not possible for most observables.Moreover, the calculations for stable Z bosons do not take into account corrections to the Z-boson decays.Finally, there are differences in the setup: The values in the literature are given for a CM energy of √ s = 14 TeV, different factorization scales and PDFs.Nevertheless, the relative EW corrections to those observables that are less sensitive to off-shell effects and corrections to the Z-boson decays can still be directly compared.This holds in particular for the Sudakov enhancement at large transverse momentum of the on-shell Z boson or the corresponding charged final-state lepton pair.

Results on integrated cross sections
The results for the integrated cross sections of the processes pp → µ + µ − e + e − + X and pp → µ For the photonic corrections, we further distinguish between the collinear-safe setup δ phot,safe qq , where the bremsstrahlung photon is recombined with any charged lepton, and the collinear-unsafe case δ phot,unsafe qq , where the muons are excluded from recombination, as described in Sec.We first analyse the inclusive scenario.The major contribution to the corrections stems from the qq channels where the purely weak correction of −4.3% has the highest impact.The photoninduced contribution matters only at the permille level, justifying that we neglect higher-order corrections to this channel.Besides the small photon flux in the proton, yet another reason for the strong suppression is that channels with two photons in the initial state involve at most one resonant Z boson, as illustrated by the sample diagram of Fig. 1(c).We count such kinematical topologies as background topologies to the dominant contribution with two possibly resonant Z-boson propagators as they appear in the qq channels. 2The qγ-induced corrections δ qγ are even further suppressed by yet another order of magnitude and are thus entirely negligible in the integrated cross section.
Summing up all contributions, we find for both the [2µ2e] and the [4µ] final state the same correction of −5.1% in the collinear-safe setup.The corrections to vector-boson pair production with on-shell Z bosons are known to be ∼ −4.5% [39], respectively ∼ −4% [40] in a slightly different setup.The differences can be attributed to the off-shell effects, including also additional virtual photon exchange, and to differences in phase-space cuts and in the employed numerical setup (cf.comments at the end of Sec.3.2).Note that there is no photon-induced contribution of the form γγ → ZZ with on-shell Z bosons in the final state.The strong suppression of the qγ channels confirms similar findings in the on-shell approximation [40].
Comparing the collinear-unsafe photonic corrections with the collinear-safe case, we observe differences in δ phot qq by ∼ 0.7% and ∼ 1.5% for [2µ2e] and [4µ] final states, respectively.In the collinear-unsafe case, final-state radiation off muons is enhanced through the mass singularity ∝ α ln(m µ ) in the phase-space integral.Since the photon is not recombined with the muons, there are systematically more events with a large energy loss in one of the muon momenta induced by final-state radiation.For this reason, less events pass the event-selection in the collinear-unsafe case, leading to more negative corrections.This "acceptance correction" scales with the number of leptons treated in a collinear-unsafe way, explaining the factor of two between the [2µ2e] and the [4µ] cases in δ safe − δ unsafe .
The amplitudes for the equal-flavour final state can be obtained from the different-flavour amplitudes by antisymmetrization with respect to a pair of equal final-state leptons.At the level of squared matrix elements, diagrams with a µ + µ − pair originating from a single vector boson interfere with diagrams where the µ + and the µ − originate from two different vector bosons.These interference terms lead to a deviation from a naive rescaling of the different-flavour cross section by a symmetry factor of two.From the ratio σ LO [2µ2e]/(2σ LO [4µ]) ≈ 1.003 we find a negative interference of about 0.3% for the integrated LO cross section.Comparing this number with the total relative correction of −5.1% for both final states, we conclude that the impact of interferences on the total cross section in the inclusive setup is a "LO effect" in the sense that the relative corrections do not modify this behaviour.The reason for the smallness of the interferences is that the LO cross section is dominated by contributions with two resonant Z-boson propagators.In the interference terms, however, in at least one diagram both Z-boson propagators are off shell.This explains also why the impact of the interferences in the γγ Born cross section alone is less suppressed with σ γγ [2µ2e]/(2σ γγ [4µ]) ≈ 1.12 since this is a background contribution with at most one resonant Z boson.
We now turn to the Higgs-specific setup.Despite the additional cuts of Eqs.(3.12)-(3.13), the cross sections for this scenario are larger than for the previously considered inclusive setup.This feature is due to the less severe cut of 6 GeV imposed on the transverse momenta of the charged leptons in the Higgs-specific setup, as compared to a cut of 15 GeV in the inclusive setup.Moreover, we observe a percent-level deviation at LO from the naive scaling factor of two between the equal-and unequal-flavour cases, σ LO [2µ2e]/(2σ LO [4µ]) ≈ 0.973.This reflects, on the one hand, an expected enhancement of the interference terms since, by construction, the whole scenario is more sensitive to the off-shell effects of the vector bosons.On the other hand, the additional invariant-mass cuts in Eq. (3.12) depend in the equal-flavour case on the chosen lepton-pairing algorithm, while they do not in the mixed-flavour case.A quantitative statement on the size of the pure interference effects would thus only be possible if the same lepton pairing was applied also in the mixed-flavour case.The same arguments apply also to the corrections in Tab. 1, i.e. the differences between the equal-flavour and the mixed-flavour final state are due to interferences and event selection.As a general pattern, we observe that the bulk of corrections to the total cross sections of around −3.5% stems from the weak corrections, while photonic corrections contribute only at the sub-percent level.

Invariant-mass and transverse-momentum distributions
In order to illustrate the impact of EW corrections on differential observables, we present several results on distributions in the inclusive setup at a proton-proton CM energy of √ s = 13 TeV.
We choose the collinear-safe setup as default and provide selected results within the collinearunsafe case subsequently.
Figure 6 shows the four-lepton invariant-mass distribution for the unequal-flavour [2µ2e] and the equal-flavour [4µ] final states.The left-hand side resolves the off-shell region with its threshold and resonance structures with a fine histogram binning, while the panels on the right-hand side show the whole range from the off-shell region over the resonances and thresholds up to the TeV regime in coarse-grained resolution.The absolute predictions of the LO and NLO distributions of both the [2µ2e] and [4µ] final states in the upper panels follow the characteristic pattern of Z-boson pair production: The first peak around M 4ℓ = M Z represents the single-resonant production of the Z boson in the s-channel according to the sample diagram in Fig. 1(b), the threshold at M 4ℓ = 2M Z stems from doubly-resonant diagrams as depicted in the diagram in Fig. 1(a).The knee above M 4ℓ = M Z + 2p T,min ≈ 120 GeV is induced by the kinematical cut of Eq. (3.9) on the lepton transverse momentum.For M 4ℓ M Z + 2p T,min , the cross section is dominated by events with one resonant Z boson (→ ℓ + 1 ℓ − 1 ) and the ℓ + 2 ℓ − 2 pair with M ℓ + 2 ℓ − 2 30 GeV, where both lepton pairs are almost at rest in the transverse plane.Since GeV is necessary for the event to pass the cuts, which drastically reduces the transition rate, since no resonance enhancement is present anymore in the diagram type illustrated in Fig. 1(a).
The panels directly below the absolute predictions for the cross sections show the relative EW corrections to the qq channels in the collinear-safe setup, comparing the purely weak contribution with the full EW contribution (EW=weak+phot).Apart from the off-shell region below the   single Z resonance, we observe that the relative EW corrections of the mixed-flavour final state and the equal-flavour final state are equal over the whole invariant-mass spectrum.This confirms at the level of differential distributions that the interference effect is mainly a LO effect, in accordance with what we have already seen for the integrated cross section.The four-lepton invariant mass in the inclusive setup is well suited to study the relative size of the interferences, as this observable does not depend on the lepton pairing.We show in the lowest panels of Fig. 6 and the following figures the ratio (dσ (N)LO [2µ2e]/dO)/(2dσ (N)LO [4µ]/dO), where O denotes the considered observable, e.g.M 4ℓ in Fig. 6.The LO and NLO curves are, as expected, almost equal.The size of the interference effect varies in the region where no lepton pair is resonant from −7% at M 4ℓ = M Z to +6% at M 4ℓ = M Z + 2p T,min .Thus, the unequal-flavour matrix elements cannot describe the equal-flavour final state there.In the region M Z + 2p T,min M 4ℓ 2M Z , where only one lepton pair can be resonant, the interference effect amounts to 2%.Above the ZZ threshold, the ratio is equal to one up to fractions of a percent, since in this region of phase space the doubly-resonant contribution dominates over any non-resonant interference effect.For higher invariant masses M 4ℓ , the overlap of the two resonance pairs becomes smaller and smaller in phase space, so that the ratio asymptotically approaches one.
We inspect the EW corrections in more detail.In the high-invariant-mass region, the correction is entirely dominated by the purely weak contribution and reaches about −20% around 1 TeV.The M 4ℓ distribution at high scales is not dominated by the Sudakov regime of ZZ production where all Mandelstam variables (s, t, u) of the 2 → 2 particle process would have to be large.Instead, Z-pair production at high energies is dominated by forward/backward-produced Z bosons, where t and u are small.At the ZZ production threshold, the weak corrections change their sign and reach up to 5% below.Note that this non-trivial sign change makes it impossible to approximate the full NLO EW results by a global rescaling factor.At the Z peak, the weak corrections are extremely suppressed.The photonic corrections remain almost constant above the ZZ threshold, at around −2% to −3%, and show the typical radiative tails: Below a threshold or close to a resonance, the LO cross section falls off steeply.Final-state radiation of a real photon, however, may shift the value of the measured invariant mass to smaller values.Since the LO cross section is small in this phase-space region, the relative correction due to the bremsstrahlung photon becomes large.The structure of the radiative tails follows precisely the thresholds and the resonances at leading order: one below 2M Z , one near M Z + 2p T,min ≈ 120 GeV, and another one below the Z resonance.
For completeness, we also show the photon-induced corrections in a separate panel (third row in Fig. 6).Over the whole range of the distribution, both the γγ and qγ contributions are strongly suppressed with respect to the qq processes.Since the γγ channel has only a single Z resonance according to the diagram in Fig. 1(c), it is strongest suppressed with respect to the qq LO cross section above the ZZ threshold and near the s-channel resonance at M Z .In the non-resonant region below M Z + 2p T,min it reaches up to 1%.Since there the LO cross section is small anyway, the overall impact remains small, in agreement with the result for the integrated cross section.Differences between the equal-flavour and the unequal-flavour final states due to interferences are only visible below M Z + 2p T,min where none of the lepton pairs is resonant.The qγ channels contribute over the whole spectrum at most at the permille level.The large correction near the phase-space boundary is phenomenologically irrelevant as the corresponding LO cross section in this region is very small anyway.Due to the negligible impact of any photon-induced corrections in the inclusive setup, we do not show them separately any more in the following plots.
Up to the details in the event selection and the corrections from the Z-boson decays, the four-lepton invariant mass M 4ℓ can be compared to the ZZ-invariant-mass distribution obtained in the NLO calculations of Refs.[39,40] with on-shell Z bosons.The relative corrections to the distribution in the invariant mass of the Z-boson pair at M ZZ = 700 GeV are given as −11% and −8% in Refs.[39] and [40], respectively, while we find for the four-lepton invariant mass M 4ℓ = 700 GeV a relative correction of −15%.We attribute the difference mainly to the final-state radiation off muons missing in the calculation with stable Z bosons.This collinearly enhanced contribution typically leads to negative corrections at the level of some percent, induced by the radiative loss in transverse momenta that can potentially shift events out of acceptance.
The µ + µ − invariant mass M µ + µ − is an example of an observable where the differential cross section in the equal-flavour case directly depends on the employed lepton pairing, as can be seen in Fig. 7.The left column compares the µ + µ − invariant mass of the [2µ2e] final state with the one of the leading µ + µ − pair of the equal-flavour case, while the right column shows the same comparison, with the subleading µ + µ − pair instead.Note that we compare two different observables, since the unequal-flavour case is binned with respect to the flavour but not with respect to the kinematic ordering described in Sec.3.2.Although this precludes us from drawing conclusions on the interference effects (as we did for the four-lepton invariant mass above), one may nevertheless learn from such a comparison to which extent the mixed-flavour case can   be used to describe the equal-flavour case.Obviously, the shape of the subleading lepton pair in the equal-flavour case widely resembles the corresponding observable of the unequal-flavour final state (though it is not equal) while the pattern of the leading lepton pair is fundamentally different.The different behaviour can be explained as follows: The µ + µ − invariant-mass distribution of the mixed-flavour case receives the largest contribution when the corresponding e + e − pair is on the Z resonance.The situation is similar when binning the subleading µ + µ − pair where the corresponding leading lepton pair is always closer to the resonance and, thus, takes over the role of the e + e − pair in the mixed-flavour case for the dominant contribution.Since interference effects are strongly suppressed, this explains the almost constant ratio (dσ[2µ2e]/dM 2ℓ )/(2dσ[4µ]/dM 2ℓ ) ≈ 0.5 above the resonance in the lowest panel in the right column of Fig. 7; it uniformly extends to larger values of M µ + µ − as can be seen in the left columns of Fig. 8.For this particular observable away from the Z resonance, the lepton pairing thus basically identifies the two µ + µ − pairs and the symmetry factor of 1/2 disappears.The large difference near the Z resonance (see lowest panel in the right column of Fig. 7) is due to the fact that the subleading muon pair is always further away from the resonance than the leading pair, in contrast to the mixed-flavour case where the invariant mass of the e + e − pair is independent of the µ + µ − pair.By contrast, when binning the leading lepton pair, the other lepton pair is forced to be further off shell with respect to the Z resonance and, hence, the distribution falls off much steeper with the distance from M Z .The steeper drop in the interval 160 GeV < M µ + µ − < 175 GeV is due to the fact that for invariant masses M lead µ + µ − > 2M Z of the leading µ + µ − pair the invariant mass of the subleading µ + µ − pair cannot be below M Z any more. 3The broad peak between 35 GeV and 60 GeV stems from the single s-channel resonance of the four-lepton invariant mass at Z /6 implying a threshold for the leading lepton pair of M lead µ + µ − > M Z / √ 6 ≈ 37 GeV in agreement with Fig. 7. On the other hand, the Z resonance in M 4ℓ contributes only for M lead µ + µ − < M Z − 2p T,min ≈ 61 GeV, since the transverse momentum of the subleading leptons cannot be lower than p T,min .The s-channel resonance gives also rise to a small bump in the invariant-mass distribution of the subleading lepton pair at somewhat smaller invariant masses (see r.h.s. of Fig. 7).The increase of the distributions towards M µ + µ − → 0 in the subleading and mixed-flavour µ + µ − invariant masses are due to the tail of the photon pole.The spectrum of the leading lepton pair is phenomenologically irrelevant at high invariant masses since it is heavily suppressed due to the lack of a resonance enhancement.
As discussed already in Ref. [44] for the [2µ2e] final state (though for the Higgs-specific setup), the EW corrections largely resemble the known structure of the photonic and weak corrections to Drell-Yan-like single-Z production, which are, e.g., discussed in Ref. [54] (cf.Fig. 12 therein).Let us first analyze the right column of Fig. 7.The weak corrections stay at the 5% level and change sign in the vicinity of the resonance.This can be understood from the fact that in the vicinity of the Z resonance there are two different types of contributions: corrections to the resonant part of the squared amplitude, and corrections to the interference of the resonant and non-resonant parts of the amplitude.The former give a constant offset of ∼ −5%, while the latter are proportional to (M 2 µ + µ − − M 2 Z ) and thus change sign at the Z resonance.This qualitatively explains the observed sign change of the purely weak corrections which is slightly shifted below the resonance due to the negative offset mentioned above.The corrections are to a large extent equal for the mixed and equal flavour case with minor deviations of ∼ 1% in the far off-shell region.Including also the photonic corrections, we observe in both cases the typical radiative tail due to final-state radiation effects below the Z resonance, similar to what has been observed in the four-lepton invariant mass.In the high-energy spectrum of the steeply falling invariant-mass distribution, shown in Fig. 8, both photonic and purely weak corrections are negative.The EW corrections for the invariant mass M lead µ + µ − of the leading µ + µ − pair differ significantly from the mixed-flavour case due to the large differences at LO.At the peak around 45 GeV, the purely weak corrections basically vanish, which is consistent with the four-lepton invariant-mass distribution in Fig. 6 where the purely weak corrections vanish at M 4ℓ = M Z .At the resonance M 2ℓ = M Z , the weak corrections are equal for the [2µ2e] and the [4µ] case because the dominant contribution where the leading and the subleading lepton pairs are both close to the resonance is not sensitive to the lepton pairing (note the ratio of (dσ[2µ2e]/dM 2ℓ )/(2dσ[4µ]/dM 2ℓ ) ≈ 0.5 at the resonance in the lowest panel and the discussion for the subleading case above the resonance).The weak corrections stay always below 5% in absolute size.The photonic corrections exhibit an additional radiative tail below the peak around 45 GeV.The radiative tail below the Z resonance is less pronounced due to the missing resonance enhancement by the subleading lepton pair.
In Fig. 8 (right-hand side) we show the distribution in the transverse momentum of the (subleading) µ + µ − pair, which can be compared with the distribution of the Z-boson transverse momentum in on-shell calculations.Since the two lepton pairs are back to back at LO, the transverse-momentum distribution depends on the choice of the lepton pair only very weakly.The interference effect of a few percent is only visible for small p T,µ + µ − .The EW corrections grow up to −45% for p T,µ + µ − = 800 GeV, while the photonic corrections stay at the level of 1%.Choosing the µ + µ − pair from [2µ2e], makes it possible to compare the p T,µ + µ − distribution of our off-shell calculation with the p T,Z distributions for ZZ production with on-shell Z bosons discussed in Refs.[39,40].However, as mentioned above, it should be kept in mind that it cannot be expected to find perfect agreement because of the differences in the event selection, which is based on final-state leptons in our calculation, and the absence of corrections to the Z-boson decays in the on-shell calculations.References [39,40] state about −40% at p T,Z = 700 GeV at √ s = 14 TeV, which agrees with our result for √ s = 13 TeV at the percent level in spite of the different setups.
Figure 9 shows the transverse-momentum distribution of the µ + .The left panels compare the leading µ + from the [4µ] final state with the µ + from the [2µ2e] final state, the panels in the right column show the corresponding comparison with the subleading µ + .Recall that our ordering of muons into "leading" and "subleading" corresponds to the ordering with respect to the distances of the virtualities of µ + µ − pairs from the Z resonance, as described in Sec.3.2, but not to the muon p T , which is frequently used as well.We observe again that the observable is very sensitive to the event selection with characteristic differences between leading and subleading leptons.Especially at high transverse momenta p T , the spectrum of the leading muon in [4µ] is suppressed with respect to the spectrum of p T in [2µ2e], while the spectrum of the subleading muon in [4µ] is enhanced.The difference can be traced back to the impact of the ZZ signal and background contributions at large transverse momenta: The leading lepton belongs to the "more resonant" Z boson, and therefore, the contribution is in general dominated by the doubly-resonant signal contributions [cf.Fig. 1(a)].The main effect of the background contribution in Fig. 1(b) for large p T arises when the µ + is back-to-back with the three other charged leptons.As already observed for the related process of pp → WW → leptons [43], the impact of background diagrams on the p T spectrum of a single lepton can be as large as the doubly-resonant contribution in the TeV range.Since there is no preselection of the µ + in the [2µ2e] final state with respect to the resonance, the p Tµ + spectrum of [2µ2e] lies between the spectra of the leading and subleading muons.This behaviour is also reflected in the size of the purely weak corrections: Since the Sudakov enhancement is larger in doubly-resonant contributions than in singly-resonant contributions, the corrections reach at high p T about −45% for the leading µ + , about −35% for the µ + of the mixed flavour case, and about −30% in case of the subleading µ + .The photonic corrections give an almost constant negative offset of −1% to −2% for the mixed-flavour final state.The shape of the photonic corrections in the equal-flavour final state is very similar to the mixed-flavour case.For the subleading lepton, they amount to −1.5% to −3%, while for the leading lepton they stay between −1% and −0.5%. 3.5 3.5

Rapidity and angular distributions
The rapidity distributions in Fig. 10 do not show any significant dependence on the lepton pairing except for a small effect at the percent level in the forward direction with rapidities The EW corrections are independent of the lepton pairing as well, and their size is almost equal for both final states.The purely photonic corrections give, in good approximation, a constant negative offset of roughly −1%, reflecting the results of the integrated cross section in Tab. 1.The impact of the weak corrections is numerically largest in the central region with about −4.7% and less negative in the forward direction with about −3%.In Fig. 11, the distribution in the azimuthal-angle distance between the muons in the µ + µ − pair is shown.We observe a maximum for ∆φ µ + µ − → π that is reached in good approximation independently of the final state and the lepton pairing.This can be explained as follows: The azimuthal-angle-distance distribution is dominated in the whole range by events in the low-energy region around the threshold at M 4ℓ = 2M Z where the cross section receives the largest contribution from doubly-resonant contributions.Moreover, the t-channel nature of the doubly-resonant diagrams favours small transverse momenta of the Z bosons.As can also be seen in Fig. 9, most of the leptons have p T M Z /2 as a result of the decay of the Z bosons that move slowly in the transverse plane, i.e. the Z bosons decay almost isotropically in the transverse plane, without a large influence of boost effects from a transverse momentum of the  Z boson.For small ∆φ µ + µ − , the behaviour of the distribution is completely different.The µ + µ − pair in [2µ2e] as well as the subleading µ + µ − pair in [4µ] show some enhancement near ∆φ µ + µ − ∼ 0.2, while this enhancement is absent for the leading µ + µ − pair.This is a result of events with small µ + µ − invariant masses M µ + µ − < M Z .Owing to the µ + µ − pairing with respect to their invariant mass, low-mass µ + µ − pairs rarely occur as leading pairs, as can be seen in Fig. 7. On the other hand, low-mass µ + µ − pairs receive a much larger contribution from virtual photon (γ ⋆ ) exchange than nearly resonant µ + µ − pairs, leading to the observed peak structures at ∆φ µ + µ − ∼ 0.2.Note that these peaks are truncated on their left side by the lepton separation cut of Eq. (3.10).
In the dominant region of large ∆φ µ + µ − , the weak corrections resemble the size observed for the integrated cross sections above.For smaller ∆φ µ + µ − the negative corrections tend to increase in size, because this region in ∆φ µ + µ − on average requires more scattering energy to yield boost effects that turn the µ + and µ − directions away from the back-to-back configuration preferred at low energies.The effect of increasing negative weak corrections for smaller ∆φ µ + µ − is to some extent balanced in the [2µ2e] case and for the subleading µ + µ − pair in [4µ], because the weak corrections to the significantly contributing Zγ ⋆ production diagrams are smaller than the ones for ZZ production (cf.results on Zγ production shown in Ref. [77]).The photonic  corrections are generically small for the ∆φ µ + µ − distribution owing to the absence of collinear enhancements, because collinear final-state radiation does not change the directions of the radiating leptons significantly.The photonic corrections are most sizeable at ∆φ µ + µ − → π with −1.2%,−1%, and −1.5% for the mixed-flavour case, the leading lepton pair, and the subleading lepton pair, respectively.They decrease towards ∆φ µ + µ − → 0 where they reach −0.3%, −0.7%, and +0.2%, respectively.Figure 12 shows the distribution in the angle φ between the two Z-boson decay planes for the unequal-flavour [2µ2e] final state and the corresponding angle defined by the leading and subleading µ + µ − pairs of the equal-flavour [4µ] final state.The angle φ is defined in the CM system of the four final-state leptons by where , and k i denote the spatial parts of the four lepton momenta k i , of the µ + µ − e + e − and µ + µ − µ + µ − final states in the CM system.The dips in the distribution for coinciding decay planes, i.e. for φ ∼ 0 and φ ∼ π, are a consequence of the lepton-separation cuts in Eq. (3.10).These cuts remove collinear lepton configurations where the decay planes tend to be coplanar.The local minima around φ = π/2 and φ = 3π/2 can be understood from the superposition of contributions with different lepton helicities: If the two equally charged final-state leptons do have the same helicity, the distribution shows a maximum at φ = 0 and a minimum at φ = π, and vice versa for opposite lepton helicities.The corresponding distribution in W + W − pair production (c.f.Fig. 16 in Ref. [61]) exhibits a maximum at φ = 0 and a minimum at φ = π, as expected for purely left-handed leptons.We observe an enhancement of the ratio [2µ2e]/ (2[4µ]) above one for small angles φ and a reduction for angles φ close to π, both at the level of a few percent.Since an exchange of the leading and subleading µ + µ − pair does not affect φ, we attribute this effect to the interference between the two different classes of diagrams in [4µ] with different Z → µ + µ − pairings.For a qualitative understanding of the interference pattern, it is instructive to consider the limit of nearly coplanar Z decays (φ → 0, π).Since photon emission effects are not dominating this observable, the majority of the events have k µ , whose direction defines the intersection line of the two decay planes.In the vicinity of the coplanar configurations this line divides the event plane into two half planes.According to the definition (3.14) of φ, the two are parallel for φ → 0 and antiparallel for φ → π, i.e. k µ + 1 and k µ + 2 lie in the same half plane for φ → 0, but in different half planes for φ → π.On average, we thus find more parallel µ + 1 µ + 2 pairs for φ → 0 and more antiparallel pairs for φ → π.Since the matrix elements are antisymmetrized with respect to the exchange of µ , leading to the observed enhancement in the [2µ2e]/ (2[4µ]) ratio.This effect is not changed by the EW corrections.The weak corrections distort the φ distribution by about 3%, while the photonic corrections only uniformly contribute by −1%.

Collinear-safe versus collinear-unsafe observables
In Fig. 13, the different recombination schemes for the muon are illustrated for the four-lepton and two-lepton invariant masses.The recombination procedure only affects the photonic corrections.As a general pattern, all the radiative tails induced by final-state radiation off the charged leptons are strongly enhanced if collinear photons are not recombined with muons.The enhancement is due to the fact that the collinear logarithms are regularized by the muon mass rather than the size of the recombination cone.The effects can be best isolated in the M 4ℓ invariant-mass distribution (left panels of Fig. 13) which is not sensitive to the lepton pairing.While the absolute prediction is only shown for the collinear-unsafe case for the mixed-and equal-flavour final states, the relative EW corrections are plotted both for the collinear-safe and -unsafe cases.The results illustrate the impact of the number of muons excluded from recombination to the distribution: The maximum of the radiative tail below the ZZ threshold increases from about +30% with full recombination to more than +50% for excluding one muon pair ([2µ2e]) up to about +70% by excluding both muon pairs ([4µ]).For [4µ], the increase is twice as large as for [2µ2e], since the recombination effect scales with the number of collinear cones that are subject to the changes in the recombination.A similar behaviour is found for the other radiative tails at smaller values of M 4ℓ .An even stronger enhancement can be seen at the invariant mass of the µ + µ − pair (right panels) below the Z resonance where the relative correction increases from almost +60% to +140%.Note that above the resonance the effect from the collinear-unsafe treatment pushes the negative collinear-safe corrections even more negative.

Invariant-mass and transverse-momentum distributions
The production of Z-boson pairs at the LHC is interesting not only per se, as a signal process, but also constitutes an important irreducible background to Higgs production in the H → ZZ ⋆ decay mode.In order to assess the impact of this background on Higgs analyses, we impose the Higgs-specific cuts of Eqs.(3.12)- (3.11) in addition to the inclusive cut of Eq. (3.10).In Ref. [44] we already presented some important results of this study, however, restricted to the unequal-flavour final-state [2µ2e] and ignoring photon-induced channels.In the following we continue the discussion started there by comparing results for the [2µ2e] and [4µ] final states and considering further observables.Figure 14 illustrates the invariant-mass distribution of the four-lepton system at LO and the corresponding NLO EW corrections for both the [2µ2e] and the [4µ] final states.In each case, we observe a steep shoulder at the Z-boson pair production threshold at about M 4ℓ = 2M Z ≈ 182 GeV, which gives rise to a large radiative tail in the photonic corrections at smaller invariant masses.Though smaller in magnitude, a similar effect can be observed at around M Z + 2p T,min ≈ 103 GeV which is due to the transverse-momentum and invariant-mass cuts we impose on the charged leptons.Like in the inclusive setup, both the purely weak and the photonic corrections exhibit a sign change at the pair production threshold around M 4ℓ = 2M Z .The pattern of the EW corrections above the ZZ threshold is very similar to the inclusive setup with at most permille level differences between the [4µ] and the [2µ2e] case.The photonic corrections decrease in absolute size from approximately −2% at the threshold to about −1% at 1 TeV.The purely weak corrections constantly increase in absolute size reaching about −20% at 1 TeV.Also in the off-shell-sensitive region below the pair production threshold, the difference between the [4µ] and the [2µ2e] cases in the purely weak corrections is below the percent level.The radiative tails in the photonic corrections are up to 5% larger in the mixedflavour case.In contrast to the inclusive setup, the phase-space cuts of the Higgs-specific setup introduce a dependence on the lepton pairing even in otherwise symmetric observables like the four-lepton invariant mass.The difference seen in the photonic corrections is thus due to both the lepton pairing and the interference effects.At the Higgs-boson mass M 4ℓ = M H , the differences of the EW corrections with respect to the final states are, however, entirely negligible.The significant differences between the [4µ] and the [2µ2e] case in the off-shell-sensitive region   are, like in the inclusive case, a priori a LO effect.Note that the non-trivial sign change of the photonic corrections leads to significant cancellations between opposite-sign contributions below and above the ZZ threshold resulting in sub-permille effects in the total cross section (cf.Tab.1), although the individual photonic corrections can be sizable in distributions.
We also show the photon-induced contribution to the four-lepton invariant-mass distribution in the third panels from above in Fig. 14.Above the ZZ production threshold the corrections are at the level of one permille.In the off-shell region, the qγ and the γγ contribution are opposite in sign, at the level of 1%, and compensate each other to a large extent.The overall impact remains at the sub-percent level.We do not show the photon-induced corrections separately in the following plots.
In Fig. 15 the invariant-mass distribution of the µ + µ − system is shown for the [2µ2e] final state, as well as the ones for the leading and the subleading µ + µ − systems of the [4µ] final state.Due to the cuts of Eq. (3.12), the invariant mass of the leading muon pair in the equal-flavour final state is restricted to the range of 40−120 GeV.This cut leads in the [2µ2e] final state to   a little bump at 40 GeV.Moreover, the local maximum near M Z /2 in the M lead µ + µ − distribution of the inclusive setup is absent for the Higgs-specific setup, because the invariant-mass cut M 4ℓ > 100 GeV in Eq. (3.13) entirely removes the s-channel resonance at M 4ℓ = M Z .Near the Z resonance the photonic and weak corrections are very similar to the results in the inclusive setup (cf.Fig. 7).The distribution peaks at the resonance, M µ + µ − ∼ M Z , and receives large photonic corrections below that are due to final-state radiation effects.The weak corrections, on the other hand, are of the order of 5% and give rise to a change in sign near the Z-boson resonance.Above the resonance the EW corrections are qualitatively similar to the ones in the inclusive setup for both the leading and subleading µ + µ − pair.Below M µ + µ − ≈ 60 GeV, on the other hand, the missing s-channel resonance at M 4ℓ = M Z leads to significant changes.The difference is, as expected, most prominent in the leading lepton pair where the local minimum of the weak corrections at 45 GeV and the entire additional radiative tail of the photonic corrections are removed.While the EW corrections show sizeable deviations between the mixed-and equalflavour final states, the main differences are LO effects that can be attributed to the cuts and the lepton pairing.
Figure 16 depicts the transverse-momentum distribution of the µ + in the [2µ2e] final state together with the leading and the subleading µ + of the [4µ] final state, respectively.We once again remind the reader that the classification of leptons as "leading" or "subleading" refers to the criteria of Eq. (3.12), i.e. the leading muon is not necessarily the muon of highest transverse    distributions in the inclusive setup shown in Fig. 11 is absent in the Higgs-specific setup.As noted above, the enhancements in the inclusive setup at ∆φ µ + µ − ∼ 0.2 are mostly due to µ + µ − pairs of low invariant mass enhanced by the photon pole.The Higgs-specific selection cuts applied in the current setup remove such contributions, leaving us with azimuthal-angle distributions that still exhibit a similar rise towards ∆φ µ + µ − → π, but no longer peak at low values of ∆φ µ + µ − .Apart from the peak structure, the impact of weak corrections on normalization and shape of the azimuthal-angle differences in the Higgs setup is similar in size as in the inclusive setup.Purely photonic corrections are even more suppressed than in the inclusive case.In the Higgs-specific scenario the fraction of events with leading muon pairs close to the Z resonance is enhanced, while the one for subleading muon pairs is reduced (compare Figs. 7 and 15).As a consequence, the distribution of the leading muon pair is enhanced compared to the unequal-flavour case for large ∆φ µ + µ − , while the one of the subleading lepton pair is reduced.For small azimuthal-angle differences the situation is reversed.We show in Fig. 19 the distribution in the angle between the two Z-boson decay planes in the four-lepton CM system in the Higgs-specific setup. 4The distribution as well as the EW corrections closely resemble those of the inclusive setup shown in Fig. 12.The ratio [2µ2e]/(2[4µ]) is qualitatively similar, showing some excess over one for φ → 0, 2π and some deficit around φ ∼ π, where the latter is, however, more pronounced than in the inclusive setup.Owing to the asymmetric treatment of the leading and subleading µ + µ − pairs in the Higgs setup, we cannot attribute the deviations of the ratio from one to interference effects only.Finally, we remark that the φ distribution in the Higgs-signal process H → ZZ * → 4 leptons looks qualitatively similar to the distribution of direct ZZ production shown in Fig. 19, but the distortions by EW corrections are quite different [78,79]

Conclusions
The production of four charged leptons in hadronic collisions at the LHC is an important process class both for the investigation of the interactions between the neutral Standard Model gauge bosons and as background process to searches for new physics and to precision studies of the Higgs boson.In the confrontation of experimental data with theory predictions precision plays a key role.In this paper we have further improved the theory prediction by calculating the next-to-leading-order electroweak corrections to the production of µ + µ − e + e − and µ + µ − µ + µ − final states without any kinematical restrictions on the intermediate states.Our results are thus accurate to next-to-leading order in all phase-space regions, no matter whether they are dominated by two, one, or zero resonant Z bosons.Our numerical discussion of the corrections focuses on two different event-selection scenarios, one based on typical lepton-identification criteria only and another one that is specifically designed for Higgs-boson analyses.Since the Higgs-boson mass of about 125 GeV lies below the Z-pair threshold, the flexibility of our calculation, allowing intermediate Z bosons to be far off shell, is essential for the study of four-lepton production as background to the Higgs-boson decay H → ZZ ⋆ .Extending our earlier study [44] of the process pp → µ + µ − e + e − + X, we have investigated further observables and channels with photons in the initial state and included the process pp → µ + µ − µ + µ − + X. Generically, the next-to-leading order electroweak corrections consist of photonic and purely weak contributions displaying rather different features.Photonic corrections can grow very large, to several tens of percent, in particular in distributions where resonances and kinematic shoulders lead to radiative tails.These effects are significantly enhanced when observables within a collinear-unsafe setup are considered.While photonic corrections might be well approximated with QED parton showers, this is not the case for the weak corrections, which are typically of the size of −5% at intermediate energies and grow to multiples of −10% in the high-energy tails of invariant-mass and transverse-momentum distributions.Moreover, the weak corrections below the ZZ threshold distort distributions that are important in Higgs-boson analyses.On the other hand, contributions induced by incoming photons, i.e. photon-photon and quark-photon channels, turn out to be phenomenologically unimportant.Comparing the results on µ + µ − e + e − and µ + µ − µ + µ − final states, we find significant differences mainly in distributions that are sensitive to the assignment of µ + µ − pairs in the µ + µ − µ + µ − final state to intermediate Z bosons.Interferences in equal-flavour-lepton final states lead to deviations of up to 10% from the mixed-flavour case in off-shell-sensitive phase-space regions.Their effect is, however, in general hidden in the effects of the selection criteria for the lepton pairing.The relative electroweak corrections are widely insensitive to details of the lepton pairing, i.e. the selection of µ + µ − pairs affects observables at leading and next-to-leading order roughly in the same way.
The full calculation is available in the form of a Monte Carlo program allowing for the evaluation of arbitrary differential cross sections.The best possible predictions for ZZ production processes can be achieved by combining the electroweak corrections of our calculation with the most accurate QCD predictions available to date.Practically, this could be achieved, e.g., by reweighting differential distributions including QCD corrections with electroweak correction factors.In this way, an overall accuracy at the percent level can be achieved for integrated cross sections that are dominated by energy scales up to a few 100 GeV, where the theoretical uncertainty is completely dominated by QCD.We estimate the contribution of missing higherorder electroweak corrections on the integrated cross section to 0.5%.The impact of missing higher-order electroweak corrections grows in the high-energy tails of transverse-momentum and invariant-mass distributions where weak Sudakov (and subleading high-energy) logarithms are known to be large.In this kinematic domain, the size of this uncertainty may be estimated by the square of the relative electroweak correction.The inclusion of the known leading two-loop effects or a resummation of logarithmically enhanced contributions could reduce these theo-retical uncertainties.At the same time, multi-photon emission effects could be systematically taken into account by structure functions or parton showers.Such improvements are, however, left to future studies.For upcoming analyses of LHC data, next-to-leading order precision in electroweak corrections is certainly sufficient, and the remaining electroweak uncertainties are negligible compared to the larger uncertainties from missing QCD corrections and from parton distribution functions.

Figure 1 :
Figure 1: Sample tree-level diagrams contributing at O(α 4 ).The dominant qq channel (a,b) defines the LO contribution, while the photon-induced γγ channel (c) is counted as a correction.

Figure 2 :
Figure 2: Sample diagrams for the virtual EW corrections.Diagram types (a)-(d) are obtained by photon and Z-boson insertions between the fermion lines of the tree-level diagrams in Fig. 1(a).The remaining diagrams may involve couplings (f)-(i) or corrections to vertices (e) that are not present at LO.

Figure 3 :
Figure 3: Illustration of the splitting of EW corrections into purely weak (a) and photonic (b)-(d) contributions for the diagram type shown in Fig. 2(d).

Figure 4 :
Figure 4: Illustration of real photon radiation off the final-state leptons and off the initial-state quarks.
energy of 13 TeV are summarized in Tab. 1. Results are given for the inclusive selection cuts of Eqs.(3.9)-(3.10)and the Higgs-specific selection cuts of Eqs.(3.10)-(3.13).Together with the LO cross sections, the NLO EW corrections to the qq contribution and the relative contributions of the photon-induced channels δ γγ = σ γγ /σ LO qq and δ qγ = σ qγ /σ LO qq are shown.The qq contribution δ EW qq = ∆σ EW qq /σ LO qq is split into the purely weak and the photonic part, δ weak qq and δ phot qq , respectively, so that δ EW qq = δ weak qq + δ phot qq .

Figure 6 :
Figure 6: Invariant-mass distribution of the four-lepton system (upper panels), corresponding EW corrections (2nd panels from above), γγ and qγ contributions (third panels from above) for the unequal-flavour [2µ2e] and the equal-flavour [4µ] final states in the inclusive setup.The panels at the bottom show the ratio of the [2µ2e] and [4µ] final states.

Figure 7 :
Figure 7: Invariant µ + µ − -mass distribution (upper panels), corresponding EW corrections (middle panels), and ratio of the [2µ2e] and [4µ] final states (lower panels) in the inclusive setup.In the left column the equal-flavour case is binned with respect to the leading lepton pair, while the right column shows results for the subleading one.

Figure 8 :
Figure 8: Invariant-mass (left) and transverse-momentum distribution (right) of the µ + µ − pair (upper panels), corresponding EW corrections (middle panels), and ratio of the [2µ2e] and [4µ] final states (lower panels) in the inclusive setup.The equal-flavour case is binned with respect to the subleading lepton pair.

Figure 9 :
Figure 9: Transverse-momentum distribution of the µ + (upper panels), corresponding EW corrections (middle panels), and ratio of the [2µ2e] and [4µ] final states (lower panels) in the inclusive setup.The left panels compare the leading µ + from the [4µ] final state with the µ + from the [2µ2e] final state, while the panels in the right column show the corresponding comparison with the subleading µ + .

Figure 10 :
Figure 10: Rapidity distribution of the µ + (upper panels), corresponding EW corrections (middle panels), and ratio of the [2µ2e] and [4µ] final states (lower panels) in the inclusive setup.The left panels compare the leading µ + from the [4µ] final state with the µ + from the [2µ2e] final state, the panels in the right column show the corresponding comparison with the subleading µ + .

Figure 11 :
Figure 11: Azimuthal-angle difference between muons within the µ + µ − pair (upper panels), corresponding EW corrections (middle panels), and ratio of the [2µ2e] and [4µ] final states (lower panels) in the inclusive setup.The left panels compare the leading µ + µ − pair from the [4µ] final state with the µ + µ − pair of the [2µ2e] final state, the panels in the right column show the corresponding comparison with the subleading lepton pair.

Figure 12 :
Figure 12: Angle between the two Z-boson decay planes in the CM system for the unequalflavour [2µ2e] and for the leading and subleading µ + µ − pair of the equal-flavour [4µ] final state (upper panels), corresponding EW corrections (middle panels), and ratio of the [2µ2e] and [4µ] final states (lower panels) in the inclusive setup.

Figure 13 :
Figure 13: Comparison of different photon recombination schemes for the four-lepton and twolepton invariant-mass distributions.The upper panels show the absolute distributions and the lower panels the relative EW corrections.

Figure 14 :
Figure14: Invariant-mass distribution of the four-lepton system (upper panels), corresponding EW corrections (2nd panels from above), photonic contributions (third panels from above) for the unequal-flavour [2µ2e] and the equal-flavour [4µ] final states in the Higgs-specific setup.The lower panels show the ratio of the [2µ2e] and [4µ] final states.

Figure 15 :
Figure15: Invariant µ + µ − -mass distribution (upper panels), corresponding EW corrections (middle panels), and ratio of the [2µ2e] and [4µ] final states (lower panels) in the Higgs-specific setup.In the left column the equal-flavour case is binned with respect to the leading lepton pair, while the right column shows results for the subleading one.

Figure 17 :
Figure 17: Rapidity distribution of the µ + (upper panels), corresponding EW corrections (middle panels), and ratio of the [2µ2e] and [4µ] final states (lower panels) in the Higgs-specific setup.The left panels compare the leading µ + from the [4µ] final state with the µ + from the [2µ2e] final state, the panels in the right column show the corresponding comparison with the subleading µ + .

Figure 18 :
Figure 18: Azimuthal-angle difference between muons within the µ + µ − pair (upper panels), corresponding EW corrections (middle panels), and ratio of the [2µ2e] and [4µ] final states (lower panels) in the Higgs-specific setup.The left panels compare the leading µ + µ − pair from the [4µ] final state with the µ + µ − pair of the [2µ2e] final state, the panels in the right column show the corresponding comparison with the subleading lepton pair. .

Figure 19 :
Figure 19: Angle between the two Z-boson decay planes in the CM system for the unequalflavour [2µ2e] and for the leading and subleading µ + µ − pair of the equal-flavour [4µ] final state (upper panels), corresponding EW corrections (middle panels), and ratio of the [2µ2e] and [4µ] final states (lower panels) in the Higgs-search setup.

Table 1 :
2.3.Since we define σ LO qq from NLO PDFs (instead of LO PDFs), the resulting relative corrections δ LO cross sections for pp → µ + µ − e + e − +X and pp → µ + µ − µ + µ − +X with the relative corrections δ i = σ i /σ LO qq for the LHC at √ s = 13 TeV.Results are given for the inclusive setup