NLO EW and QCD corrections to polarized ZZ production in the four-charged-lepton channel at the LHC

Measuring the polarization of electroweak bosons at the LHC allows for important tests of the electroweak-symmetry-breaking mechanism that is realized in nature. Therefore, precise Standard Model predictions are needed for the production of polarized bosons in the presence of realistic kinematic selections. We formulate a method for the calculation of polarized cross-sections at NLO that relies on the pole approximation and the separation of polarized matrix elements at the amplitude level. In this framework, we compute NLO-accurate cross-sections for the production of two polarized Z bosons at the LHC, including for the first time NLO EW corrections and combining them with NLO QCD corrections and contributions from the gluon-induced process.

structure of the underlying resonances, both in the case of the Higgs and single-Z decay into four leptons, and in the case of the leptonic decay of two Z bosons. In fact, the Higgs decay into four leptons has been widely studied [38][39][40][41][42][43][44][45][46][47] and then measured with Run-2 data [1][2][3][4][48][49][50][51] with the purpose of determining the spin and parity of the Higgs boson as well as its CP properties via angular and energy correlations.
For the process pp → ZZ → e + e − µ + µ − , the possibility to fully reconstruct the final state allows to access observables that are sensitive to the polarization of Z bosons. This can be done in the Higgs-resonance region [52][53][54], where the polarization of the weak bosons may give hints of newphysics effects in the coupling of the Higgs to longitudinally-or transversely-polarized bosons. The polarization structure can also be studied in the so-called on-shell region [17], where two Z bosons are produced and then decay into charged lepton pairs: this is the target of this paper.
The polarization of W and Z bosons has been measured at the LHC in W + jet and Z + jet production [55][56][57][58] and in the decay of top quarks [59][60][61] with Run-1 data, and more recently also in WZ production [62] and in same-sign WW scattering [63] with the 13-TeV data. The upcoming high-luminosity LHC run is expected to enhance the sensitivity to the polarization of Z bosons produced in ZZ inclusive production and scattering [64,65].
The increasing experimental interest in measuring the weak-boson polarization in multi-boson processes has triggered a number of phenomenological studies for the LHC environment. The seminal works concerned the production of weak bosons in association with jets [66] and a number of other processes [67]. The effect of lepton cuts on polarization extraction was highlighted in Refs. [66,68]. More recently, a study of polarization-sensitive coefficients has been performed for W ± Z inclusive production [69,70]. The separation of polarization states at the amplitude level at leading order (LO) has been automated and investigated in vector-boson scattering [71][72][73] with the PHANTOM Monte Carlo [74] and in Higgs decays to four leptons [53,54]. The generation of polarized events in general multi-boson processes is available at LO in the MadGraph framework [75] via the spin-correlated narrow-width approximation [76][77][78]. The separation of polarized signals at the amplitude level has been extended to NLO QCD [79,80] and to NNLO QCD [81] for diboson production in the fully-leptonic channel. The growing interest in the polarization structure of di-boson production and the increased sensitivity to polarization observables have motivated us to extend these polarization studies to NLO EW.
The aim of this work is twofold. On the one hand, we present a general, well-behaved definition of polarized ZZ signals (both bosons have definite polarization state) including the decay modelling at NLO EW accuracy, which is achieved via the extension of the double-pole approximation (DPA) used in Refs. [79,80] to real radiation from leptonic decay products, and via the separation of polarization at the amplitude level in all parts of the calculation. On the other hand, we investigate the LHC phenomenology of polarized ZZ production including all LO and NLO contributions (both of QCD and of EW origin) in the presence of realistic kinematic selections, in order to address possible experimental analyses in this di-boson channel.
This paper is organized as follows. In Section 2, we point out the main difficulties that arise when separating resonant ZZ contributions from the non-doubly-resonant ones at NLO, and we describe in detail (Sections 2.1 and 2.2) our method to achieve this goal. The separation of polarizations is described in Section 2.3. The SM input parameters and the kinematic setups that are used in the numerical simulations are described in Sections 2.4 and 2.5, respectively. The integrated and differential results for polarized ZZ production at NLO EW are discussed in Section 3.1 for an inclusive setup. The combined NLO QCD and EW corrections, as well as the loop-induced gluonic corrections to the polarization structure of ZZ production are presented in Section 3.2 for a fiducial region that mimics the latest experimental selections for this process. In Section 4 we draw our conclusions.

Calculation details
We consider the production of a ZZ pair in the fully-leptonic channel, pp → e + e − µ + µ − + X , (2.1) at NLO accuracy, including both QCD and EW radiative corrections. At tree-level [O(α 4 )] only quark-antiquark and photon-photon partonic channels are present, while at NLO both photonquark and gluon-quark channels contribute to the real radiation at O(α 5 ) and O(α 4 α s ), respectively. Since the final state is charge neutral, the loop-induced gluon-gluon partonic channel [formally of order O(α 4 α 2 s )] gives a non-negligible contribution to this process thanks to large gluon luminosity in the proton. Sample LO diagrams are shown in Figure 1. Part of LO contributions to the qq and gg channels involve two Z bosons that are produced and then separately decay into charged lepton pairs: we refer to these contributions as doubly resonant or ZZ resonant. However, the full SM amplitude receives contributions also from diagrams that do not allow for the interpretation in terms of production and decay of two Z bosons: we refer to them as non-doubly-resonant. In order to separate the polarization modes of Z bosons, only ZZ-resonant diagrams must be selected. Nonetheless, this selection cannot be performed simply dropping non-doubly-resonant diagrams, because the EW gauge invariance would be lost. In order to recover gauge invariance, we use a pole-approximation technique [82][83][84][85][86][87], whose details are discussed in Sections 2.1 and 2.2. Another method that is often used to simulate the production and decay of unstable particles is provided by the narrow-width approximation [76][77][78]81].
In the γγ partonic channel, at most one s-channel Z propagator can be present in the SM amplitude. Therefore, this partonic process cannot contribute to the ZZ-resonant production neither at LO nor at NLO EW. It contributes to the full calculation only, and we have computed it including also its complete NLO EW corrections.
The selection of ZZ-resonant diagrams and a prescription to recover gauge invariance is needed also for NLO radiative corrections of virtual and real origin. Since we consider a final state with colour-less particles, both real and virtual QCD corrections of order O(α 4 α s ) only modify the structure of the initial state. This enables a simple separation between doubly-resonant and non-doublyresonant contributions, exactly in the same fashion as at LO. Sample real and virtual diagrams for QCD corrections to resonant ZZ production are shown in Figure 2. The details of the DPA used for NLO QCD corrections are described in Section 2.1.
The NLO EW corrections of order O(α 5 ) to the production of four charged leptons at the LHC are more involved than the QCD ones, as both real and virtual contributions modify the structure of the initial state as well as the structure of the final state, allowing for additional EW propagators (either on-shell or off-shell) connecting the initial and the final state. This non-factorized structure implies a more involved selection of ZZ-resonant diagrams, which is needed for the separation of polarization states. Sample tree-level and one-loop contributions to NLO QCD real and virtual corrections to resonant ZZ production. The shaded circle denotes arbitrary tree-level sub-diagrams.  Sample one-loop diagrams are shown in Figure 3. In order to extract ZZ-resonant contributions, only loop diagrams that can be interpreted as production and decay of two Z bosons, called factorizable, must be selected. They include both corrections to the ZZ production process, like the one of Figure 3(a), and corrections to the two decays Z → + − , = e, µ, like the one in Figure 3(b).
Also contributions with a virtual soft photon connecting either initial-and final-state particles [ Figure 3(c)], or two final-state particles coming from different Z-boson decays [ Figure 3(d)], called non-factorizable, are doubly resonant. In fact, such corrections are characterized by a universal structure that factorizes the doubly-resonant LO squared amplitude [88][89][90]. However, the impact of these corrections is expected to be small thanks to cancellations against the corresponding nonfactorizable real-radiation contributions [88,89]. Therefore, we exclude both virtual and real nonfactorizable corrections from the DPA calculation.
All other one-loop, non-doubly-resonant contributions need to be dropped. The virtual EW corrections to the γγ partonic channel involve doubly-resonant contributions (box and triangle fermion-loop diagrams), which however need to be interfered with tree-level diagrams which can only be non doubly resonant. Therefore, these loop diagrams with a γγ initial state enter the DPA calculation only at order O(α 6 ), in the same fashion as the gg channel gives doubly-resonant contributions at order O(α 4 α 2 s ). The real-radiation contributions to resonant ZZ production in the quark-antiquark channel, shown in Figure 4, can be divided in three classes. The first one concerns the radiation off the initial state [ISR, Figure 4(a)], the second and the third ones are characterized by the final-state radiation (FSR1, FSR2) of a photon off the e + e − and µ + µ − pair, respectively [see Figures 4(b)-4(c)]. It is important to note that squaring the sum of real-radiation contributions as in Figure 4 returns also non-factorizable corrections, arising from the interference of diagrams of different types. The soft singularities embedded in such contributions are cancelled by virtual non-factorizable contributions as those shown in Figures 3(c)-3(d). After combining real and virtual contributions, the non-factorizable corrections are expected to be smaller than the factorizable ones, given by the Figure 4. Real-radiation contributions to resonant ZZ production in the qq partonic channel at NLO EW. Figure 5. Real-radiation contributions to four-charged-lepton production in the γq partonic channel at NLO EW.
incoherent sum of squared ISR, FSR1 and FSR2 corrections [88,89]. In order to cancel infrared (IR) singularities properly, a careful selection of subtraction counterterms is needed, retaining only those (both integrated and unintegrated) that cure the singularities of doubly-resonant real corrections. In the full computation, also non-doubly-resonant structures contribute to real and virtual corrections, leading to a larger number of singular configurations that need to be subtracted. We perform the subtraction of IR divergences in the Catani-Seymour (CS) dipole formalism [91][92][93].
Photon-induced partonic channels open up at NLO EW. Sample diagrams for these real corrections are shown in Figure 5. The selection of ZZ-resonant diagrams in this case is much easier, as the final state (anti)quark cannot come from the decay of a Z boson and the doubly-resonant contributions are of ISR type only. The diagram topologies shown in Figure 5(a) and 5(b) can embed singular configurations in which the emitted (anti)quark is collinear either to the initial-state photon [ Figure 5(a)] or to the initial-state (anti)quark [ Figure 5(b)]. It is important to select only the subtraction counterterms that absorb the first kind of configurations, since the underlying LO contribution of the second configuration is non doubly resonant (γγ induced).
The selection of doubly-resonant diagrams in Born-level, virtual, and real-radiation contributions is essential to enable the separation of polarizations in Z-boson propagators. In the DPA [85,86], the two resonant bosons must be projected on their mass shell in order to restore EW gauge invariance. This is achieved via suitable on-shell projections that are described in Sections 2.1 and 2.2.
Before going into the details of the DPA used for this calculation, we point out that the neutral electric charge of the Z boson simplifies the selection of doubly-resonant structures both in real and in virtual corrections, as photons cannot be radiated off the resonant Z-bosons. In the case of W bosons, additional care must be taken in the separation of initial-state and final-state radiation contributions, in particular concerning the definition of subtraction counterterms. We postpone the treatment of charged resonances at NLO EW to future work.

A double-pole approximation for Born-like doubly-resonant structures
After selecting doubly-resonant diagrams, it is crucial to apply an on-shell projection to all LO and NLO matrix elements for ZZ production, retaining the off-shell-ness of the Breit-Wigner modulation in the denominator of Z-boson propagators. This has to be pursued with special attention to the subtraction counterterms, devising suitable projections that do not spoil the cancellation of IR singularities. A general NLO (either QCD or EW) observable receives the following contributions in the dipole formalism [91][92][93], where B, V, R are Born, virtual and real squared matrix elements, respectively, including spin and colour averages as well as symmetry factors, while 1/(2s) is the flux factor. The D term is the sum of subtraction dipoles that appears both in its 4-dimensional, unintegrated version and in its D-dimensional, integrated version. The quantity C represents the collinear counterterm that cancels left-over collinear singularities of initial-state type. The arguments of the δ symbols indicate whether the observable ξ is evaluated with n-body or with (n + 1)-body kinematics (with D = 4). Note that the subtraction dipoles depend on the mapped n-body kinematicsΦ n , which is derived from (n+1)-body kinematics via CS mappings [91], and on the radiation phase-space measure Φ rad . While we keep n variable in the formulas, n = 4 for ZZ production with leptonic Z decays. Let us consider all NLO contributions which embed a doubly-resonant sub-amplitude with two two-body decays, namely Z (→ e + e − ) Z (→ µ + µ − ). This is the case for Born matrix elements, virtual corrections, integrated dipoles, and real-subtracted initial-state radiation contributions.
We label the pole approximation applied to these contributions with DPA (2,2) . Let us call q 1 , q 2 the two off-shell Z-boson momenta in the laboratory system, andq 1 ,q 2 the corresponding on-shell-projected momenta. We chooseq 1 ,q 2 such that and the spatial direction of the two bosons in the di-boson centre-of-mass frame is preserved [86]. This procedure is valid for M e + e − µ + µ − > 2M Z . In the four-charged-lepton channel, this constraint can be achieved also experimentally thanks to the reconstructable final state [10]. We are left with the modification of the momenta of the four final-state leptons, according to the on-shell kinematics of the two projected bosons. Considering the Z (q 1 ) → e + (k 1 ) e − (k 2 ) decay, we choose to preserve the energy fraction and the spatial direction of the electron (and of the positron) in the corresponding Z rest frame: -we boost k 1 , k 2 into the rest frame of the off-shell Z boson (momentum q 1 ), where the two lepton momenta haven k1 andn k2 as spatial directions; -we rescale the lepton energies according to the on-shell-ness of the projected Z boson, i.e. we multiply the off-shell energies by M Z / q 2 1 ; -we set the spatial directions in the rest frame of the on-shell Z bosons equal to the off-shell ones (n k1 ,n k2 ); it is easy to see that in this reference frame the two projected momenta of the leptons readk 1 = (1,n k1 )M Z /2 andk 2 = (1, −n k1 )M Z /2; -we boost backk 1 ,k 2 into the laboratory frame according to the projected Z-boson kinematics (momentumq 1 ).
The same procedure is then applied to the decay products of the second Z boson (Z → µ + µ − ). The described DPA can be safely applied to factorizable virtual corrections and integrated dipoles, as well as to the real-subtracted corrections of ISR type. Since the QCD corrections only modify the ZZ production process, all contributions at order O(α 4 α s ) can be treated with DPA (2,2) . Note that only two subtraction dipoles are needed to cancel IR singularities of QCD origin, and they are both of initial-initial type, therefore the modification of the momenta of finalstate leptons (colourless) does not interfere with the phase-space mappings used in subtraction counterterms. Both subtraction terms D and initial-state collinear counterterms C [in Eq. (2.2)] need to be computed with DPA (2,2) . This DPA has already been used to compute NLO QCD corrections to polarized W + W − and W + Z production [79,80].
In the computation of real EW corrections, only ISR contributions [see Figure 4(a)] can be computed with DPA (2,2) , following the same procedure that has been used in the NLO QCD case, but taking care of one additional subtlety. Differently from QCD corrections, the spectator for the ISR subtraction dipoles may be an initial-state or a final-state particle. However, only an initialstate spectator must be selected to properly cancel the IR singularities of ISR type only in DPA (2,2) . Furthermore, treating the (n + 1)-body kinematics and the mapped kinematics (in the sense of CS mappings) with the same DPA technique ensures that the subtraction is not spoiled. Adapting the notation of Eq. (2.2), the real-subtracted contribution reads (with off-shell kinematics), whereΦ n is the mapped n-body kinematics, according to an initial-initial CS mapping, while k i is the radiated-photon momentum. Both the emitter and the spectator are initial-state particles (I 1 , I 2 ). After the application of DPA (2,2) , whereΦ n andΦ n+1 represent the on-shell-projected momenta of the n-and (n + 1)-body phase space. The DPA (2,2) does not modify the photon radiation, neither in the (n + 1)-body kinematics, nor in the radiation phase-space that appears in the dipoles. Since the initial-initial CS mappings modify all final-state particles with the same boost, the application of DPA (2,2) to the mapped kinematics does not interfere with the CS-mapping procedure, as the total momentum of the four leptons is preserved in the DPA, i.e. the on-shell projection of DPA (2,2) and the initial-initial CS mappings commute.
Note that the DPA (2,2) can also be used to treat the photon-induced real contributions, since the corresponding doubly-resonant contributions do not allow for final-state radiation.
We stress that, since also integrated dipoles are treated with DPA (2,2) , the needed correspondence between subtraction and integrated dipoles is maintained for ISR contributions.

A double-pole approximation for final-state radiation
The treatment of final-state-radiation contributions to EW real corrections is more delicate than the treatment of ISR ones. Let us consider the case in which a Z boson decays into e + (k 1 ) e − (k 2 ) γ(k 5 ), while the other one decays into µ + (k 3 ) µ − (k 4 ), corresponding to the diagram topology shown in Figure 4(b). For such corrections, we have devised an on-shell projection (labelled DPA (3,2) ) that differs from DPA (2,2) but follows the same driving idea in terms of preserved quantities. For the projection of the momenta of Z bosons, we apply the same method as for DPA (2,2) , with the difference that the momentum of the first Z boson is now given by the sum of three momenta, After this step, we haveq 2 1 =q 2 2 = M 2 Z andq 1 +q 2 = q 1 + q 2 . The momenta k 3 , k 4 can be projected following the prescription of DPA (2,2) (see Section 2.1), as a two-body decay is understood for the second Z boson. For the three-body decay, -we boost k 1 , k 2 , k 5 into the rest frame of the off-shell Z boson (momentum q 1 ), where the three particles haven k1 ,n k2 ,n k5 as spatial directions; -we rescale the lepton and photon energies according to the on-shell-ness of the projected Z boson, i.e. we multiply the off-shell energies by M Z / q 2 1 ; -we set the spatial directions of the final-state momenta in the rest frame of the on-shell Z boson equal to the off-shell ones (n k1 ,n k2 ,n k5 ), which allows to preserve the relative angles among the three particles; -we boost backk 1 ,k 2 ,k 5 into the laboratory frame according to the projected Z-boson kinematics (momentumq 1 ).
It is easy to check that this prescription does not violate momentum-conservation and preserves the ratio of Lorentz invariants constructed with the three decay products, It is crucial that the energy fractions and the relative angles are preserved in the three-body decay, as the radiated photon may be soft and/or collinear to the electron or to the positron. The DPA (3,2) guarantees that the soft-ness and/or collinearity of the photon is preserved. The subtracted FSR1 real corrections read [with off-shell (n + 1)-body kinematics]: whereΦ n is the mapped phase-space according to a final-final CS mapping. The contributing dipoles involve as emitter and spectator either the electron or the positron, which form with the photon the decay products of the first Z boson. At variance with the DPA (2,2) , the photon momentum k 5 is modified by the DPA (3,2) . Therefore the same DPA must be applied to the radiation phase-space and to the factorized phase-space in the CS-mapped process to guarantee the proper subtraction of IR singularities. This could spoil the correspondence between the subtraction dipoles appearing in Eq. (2.8) and their integrated counterparts, as the former are projected with DPA (3,2) , while the latter can only be projected with DPA (2,2) (as they feature n-body kinematics). We will see in the following that this is not the case. If no DPA is applied, the CS mapping for D 15,2 reads, where s ab = 2k a · k b and s 125 = s 12 + s 25 + s 15 . In order to use the DPA for the subtraction dipoles, we cannot simply apply the DPA (2,2) to the mapped momenta, otherwise we would implicitly assume that the photon momentum is untouched, which is not the case when applying DPA (3,2) to FSR1 real contributions. A consistent procedure is to apply the CS mapping to the (n + 1)-body kinematics that has already been projected with the DPA (3,2) , k 1 =k 1 +k 5 −s 15 s 12 +s 25k 2 ,k 2 =s 125 s 12 +s 25k 2 ,k 3 =k 3 ,k 4 =k 4 . (2.10) The final-final CS mapping does not modify the virtuality of the Z boson [q 2 1 = (k 1 +k 2 +k 5 ) 2 = M 2 Z ]. Furthermore, according to Eq. (2.7), the Lorentz invariants that appear in Eq. (2.10) coincide with those of Eq. (2.9).
The more involved character of the DPA (3,2) with respect to DPA (2,2) originates from the fact that the on-shell projection of DPA (3,2) and the final-final CS mappings are not commuting, unlike the case of DPA (2,2) application and initial-initial CS mappings. In order to further prove the decent behaviour of the DPA (3,2) , we check that its application to the subtracted final-stateradiation corrections does not spoil the correspondence between unintegrated and integrated dipoles. Let us consider the D-dimensional structure of a dipole with emissus i (photon), emitter j (e + or e − ) and spectator k (e − or e + ), recalling that the three particles are decay products of the same Z resonance, where Q j , Q k are the electric charges (in unit of e) of the emitter and spectator, respectively, and the phase-space measure depends on the variables (2.12) The kernel that has to be integrated is [91,92] s|V ijk (y, z)|s = 8πα µ 2 Q 2 We can then re-write Eq. (2.11) as where we have used Q j Q k = −1 and Q 2 j = 1 for the e + e − final state. The 4-dimensional version of the subtraction dipole reads while its integration in D dimensions (D = 4 − 2 ) gives dΦ n α 2π For ZZ production, we have n = 4 and the LO doubly-resonant structure can be written as where the barred kinematics depends on the choice of the emitter j and spectator k that is used in the dipole phase-space mapping. The mapped momenta are labelled by the indices of the emitteremissus pair j5 and the spectator k. The (n + 1)-body part of the NLO correction for the doublyresonant FSR1 contribution reads (k 5 is the photon momentum): where y , z are defined as y, z up to the exchange of the momenta k 1 and k 2 . Note also that, according to the final-final CS mapping,s (15,2) =s (25, 1) = s 125 ands (15,2) 34 =s (25, 1) 34 = s 34 . We observe that the action of DPA (2,2) on the mapped n-body kinematics would result in the modification of the radiation measure dΦ rad (k i ) as well as of the subtraction dipole, because the radiation momentum would be left untouched by the on-shell projection. On the contrary, if DPA (3,2) is applied to the factorized (n + 1)-body kinematics, both the integration measure and the dipoles do not change thanks to Eq. (2.7). Furthermore, this way the three partons are correctly projected to reconstruct an on-shell Z boson. The application of DPA (3,2) to both the complete and the factorized (n + 1)-body kinematics gives the DPA-regulated version of Eq. (2.17): where we have used the fact thatỹ = y andz = z, since these variables only depend on ratios of (n+1)-body invariants that are preserved according to Eq. (2.7). We stress again that the momenta in the numerators N B are obtained from (n+1)-body kinematics performing first the DPA and then the CS mapping (the procedures do not commute). Owing to Eq. (2.7) the integration measure dy dz is left untouched by the DPA, and the integrated dipoles do not need to be modified but give exactly the expression of Eq. (2.15). This is the case because the DPA (2,2) applied to the n-body kinematics of the integrated dipoles proceeds in the same fashion as DPA (3,2) . The analogous procedure, which we label DPA (2,3) , is applied to FSR2 real contributions shown in Figure 4(c) and the corresponding two subtraction dipoles.
The entire DPA formalism that we have presented for ZZ inclusive production, in particular the on-shell projections for Born-like and final-state radiation structures, can be straightforwardly extended to more complicated processes involving Z bosons, e.g. ZZ scattering, or other neutral resonances.
As a last comment of this section, we stress that the DPA can be applied rather straightforwardly to same-flavour leptonic decays, e.g. pp → Z(e + e − )Z(e + e − ) + X. The intrinsic ambiguity in identifying the resonances can be solved by associating the identical final-state leptons to specific Z bosons and accounting for the proper symmetry factor in the resonant amplitudes. A DPA for a process with such an ambiguity has, for instance, been used in Ref. [94] (see discussion in Section 2.2 there).

Separating polarizations
The application of the DPA to all contributions of the NLO (EW or QCD) computation enables us to isolate in a gauge-invariant way the resonant contributions to ZZ production from the nondoubly-resonant background, which can be estimated as the difference between the full result and the doubly-resonant contributions in the DPA, and is expected to Thanks to the DPA, it is possible to separate the polarization of Z bosons at the level of doubly-resonant amplitudes, following the procedure that has already been applied to vector-boson scattering [71][72][73], to W + Z and W + W − production at NLO QCD [79,80], and to Higgs decays [53,54]. After selecting doubly-resonant diagrams, the following structure (written in the 't Hooft-Feynman gauge) characterizes the SM amplitude, where P is the production part of the amplitude, while the D terms (not to be confused with the subtraction terms in Sections 2.1 and 2.2) are the parts corresponding to the two-or three-body decays of the Z bosons. Writing the metric tensors in terms of polarization vectors, Eq. (2.19) becomes where for both Z bosons the polarization sum runs over the longitudinal (L), left-handed (−), and right-handed states (+). Since we consider massless final-state leptons, the gauge terms that are proportional to the boson momenta vanish for all LO and NLO contributions. In fact, it is easy to see that the Z-boson propagators are always contracted with massless fermionic currents, both in tree-level and in one-loop doubly-resonant amplitudes. In the general case, the contributions of the terms proportional to boson momenta in the Z-boson propagators cancel against diagrams where the Z bosons are replaced by the corresponding would-be Goldstone bosons. From Eq. (2.20) it is clear that applying the DPA and squaring the doubly-resonant (unpolarized) amplitude, one obtains well-defined polarized squared matrix elements and off-diagonal contributions which are called interference terms: Replacing the unpolarized squared matrix element with a single {λ, λ } term in the sum of Eq. (2.21), we can generate events for the production and decay of two Z bosons with definite polarization states (λ and λ , respectively).
In order to reduce the number of different contributions and to minimize the impact of interference terms, we have defined the transverse mode (T) as the coherent sum of the left-and right-handed modes (including the left-right interference).
The polarization vectors in Eq. (2.20) must be defined in a certain reference frame. Common choices for di-boson production are the laboratory (LAB) frame [69,79] and the centre-of-mass (CM) frame of the two bosons [70,80]. The CM-frame choice has been shown [80] to be more natural than the LAB one for di-boson production and has been used in the latest ATLAS analysis of W ± Z production [62]. It is essential that the same definition of polarization vectors is used in all parts of the NLO calculation. In the case of the CM-frame definition, for final-state-radiation contributions to NLO EW corrections, the CM frame is the rest frame of the system formed by the four charged leptons and the photon, which is crucial to have a functioning subtraction of IR singularities. In Section 3.1 we compare the results in both the CM-and the LAB-frame definitions, while in Section 3.2 we only consider polarizations defined in the CM frame.
Before presenting the numerical results of our calculation, it is worth recalling that the polarization structure of ZZ production at the LHC is expected to be rather different from the one of W + W − [79,81], in spite of the same contributing partonic channels. At LO, the differential crosssections for on-shell ZZ production in qq annihilation can be written in terms of the CM energy squared s and the scattering angle θ in the CM frame, where c L,q , c R,q are the left-and right-handed couplings of the Z boson to quarks q, and s w , c w are the sine and cosine of the EW mixing angle. At a given scattering angle, the expansion in powers of s gives a LL cross-section that is suppressed by 1/s 2 w.r.t. the TT one and by 1/s w.r.t. the mixed one. The suppression of the LL signal can be explained using the Goldstone-boson equivalence theorem [95][96][97][98], which relates the amplitudes for longitudinal vector-boson production to those for would-be-Goldstone-boson production at high energies of the vector bosons. Since the couplings of the Goldstone bosons to massless fermions vanish, there are no Feynman diagrams for the production of a pair of neutral would-be Goldstone bosons at LO. Consequently, the amplitude for the production of a pair of longitudinal on-shell Z bosons vanishes in the high-energy limit as 1/s, and the corresponding cross-section is suppressed w.r.t. the one for purely-transverse and mixed polarization states at large energy [99,100]. In contrast, the amplitude for longitudinal W + W − production is not suppressed in the same way, since the amplitude for the production of a pair of charged Goldstone bosons is finite, thanks to the diagram with an s-channel exchange of a neutral EW boson. These on-shell results indicate a very small purely-longitudinal ZZ signal, which is expected also when including off-shell effects and radiative corrections.

Numerical tools and input parameters
We study the production of e + e − µ + µ − in LHC proton-proton collisions at a centre-of-mass energy of 13 TeV with NLO QCD and EW accuracy within the SM.
The computation is performed with MoCaNLO, a multi-channel Monte Carlo integrator that has already been used for several multi-boson LHC processes, and in particular for the NLO QCD corrections to polarized W + W − and W + Z production [79,80] via the well-known expressions [82], (2.24) The EW coupling α is fixed in the G µ scheme [86], The top-quark and Higgs-boson masses are set to The complex-mass scheme is utilized for the treatment of unstable particles [85][86][87]105]. The subtraction of infrared singularities is achieved in the dipole formalism [91][92][93] both for QCD and EW corrections. The values of parton distribution functions (PDFs) and the running of α s are passed to MoCaNLO via the LHAPDF6 interface [106]. We choose the NNPDF31 nlo as 0118 luxQED PDF set [107,108], which understands α s (M Z ) = 0.118 and includes also the modelling of the photon luminosity in the proton [108]. The same PDF set is used for both LO-and NLO-accurate results. Initial-state collinear singularities are absorbed in PDFs via the MS scheme. The factorization and renormalization scales are simultaneously set to the Z-boson pole mass, i.e. µ R = µ F = M Z .

Kinematic selections
Photon-and jet-clustering is performed with the anti-k T algorithm [109]. Recombination is done for particles with a maximum rapidity of 5. In particular, photons are recombined with charged leptons = e, µ if QCD jets can only arise from real initial-state radiation and cannot generate singular configurations with the final-state leptons. Their kinematics is not constrained by any angular-distance cut. Following a similar approach as in Ref. [79], we have considered two sets of cuts, a first one (label INC) that is rather inclusive in the kinematics of final-state leptons and a second one (label FID) that mimics the fiducial region used in the latest ATLAS measurement [10]. The INC setup just involves • a technical cut on the transverse momentum of each charged lepton, p T, > 0.001 GeV; • an invariant-mass cut on the two pairs of opposite-sign, same-flavour leptons, The second cut is essential to enhance the on-shell production of two Z bosons, i.e. to reduce the photon contamination and other non-doubly-resonant contributions. The FID setup involves • a minimum transverse-momentum and a maximum-rapidity cut for the electron and the positron, p T,e ± > 7 GeV, |η e ± | < 2.47; • a minimum transverse-momentum and a maximum-rapidity cut for the muon and the antimuon, p T,µ ± > 5 GeV, |η µ ± | < 2.7; • a minimum transverse-momentum cut on the leading, p T, 1 > 20 GeV, and on the sub-leading charged lepton, p T, 2 > 10 GeV (sorted according to transverse momentum); • a rapidity-azimuthal-angle separation between any two leptons, ∆R > 0.05; • an invariant-mass cut on the two pairs of opposite-sign, same-flavour leptons, |M + − − M Z | < 10 GeV; • a minimum invariant-mass cut on the four-lepton system, M 4 > 180 GeV.
We have checked numerically that the effect of the M 4 > 180 GeV cut is to enhance the on-shell ZZ production even in the absence of the |M + − − M Z | < 10 GeV cut: under these conditions the non-doubly-resonant background would account for 15% of the full result. The joint effect of M 4 > 180 GeV and |M + − − M Z | < 10 GeV cuts reduces the non-doubly-resonant background to less than 1% of the full cross-section.

Results
We now present the numerical results of our computation. In Section 3.1 we consider the INC setup and NLO EW accuracy. These results provide both a validation of the calculation and a first assessment of how the polarization structure of ZZ production is modified by EW corrections.

Inclusive phase-space: NLO EW results
The results presented in this section have been obtained in the INC setup described in Section 2.5. Before tackling the polarized signals, we discuss the differences between the full calculation and the DPA-approximated one for unpolarized production. The proper functioning of the subtraction of EW IR singularities in the DPA calculation has been tested by varying the subtraction parameter α dip [110] between 10 −2 and 1. The sum of real-subtracted and integrated-dipole contributions has been found to be independent of α dip , demonstrating that the combined application of DPA techniques presented in Section 2.1 and Section 2.2 works in the desired manner. In Table 1 we show the various contributions to the NLO EW cross-section, separating them into contributions of partonic channels. The LO DPA result is roughly 1.2% smaller than the full one, in agreement with the expected difference of order Γ Z /M Z = 2.7% (intrinsic accuracy of the DPA). Note that when performing the DPA for Z bosons, an on-shell constraint (|M + − − M Z | < 10 GeV in our case) is crucial to suppress the photon contamination to the leptonic decays, which is included in the full but not in the DPA result [72]. At LO the full result receives a very small contribution from the γγ-induced partonic process.
The EW corrections in the qq partonic channel dominate the total NLO EW correction, and the DPA underestimates (in size) the full result by 2.2%. This difference comes mostly from the large virtual and integrated-dipole contributions (−2.610 fb in the full calculation) and is enhanced by the relative sign between virtual and real corrections. Note that the real contributions are small and positive (0.337 fb in the full calculation) and are approximated by the DPA within 0.9%.
The γq and γq partonic channels participate in the calculations as real corrections that involve only initial-state singularities. A marked difference is found between the DPA numerical result and the full one for such partonic channels, since the real-subtracted contributions in the full case embed additional dipoles compared to those that are included in the DPA, corresponding to the non-doubly-resonant contributions with underlying γγ as Born partonic channel [see Figure 5 (b)]. This causes the presence of additional integrated dipoles from the γγ Born contribution that are absent in the DPA calculation. Note that the EW corrections to the γγ-induced process account for less than 0.01% of the LO cross-section. Summing all contributions, the NLO EW cross-section in the DPA is 1.1% smaller than the full one, indicating a small non-doubly-resonant background from photon contamination and missing off-shell effects, in particular in the region between 162 GeV < M 4 [36]. If a looser invariant-mass cut is applied, the impact of EW corrections is smaller [18,19]. The difference between the full and the unpolarized DPA calculations implies that the nondoubly-resonant background is very small (≈ 1%), and on-shell ZZ production dominates the four-lepton production in the considered setup. Larger effects may arise in differential distributions, as observed in W + W − and WZ production [79,80] where at large transverse momentum single-resonant contributions become sizeable, generating large non-doubly-resonant backgrounds. However, in the four-charged-lepton channel the single-resonant contributions are much more suppressed than in the presence of final-state neutrinos, thanks to the possibility of constraining both Z bosons to be almost on-shell via physical invariant-mass cuts.
We now consider the ZZ signals where both bosons have definite polarization states (L or T). In Table 2 we show the LO and NLO EW integrated cross-sections for the unpolarized and polarized process. Polarizations are defined either in the CM or in the LAB frame. The interference among polarization states is evaluated as the difference between the unpolarized DPA result (which, by definition, contains also interference terms) and the sum over doubly-polarized results. For integrated cross-sections, interferences are found to be at the 0.2% level both at LO and at NLO EW.
It has been shown [80] that the spin correlations between the two bosons can be sizeable in the CM-frame definition of polarization and enhanced if both bosons are in a definite polarization state. Therefore, the sum over doubly-polarized modes is expected to differ from the unpolarized DPA result more than the sum over singly-polarized modes, i.e. larger interference effects are present.
Up to the correlation effects mentioned above, vanishing interferences are expected for the inclusive case at LO, where the two-body decay of each Z boson enables an exact analytic expression of the unpolarized DPA cross-section in terms of a sum over polarized contributions [66,67]: where θ * + is the decay angle of the antilepton in the rest frame of the corresponding Z boson (with momentum equal to the sum of the two lepton momenta) computed w.r.t. the Z direction in the CM (LAB) frame for polarization vectors defined in the CM (LAB) frame. The coefficients f L , f + , f − are polarization fractions, and c L, , c R, are the left-and right-handed couplings of the Z boson to leptons. This expression can be derived assuming that the complete phase-space of decay products is available, such that interference terms vanish upon integration over the azimuthal decay angle [71]. Projecting the cos θ * + distribution onto suitable polynomials in cos θ * + [66,67,71], it is possible to extract the polarization fractions for the Z boson and compare the result with the one obtained by directly simulating polarized Z bosons with the Monte Carlo. In the INC setup, this comparison yields a very good agreement for the LO predictions as well as for other n-body contributions (virtual, integrated dipoles) to the NLO EW cross-section. While doubly-polarized polarization fractions cannot be obtained from projections on parametrizations of unpolarized cross-sections like Eq. (3.1), their definition is straightforward based on Eq. (2.21). The LO interference estimated from the sum of the thus defined doubly-polarized Monte Carlo predictions is not identically zero (as it would be summing singly-polarized signals) due to the above-mentioned correlation effects.
In the presence of final-state real corrections, Eq. (3.1) does not hold anymore, as one Z boson decays into three particles, therefore projecting the cos θ * + distributions in the same way as at LO gives results that have nothing to do with polarization fractions [69]. Furthermore, in the presence of photons radiated off decay leptons, the interferences among polarization states are not vanishing anymore, even if the complete phase-space is available for the three-body decay. For example, if we consider a polarized Z boson that decays into e + e − γ, the polarization fractions from Monte Carlo simulations (f MC The discrepancy between the two evaluations of polarization fractions is large, indicating that real contributions to polarized signals can only be computed with polarized amplitudes in the Monte Carlo, even in a very inclusive setup like the one that is understood here. Note that the sum of polarization fractions extracted with cos θ * e + projections gives 1 by definition, while the sum of polarization fractions (ratios of polarized cross-sections over the DPA unpolarized one) computed with the Monte Carlo suggests that a 2% effect comes from interferences. However, although interferences in FSR1/2 corrections are at the 2% level, when summing all EW corrections to ZZ production this effect is hardly visible in the NLO EW total cross-section, since the real corrections are small compared to the virtual ones. This explains why even at NLO EW the overall interferences are small (0.16% and 0.20% for polarizations defined in the CM and the LAB frame, respectively).
It is worth noting that in both polarization definitions the relative NLO EW corrections stay between −10% and −12% for all modes, with slighly more sizeable corrections for the LL and TT modes than for the mixed ones. In both definitions, the polarization fractions are hardly modified by NLO EW effects: the modes with at least one longitudinal boson are slightly enhanced, while the TT mode is mildly diminished. In the CM-frame definition, the LL and the mixed modes are of the same order of magnitude (5.8% and 11.6% respectively), while the TT mode is dominant. This situation is completely changed in the LAB definition, as the LL mode accounts for just 1.6%, while the sum of mixed contributions gives 57.8% of the total. This difference, in a similar fashion as in W + Z production [80], is due to the fact that part of the LL contributions defined in the CM frame (where the polarization vectors just depend on one Z-boson momentum, thanks to the back-to-back kinematics) become transverse when boosting to the LAB frame, giving a more involved dependence on kinematic quantities. This has also the effect of reducing the TT component in favour of the mixed ones.
The polarized results for differential cross-sections can embed effects that are rather different from those obtained at the integrated level, in particular concerning the interference effects. In Figures 6-7, we study the distributions in a number of LHC observables, focusing on the shapes and normalizations of differential cross-sections in top panels, the impact of EW corrections in middle panels, and the impact of non-doubly-resonant (grey curves) and interference effects (grey vs. magenta curves) in bottom panels.
In Figure 6 we consider the decay angle of the positron in the Z rest frame, which has been introduced in Eq. (3.1). This angular variable is expected to be highly sensitive to the polarization of the first Z boson. At NLO EW the Z boson is understood as the system of the two dressed leptons (e + e − ), i.e. after photon recombination, in order to ensure the IR safety of the variable. The angle θ * e + can be computed w.r.t. the Z-boson direction in the CM frame or w.r.t. the Z-boson direction in the LAB frame. For both definitions of θ * e + , the non-doubly-resonant background effects are flat and equal to those found at the integrated level. If θ * e + is computed w.r.t. the Z-boson direction in the same frame where polarization vectors are defined, the analytic expression for the differential distribution is known at LO [see Eq. (3.1)] and interferences are expected to vanish. This is the case in Figures 6(a) and 6(d), where also at NLO EW the interferences are at the few-permille level owing to small FSR1/2 real contributions. Both in Figure 6(a) and in Figure 6(d), the LL and LT distributions display the same shape, which has a sin 2 θ * e + functional dependence on the angle. In Figure 6(a), the TL and TT distributions exhibit the same symmetric shape, because of left-right symmetry of Z-boson production in the CM frame. 1 On the contrary, in Figure 6(d) the TL distribution is not symmetric, due to a different left-and right-polarization content for the first boson in the presence of a recoiling longitudinal boson. The difference in shape between the TL and the TT curves indicates that the polarization of the first boson is affected differently by different polarization states of the second boson.
If both the polarization vectors and the θ * e + angle are defined in the CM frame [ Figure 6(a)], the relative EW corrections are flat for both polarized and unpolarized signals. On the contrary, for the LAB-frame definition [ Figure 6(d)], the EW corrections enhance the various polarized signals differently. They are rather uniform for the LL, LT and TT components, while for the TL component they increase by about 6% going from the anti-collinear configuration towards the collinear one. Therefore, differences between the two polarization definitions are found at the level of EW corrections, even for very inclusive variables like decay angles. Thanks to the flat behaviour of EW    corrections, the CM-frame definition seems preferable over the LAB-frame one. If polarizations are defined in the CM (LAB) frame, the distributions in the angle θ * e + computed w.r.t. the Z direction in the LAB (CM) frame receive large interference effects, as can be appreciated in Figures 6(b)-6(c). The interference pattern is different in the two cases: interferences are in fact negative in the (anti)collinear configuration in Figure 6(c) but positive in Figure 6(b). A difference in size is also evident: if polarizations are defined in the CM frame and the decay angle is defined w.r.t. the Z direction in the LAB frame [ Figure 6(c)], the interferences give at most a 3% effect, while in the opposite case they give up to 6% effects.
The distribution in the decay angle θ * is highly sensitive to the polarization state of the corresponding Z boson. However, in order to avoid noticeable interference effects, it is essential that the definition of the angle is consistent with the frame where polarizations are defined, as in the case of Figures 6(a) and 6(d).
In Figure 7 we consider the azimuthal distance between the positron and the muon computed w.r.t. to the beam axis. The two particles are decay products of different bosons, and tend to be produced with ∆φ e + µ − = π. For both the unpolarized and the polarized distributions the EW corrections monotonically decrease with increasing ∆φ e + µ − . For unpolarized cross-sections, EW corrections range between −9% and −13%. For polarizations defined in the CM frame, all polarized cross-sections receive the same flat corrections for ∆φ e + µ − < π/2. For ∆φ e + µ − > π/2, the EW corrections for TT decrease along with the unpolarized one, while those for other polarizations remain rather flat. A similar behaviour is found in the LAB frame, where all polarized curves (except from the TT one) are pretty flat in the whole angular range.
The two definitions of polarizations lead to very different interference patterns to this angular variable. In the CM-frame definition, the interferences never exceed 0.8%, with mild differences between the LO and the NLO EW results. In the LAB definition, the NLO EW interferences are positive for ∆φ e + µ − < π/2 and amount to 2.5% at ∆φ e + µ − = 0, they are negative but less sizeable in the rest of the spectrum. In Figure 8 we consider differential cross-sections in the positron transverse momentum. The purely-longitudinal signal is strongly suppressed at large transverse momentum already at LO [100] and is also characterized by huge and negatively-increasing EW corrections for both polarization definitions. Around 280 GeV the negative virtual corrections exceed the sum of Born and real contributions, implying that the truncation of the perturbative series at order O(α 5 ) is not sufficient at large transverse momentum. We have checked numerically that the same effect can be reproduced simulating the production of longitudinal on-shell bosons (no decays), which was also addressed in Ref. [17]. The longitudinal signal is suppressed by 1/s and 1/s 2 w.r.t. the mixed and the TT contributions, respectively, both at LO and in NLO real corrections, which manifests itself by vanishing amplitudes for neutral Goldstone-boson pair production at tree level. On the contrary, some one-loop diagrams give non-vanishing contributions to Goldstone-boson amplitudes, which at NLO are interfered with tree-level amplitudes. Including the squares of such weak one-loop contributions [formally at order O(α 6 )] makes the cross-section positive even at large transverse momentum [17]. However, the smallness of polarized signals at large p T,Z renders this region definitely out of reach even with increased LHC luminosity. Similar effects can be found in the distributions in the lepton-pair transverse momentum, as well as in the four-lepton invariant mass.
The TT polarization state dominates already at moderate transverse momentum: in the CMframe definition, the mixed states are one order of magnitude smaller than the TT one already at p T,e + ≈ 100 GeV, despite less sizeable EW corrections. In the LAB-frame definition, the mixed states and the TT one have similar size in the soft region of the spectrum. It is interesting to  observe the differences in the interference pattern between the two polarization definitions. In the CM-frame definition they amount to ±3% in the soft region, while they are almost negligible in the tail. In the LAB-frame definition, interferences are positive and more sizeable (at most 6%) even at moderate-p T . Although the interference effects account for 0.2% of the integrated cross-section, they can give more sizeable effects at the differential level, even in the case of inclusive variables. This holds in particular when polarizations are defined in the LAB frame. The impact of the non-doublyresonant background on the considered distributions is not so different from the one found for the total cross-section. The NLO EW corrections significantly distort the distribution shapes for both polarized and unpolarized signals, with moderate effects in angular variables and dramatic effects in the tails of transverse-momentum and invariant-mass variables.

Fiducial phase-space: complete NLO results
Thanks to the back-to-back kinematics at LO, the definition of polarizations in the CM frame represents a natural choice for di-boson production. It gives an enhanced LL contribution, rather small mixed contributions, and less sizeable interference effects than the polarization definition in the LAB frame. Therefore, we only present results for the CM-frame definition in the fiducial setup.
The realistic cuts applied in the FID setup are rather loose in the lepton transverse momentum and R distance. The decrease of the inclusive cross-section due to the fiducial cuts is mostly due to the rapidity cuts, which also introduce an asymmetry between the two Z bosons (|y e ± | < 2.47 and |y µ ± | < 2.7). Once the invariant-mass cut on the four-lepton system, M 4 > 180 GeV, is applied, the cuts on the leading and sub-leading lepton do not modify much the event rate. The presence of lepton cuts is expected to give larger interference effects than in the inclusive case [67,68,71].
In Table 3 we show the integrated cross-sections in the FID setup, both for polarized and unpolarized signals. We have included NLO effects of EW and QCD type, as well as the loopinduced contribution from the gg partonic channel. Starting from the LO results, we have combined the higher-order corrections both with an additive and with a multiplicative approach [31]: where we have defined the relative corrections in the following way: The cross-section σ LO includes LO contributions in the qq and γγ partonic channels. The quantity ∆σ EW embeds all EW corrections in the qq, γγ, γq and γq channels, while ∆σ QCD furnishes the QCD corrections in the qq, gq and gq channels. Loop-induced contributions from the gg partonic channel are included in σ gg . The combination of corrections according to Eqs. (3.3) and (3.4) is performed on a bin-by-bin basis. The difference between the multiplicative and the additive results (= δ EW δ QCD dσ LO ) gives an approximated estimate of the mixed EW-QCD corrections of order O(α 5 α s ). The multiplicative approach is generally preferable at large partonic energies (ŝ M 2 Z ) as the dominant EW corrections factorize w.r.t. the QCD effects. In principle, for the factorization properties of EW and QCD corrections, the multiplicative approach should be applied to the qq contributions only. However, we have checked numerically that, due the smallness of photon-induced contributions, this further refinement of the multiplicative prescription agrees at the permille level with results obtained applying Eq. (3.4), both for integrated cross-sections and for differential distributions. This conclusion was also found when combining NNLO QCD and NLO EW corrections to ZZ production [31]. Therefore, in the following we stick to the simple multiplicative prescription of Eq. (3.4).
As in the inclusive setup (see Table 2), the EW corrections are negative and sizeable, accounting for about −10% of the LO results for all polarized signals with at least one longitudinal boson, about −11% for the purely-transverse and unpolarized signals. The QCD corrections give the largest correction to LO cross-sections, ranging between +30% and +45% for the various polarized and unpolarized signals. Such corrections are particularly large for the mixed states, an effect that characterizes also W + Z production [80]. The loop-induced contribution from the gluon-gluon partonic channel gives a sizeable enhancement to the LL and TT signals of 15% and 20%, respectively. For mixed polarization modes, this gluon-initiated contribution is only at the level of 3%.
The additive combination of NLO corrections and loop-induced contributions results in an enhancement of the LO cross-section of about 40% for all signals, while the multiplicative combination gives a smaller overall correction of about 35% for all signals. At the level of total cross-sections the difference between the additive and multiplicative results is about 2.8%, which is smaller than NLO scale uncertainties.
The uncertainties of the fiducial NLO cross-section from 7-point scale variations are strongly affected by the gluon-induced partonic channel. In fact, the noticeable reduction of scale uncertainties from LO to NLO QCD (from order 5% to order 2%) is counterbalanced by 25% scale uncertainties of the gluon-induced contribution, resulting in only a mild reduction of the combined scale uncertainties. The most noticeable scale-uncertainty reduction is found for the mixed polarization states (from 7% at LO to 4% at NLO), owing to the small size of loop-induced contributions.
Thanks to the joint effect of the four-lepton and lepton-pair invariant mass cuts, the non-doublyresonant background accounts for just 0.8% of the full result at LO and even less when including radiative corrections (0.5% for NLO × ). The interferences are larger than in the inclusive setup. The polarization fractions and the impact of the interferences are shown in Table 4, as relative percentages of the DPA unpolarized result. Interference effects are at the level of 1% both at LO and at NLO, but larger effects could be present in differential results. They are particularly small for the gluon-induced process. The LO polarization fractions are roughly conserved when including NLO EW corrections, with the LL contribution accounting for 6% of the total, which is dominated by purely-transverse states (70%). The NLO QCD fractions show a 0.8% increase in the mixed cases, and correspondingly a 3% decrease in the TT one. The polarization fractions for the loop-induced contribution differ greatly from those found in qq channels, with very small mixed contributions of about 2% each and a longitudinal fraction of about 5%, which mostly comes from box diagrams involving top loops, and a 90% TT contribution. The combined NLO fractions are roughly the same as the LO ones, with basically no difference between the additive and multiplicative approaches. We now present differential results in the fiducial region for observables that are relevant for the discrimination among polarization states of the Z bosons using the multiplicative combination of EW and QCD corrections [Eq. (3.4)]. We have checked numerically that the corresponding differential results obtained with the additive prescription [Eq. (3.3)] are very similar. Non-negligible differences between the two approaches only appear in phase-space regions where the cross-section for ZZ production is suppressed. While the multiplicative combination of NLO corrections is theoretically preferable over the additive one, it causes artificial effects in phase-space regions where the LO cross-section is suppressed, particularly for the LL polarization state. In Figures 9-17 we focus on the size of interferences, on the relative impact of higher-order corrections, and on the differences among the shapes of polarized distributions.
As first observable we consider the transverse momentum of the electron-positron system in Figure 9. In the case of DPA calculations, this observable coincides with the transverse momentum of the Z boson that decays into the two dressed leptons e + and e − . In the soft region of the spectrum, the non-doubly-resonant background (the deviation of the grey curves from 1 in the left middle panel) contributes 6% (4%) to the full result at LO (NLO × ). In the same kinematic region, the interferences (the difference between the grey and magenta curves in the left middle panel) range between 2% and 3%, while they are below the 1% level at moderate and large transverse momenta. The unpolarized distribution is dominated by the TT polarization mode both at LO and when including radiative corrections. The LL and LT polarized distributions are strongly suppressed at large transverse momenta, where they are two orders of magnitude smaller than the unpolarized and purely-transverse ones, while the TL mode is only suppressed by one order of magnitude. The LL polarization state, which is strongly suppressed at LO w.r.t. to all other polarization states, receives a sizeable enhancement from QCD real corrections and loop-induced contributions, already at moderate transverse momentum. The two mixed-polarization distributions have similar behaviours for p T,e + e − 100 GeV, while at larger values the TL mode receives huge (up to 1000%) corrections from real QCD radiation, mostly coming from the gq, gq partonic processes that open up at NLO QCD. In the TL polarization configuration, a large p T of the transverse boson implies that also the longitudinal boson has large p T at LO, giving a suppression in the cross-section by 1/p 2 T w.r.t. to the TT polarization state. In the presence of QCD radiation, the transverse momentum needed to balance a large p T of the transverse boson is shared by the longitudinal boson and the radiated gluon or quark. Therefore, the longitudinal Z boson carries smaller p T than at LO giving  Figure 9.
Distributions in the transverse momentum of the e + e − system. Unpolarized and doublypolarized results are shown for ZZ production at the LHC in the fiducial setup described in Section 2.5 with polarizations defined the CM frame. Left panel from top down: NLO× differential cross-sections, ratios of DPA over the full result, normalized NLO× shapes (unit integral). Right panel from top down: NLO× K-factor, NLO EW corrections δEW, NLO QCD correction δQCD, gg contribution δgg. All curves in the right panel are relative to the LO cross-sections. a much less severe suppression to the signal. For the LT polarization configuration, on the other hand, the presence of real radiation does not influence the transverse momentum of the longitudinal boson in the considered distribution. The large transverse momentum of the longitudinal Z boson suppresses the contributions at LO and NLO QCD alike. This argument explains a LT signal that is 10 times smaller than the TL one for p T,e + e − 500 GeV and the enhancement of the LL mode by real QCD radiation. An analogous effect is not present at NLO EW, because the EW virtual corrections are noticeably larger than real ones, the EW coupling is suppressed w.r.t. the QCD case, and the photon-induced partonic processes (γq, γq) are less luminous than those with a gluon in the initial state. The loop-induced gluon-gluon contributions are particularly sizeable for the TT signal in the very soft region of the spectrum, where the quark-induced LO TT cross-section is suppressed owing to angular-momentum conservation (the soft-p T region corresponds to a scattering angle that is close to zero or π). Although all polarized distributions are peaked at small transverse momenta, the transverse momentum of a single Z boson is sensitive to the polarization state of both Z bosons thanks to the strong kinematic and spin correlation between the two bosons. In fact, the normalized shape of the LL distribution is much narrower than the TT one, while the normalized shapes of mixed-polarization configurations show an intermediate behaviour between the LL and TT ones. This indicates that the transverse momentum of a single Z boson provides a handle to separate the LL polarization state from the others.
In Figure 10, we display the distributions in the invariant mass of the e + µ + pair. This variable shares strong similarities with the invariant mass of the four-lepton system (not shown) concerning relative NLO corrections, but is characterized by an enhanced discrimination power among  Figure 10. Distributions in the invariant mass of the e + µ − system. Same structure as in Figure 9.
polarizations. It has been found to be sensitive to polarizations of Z bosons also in Higgs decays [53,54]. The non-doubly-resonant background is as small as at the integrated level, and the interferences amount to at most 2% of the unpolarized result in the considered range. Although the LL and mixed modes are strongly suppressed w.r.t. the TT mode already at moderate invariant mass, in the soft region of the spectrum the shapes of polarized distributions are different. The TT distribution is characterized by a maximum around 50 GeV and a rather large width, while the maximum of the LL curve is around 80 GeV and the shape is narrower. The results for the mixed states lie in between those of the TT and LL ones. The suppression of the LO distributions for the polarization states with one or two longitudinal bosons at high M e + µ + entails a strong impact of NLO and loop-induced contributions. As observed in invariant-mass and transverse-momentum distributions in the inclusive setup [see Figure 8(a)], the EW corrections grow large and negative for the purely-longitudinal signal, accounting for −40% of the LO prediction at M e + µ + ≈ 400 GeV, while milder EW corrections are found for other polarization states, ranging between −10% and −20% in the considered range. The QCD corrections are similar to those at the integrated level in the soft region of the spectrum, while for M e + µ + 100 GeV the LL and mixed signals are enhanced by very large real-radiation contributions both in the qq and in the gq, gq partonic channels. The QCD radiative effects are particularly large for the mixed states, resulting in NLO QCD crosssections that are five times larger than the LO ones at M e + µ + ≈ 600 GeV. This effect is due to huge real contributions from the gluon-(anti)quark channels that enhance the mixed signals more sizeably than the LL one. 2 The purely-longitudinal signal is affected by large gluon-induced contributions that grow fast towards larger invariant masses, where the LO qq contributions are highly suppressed. The multiplicative combination of corrections gives huge and monotonically increasing K-factors for the LL state, slightly smaller but still very large K-factors for the mixed states, while the TT and unpolarized results feature decreasing K-factors (+40% in the soft region, +25% at M e + µ + ≈ 600 GeV) owing to the negative EW corrections that are also applied to the flat QCD  Figure 11. Distributions in the scattering angle of the Z boson decaying into e + e − in the four-lepton CM frame computed w.r.t. the four-lepton-system direction in the LAB frame. Same structure as in Figure 9.
correction, and to the decreasing gluon-induced contribution in the high-energy regime. Very large effects of higher-order corrections to longitudinal signals are not only found in the tails of invariant-mass and transverse-momentum distributions, but also in inclusive angular variables like the cosine of the scattering angle between the two Z bosons presented in Figure 11. The scattering angle θ CM V1 is defined as the angle between the momentum of the e + e − pair computed in the rest frame of the four-lepton system and the momentum of the four-lepton system computed in the LAB frame. At LO it coincides with the scattering angle between the two Z bosons that is used (for on-shell bosons) in Eq. (2.22). At NLO it receives substantial modifications by realradiation contributions. The interferences are almost vanishing for θ CM V1 = π/2 , while they reach the 3% level in the forward and backward regions, where the unpolarized, transverse-transverse and mixed-polarization distributions have their maxima. The peaks in these regions are particularly pronounced for the mixed signals. The LL signal has maxima at slightly more central configurations (cos θ CM V1 ≈ ±0.9). At cos θ CM V1 = 0 the polarized signals display a rather different behaviour. The TT signal has a minimum, the mixed states have a local maximum, while the purely-longitudinal distribution is characterized by an impressive enhancement, which is an artificial effect of the multiplicative combination of the NLO corrections. In fact, the LL signal vanishes at cos θ CM V1 = 0, while the NLO corrections do not, i.e. the additional term w.r.t. the additive combination (∆σ EW ∆σ QCD /σ LO ) becomes huge. We have checked that the additive combination gives a smooth LL curve around cos θ CM V1 = 0. As a consequence, the multiplicative combination of NLO corrections should be applied with a grain of salt and compared against the results of the additive combination, particularly in those regions where the LO is very small or vanishing. The fake effects of the multiplicative combination at cos θ CM V1 = 0 are also present in the K-factor shown in the right top panel of Figure 11, while the huge effects in the relative corrections δ EW , δ QCD , and δ gg simply result from the normalisation to the vanishing LO cross-section. Also EW corrections to the LL distribution become large and positive in the central region, while they are negative and in its corresponding Z-boson rest frame computed w.r.t. the Z-boson spatial direction in the CM frame. Same structure as in Figure 9.
flat towards cos θ CM V1 = ±1. The combined NLO × corrections to the TT distribution are relatively flat, apart from an enhancement in the forward and backward regions due to the gluon-induced contribution. Mixed-polarization distributions receive NLO QCD corrections and gluon-induced contributions that are maximal around cos θ CM V1 = ±0.75. A similar shape in the relative correction is found also for EW corrections. Up to the effects of the multiplicative combination, the different shapes for the various polarized signals and the rather small interference effects make this angular variable very suitable for the discrimination among polarization states.
The possibility to fully reconstruct the rest frame of each Z boson gives access to the lepton decay angles, which are the most direct probes of the polarization states of decayed bosons. We consider in Figure 12 the cosine of the polar decay angle of the positron in the corresponding Zboson rest frame. The NLO EW distributions in the inclusive setup for the same variable are shown in Figure 6(a). The effect of transverse-momentum and rapidity cuts on this observable is mild and restricted to cos θ * ,CM e + = ±1, where the TT and TL distributions do not have a maximum anymore. The size of interferences is similar to the one found at the integrated level, ranging from 0.3% of the full result in the central region to 2% around cos θ * ,CM e + ≈ ±0.7. This observable is highly sensitive to the polarization of the Z boson that decays to e + e − , as can be appreciated from the comparison of normalized shapes for the various polarized signals. There is no shape difference between the LL and LT curves, despite the different cuts applied on the e ± and µ ± kinematics, and very small shape differences are found between the TL and TT curves. Obviously, in order to discriminate between the longitudinal and transverse polarizations of the Z boson decaying into µ + µ − the cos θ * ,CM e + variable is not helpful. The distributions are not distorted sizeably by the inclusion of higher-order corrections neither for polarized signals nor for the unpolarized ones. In fact, the NLO effects and the loop-induced contribution entail rather flat corrections to this observable, resulting in combined K-factors that are comparable with those of the integrated results in the whole angular  Figure 13. Distributions in the difference between azimuthal decay angles φ * ,CM e + and φ * ,CM µ + , where φ * ,CM + is the decay angle of the positively-charged lepton in its corresponding Z-boson rest frame computed w.r.t. the Z-boson spatial direction in the CM frame. Same structure as in Figure 9.
range. Owing to lepton universality, the results for the cos θ * ,CM µ + variable are very similar to those for cos θ * ,CM e + up to small effects due to the different transverse-momentum and rapidity cuts for the two lepton flavours.
It is clear that any experimental analysis that targets the polarization of Z bosons should rely on decay angles as discriminating variables. In the rest frame of a single Z boson, the leptonic-decay kinematics is fully described by the polar decay angle, θ * ,CM + , and the azimuthal one, φ * ,CM + . If two Z bosons are produced, it is interesting to investigate the spin correlation between them by means of the difference between azimuthal decay angles of the positively-charged leptons, each computed in the corresponding Z-boson rest frame. This variable has been proved to be sensitive to the polarization of weak bosons that come from a Higgs decay [53,54]. Since polarizations are defined in the CM frame, the reference direction for the angle φ * ,CM e + coincides with the one for φ * ,CM µ + , up to a minus sign, and both angles are measured relative to the same plane. The azimuthal difference is defined as ∆φ * ,CM e + µ + = min φ * ,CM and the distributions in this variable are shown in Figure 13. The polarized and unpolarized distributions are symmetric about ∆φ * ,CM e + µ + = π/2, and the non-doubly-resonant background reflects the integrated results with no shape distortion. The interferences vary between 0.5% and 3% at LO, showing a sizeable increase towards ∆φ * ,CM e + µ + = π. On the contrary, the combined NLO results give rather constant interference effects over the whole angular range. The QCD corrections are pretty flat for all polarized and unpolarized signals, while a noticeable shape difference is found between LO and loop-induced results for the TT signal. The EW radiative corrections are constant for polarized signals, while a 1.5% variation shows up in the relative correction for the unpolarized case, which is due to the different interference patterns at LO and at NLO EW. While an 8% variation is found in the combined K-factor for the TT signal, the variations for signals with one or  Figure 14.
Distributions in the azimuthal separation between the positron and the electron. Same structure as in Figure 9.
two longitudinal bosons are much smaller. The comparison of normalized shapes demonstrates that this variable is particularly suitable for the separation of the TT polarization mode. In fact, the distribution for this mode has a maximum at ∆φ * ,CM e + µ + = π/2 and minima at ∆φ * ,CM e + µ + = 0, π, while the other modes exhibit the opposite behaviour. Furthermore, the LL and mixed distributions are somewhat flatter than the TT one. Separating the purely-longitudinal state from the mixed ones is less straightforward with this observable, as the shape differences among these modes are mild.
Also other angular variables are sensitive to the polarization of Z bosons. In Figure 14 we present the differential distributions in the azimuthal separation (computed in the LAB frame) between the positron and the electron. In the DPA calculations, either polarized or unpolarized, the two leptons are decay products of a Z boson and are typically produced in opposite hemispheres. The non-doubly-resonant effects are compatible with zero in most of the spectrum of this observable, while they get slightly larger in the most populated region, reaching 2.5% and 2% at LO and at NLO × , respectively. Around ∆φ e + e − = π also large interferences show up, contributing up to 5% to the full results, both in LO and in NLO × predictions. The NLO EW corrections are at the −10% level in the most populated region, while they become more sizeable at ∆φ e + e − = 0 (−16% for diagonal polarization states, −12% for mixed ones). Huge QCD corrections (more than 300%) are found for the TL signal in the regime of small azimuthal separation, which result from real corrections, in particular from those with initial-state gluons. The QCD corrections are smaller, but still large, for the LL signal (around 100%), while for other polarized and for the unpolarized signals QCD corrections never exceed 60% in the whole range. Large gluon-induced contributions to the purely-longitudinal signal are found for ∆φ e + e − < π/2. The combination of all corrections gives large K-factors for the TL signal (up to 4.5) and to the LL one (up to 2.7). The comparison of polarized signals reveals different distribution shapes, with monotonically-increasing curves for the mixed states, a plateau in the region ∆φ e + e − > 5π/6 for the TT state, and a maximum near ∆φ e + e − = 8π/9 for the LL state. This renders this observable suitable for the discrimination  Figure 15. Distributions in the azimuthal separation between the positron and the muon. Same structure as in Figure 9.
among polarization modes in LHC data, provided that also the interference contribution and the irreducible background are properly accounted for in a SM-template fit of the data. Another variable that is suitable for the polarization discrimination is the azimuthal separation between the positron and the muon considered in Figure 15. Since the two particles originate from different Z-boson decays, the distributions are less peaked at π than the ∆φ e + e − ones. Both LO and combined NLO results give positive interferences ranging between 1% and 2% of the full distribution in the whole angular range. The NLO EW corrections cause sizeable distortions only to the TT polarization mode, as already observed in the inclusive setup [see Figure 7(a)]. The impact of QCD corrections diminishes monotonically towards the most populated region, with a 10% variation for the LL and mixed signals and a slightly smaller for the TT one. The gluon-induced process enhances the LO distributions in a rather uniform way for most polarization modes, apart from a 10% increase for the LL mode. The combined NLO distributions suggest that this angular variable is suitable for the separation of the TT signal from the others. In fact, at the level of normalized shapes, the TT curve increases by 100% going from ∆φ e + µ − = 0 to ∆φ e + µ − = π, while other polarized curves increase by roughly 50% only.
In Figure 16 we present the distributions in the rapidity of the electron-positron system. Owing to the selection cut on the rapidities of both e + and e − , |y e ± | < 2.47, also the combination of the two momenta, i.e. the reconstructed Z boson, is central. The non-doubly-resonant effects are flat and equal to those in the total cross-section. A positive 3.5% interference contribution appears in the most forward and backward configurations allowed by the selections, while in the central region, which is the most populated one, the interferences account for less than 1%. The EW corrections are pretty flat for all polarization modes, while some distortions of the LO curves are introduced by QCD corrections and by loop-induced contributions. The NLO QCD corrections enhance in particular the TL distribution by up to 70% in the most forward and backward regions, while milder effects are found for other polarization states. The loop-induced contributions sizeably  Figure 16. Distributions in the rapidity of the e + e − pair (single Z boson). Same structure as in Figure 9.
enhance the LL and the TT distributions, with maximal effects in the central region (23% for TT, 15% for LL, at y e + e − = 0). The multiplicative combination of corrections results in K-factors varying between 1.3 and 1.4 (with a similar pattern) for all polarized signals but the TL one, which varies from 1.3 in the central region to 1.55 in the most forward and backward regions. All combined distributions have a maximum at zero rapidity, and are characterized by a negative convexity in the whole available range, with very small shape differences among different polarization states. This shows that the rapidity of a single Z boson is hardly sensitive to the polarization of the decayed boson, as it has been already found in WZ production with polarizations defined in the di-boson CM frame [80]. Finally, we consider distributions in the absolute value of the rapidity separation between the two Z bosons in Figure 17. The momenta of the Z bosons are determined from the momenta of their decay products after possible photon recombination. Owing to the rapidity cuts on all charged leptons, the rapidity distance between the two opposite-flavour lepton pairs is implicitly cut such that |∆y ZZ | 5, but already at |∆y ZZ | ≈ 4 all cross-sections are more than 3 orders of magnitude lower than those at the maximum of the distributions. The region of small rapidity separation is the most populated one for all polarized and unpolarized signals. In this configuration, the events are characterized by slightly larger non-doubly-resonant effects (2%) than those found at the integrated level, and by almost vanishing interferences, both at LO and at NLO. For larger rapidity separations, the non-doubly-resonant background becomes smaller than 0.5% of the full result, while interferences increase up to 3%. The NLO EW corrections, despite the different patterns for various polarization modes, induce only mild shape distortions of the distributions with at most 3% variations in the considered range for the LL and mixed polarization states. Even flatter EW corrections are found for the purely-transverse and for the unpolarized distributions. At variance with the EW effects, the NLO QCD corrections sizeably increase from small to large rapidity separations for the LL and mixed polarization states reaching 120% and 180%, respectively, at |∆y ZZ | = 3. A 50% enhancement results for the LL state at zero rapidity separation, due to the fact  Figure 17. Distributions in the absolute value of the rapidity separation between the two Z bosons. Same structure as in Figure 9.
that this polarization state is strongly suppressed in such configuration at LO. The QCD corrections for the purely-transverse state turn out to be relatively flat. The loop-induced process enhances the suppressed LO cross-section with both longitudinally-polarized bosons by 150% at |∆y ZZ | ≈ 0 and causes effects of up to 50% for large rapidity separations for the same polarization mode. About 20% loop-induced contributions characterize the TT signal at small rapidity separation. The normalized TT distribution has a flat maximum at |∆y ZZ | ≈ 0.2, while mixed contributions have pronounced maxima at zero rapidity separation. The LL distribution is maximal at |∆y ZZ | ≈ 0.5 with a shape that is noticeably different from the other three polarization states. Therefore, using this observable in a fit to LHC data is expected to improve the capability of isolating the purely-longitudinal signal from the other ones.

Conclusion
Extending the methods of Refs. [79,80] to NLO EW corrections and in particular to the emission of photons off the decay products of Z bosons, we have presented a general procedure to compute crosssections for the production of Z bosons in definite polarization states valid up to NLO accuracy both in the EW and in the QCD coupling. The definition of polarized signals relies on the application of a double-pole approximation (DPA) to contributions with two resonant Z bosons and on the separation of polarization states at the amplitude level. While the method is applicable to any LHC process with neutral resonances, even including identical particles in the final state, we have focused on inclusive Z-pair production in the different-flavour four-charged-lepton channel at complete NLO accuracy, including also gluon-initiated, loop-induced contributions. We have first applied the newly-introduced DPA techniques to the unpolarized process at NLO EW for very inclusive selection cuts, finding a small contribution (≈ 1%) of the non-doubly-resonant background to the full computation. Beyond its phenomenological importance, this result provides a validation of the DPA method for complete NLO corrections, and in particular of the well-behaved definition of subtraction counterterms that are needed for the cancellation of IR singularities in the DPA calculation.
In the same inclusive setup, the separation of polarizations for both Z bosons has been investigated, comparing the results of two different definitions of polarization vectors, one in the di-boson centre-of-mass (CM) frame and one in the laboratory (LAB) frame. Despite the absence of cuts on single leptons, non-negligible interference patterns are found in some differential distributions both at LO and at NLO EW. Sizeable differences show up between the CM-and LAB-frame definitions of polarizations, but smaller interference effects are present when defining the polarizations in the CM frame, which is also more natural for ZZ production.
In the presence of realistic fiducial cuts, inspired by the latest ATLAS analysis, we have combined NLO EW and QCD corrections and gluon-induced contributions (formally of NNLO QCD accuracy) both in a multiplicative and in an additive way. While the multiplicative combination is better motivated from the theory point of view in general, it introduces artificial effects in regions of phase space where the purely-longitudinal signal is suppressed at LO.
Scale uncertainties are at the 5% level in the combined NLO predictions both for unpolarized and polarized signals. A substantial fraction of the NLO scale uncertainties is due to the gluon-induced process, in particular for the polarization modes with equally polarized Z bosons (transverse-transverse and longitudinal-longitudinal).
In spite of sizeable and negative EW corrections of about −10%, the largest radiative effects originate from NLO QCD corrections which cause an increase of the LO cross-section of about 35%. Also the gluon-induced process contributes noticeably to ZZ production (≈ 15-20% for unpolarized ZZ pairs and those with equal polarizations). Polarization fractions are roughly conserved going from LO to NLO accuracy, but sizeable distortions are found at the differential level with rather different K-factors for the various polarization states.
At variance with di-boson processes with neutrinos as decay products, the possibility to fully reconstruct the final state allows for access to a large set of observables. We have presented a number of differential results for fiducial event selections, focusing on those variables that are particularly suited for the separation of polarized signals. The non-doubly-resonant background (0.5% at the integrated level) never exceeds 3% of the full result for all presented differential distributions. Interferences, which are at the 1% level in the integrated cross-section, can reach up to 5% in more exclusive phase-space regions. Several angular variables, including decay angles as well as azimuthal and rapidity separations are well suited for polarization discrimination. However, also some invariant-mass distributions show a fair discrimination power among different polarization states. A combination of the variables considered in this paper, either in a multivariate analysis or as input to machine-learning techniques, is expected to allow for the separation of polarized states in the LHC data despite the small rate of the unpolarized process and the strongly-suppressed cross-section for longitudinal bosons.