Evidence of pair production of longitudinally polarised vector bosons and study of CP properties in 𝒁 𝒁 → 4 ℓ events with the ATLAS detector at √

A study of the polarisation and CP properties in 𝑍𝑍 production is presented. The used data set corresponds to an integrated luminosity of 140 fb − 1 of proton–proton collisions at a centre-of-mass energy of 13 TeV recorded by the ATLAS detector at the Large Hadron Collider. The 𝑍𝑍 candidate events are reconstructed using two same-flavour opposite-charge electron or muon pairs. The production of two longitudinally polarised 𝑍 bosons is measured with a significance of 4.3 standard deviations, and its cross-section is measured in a fiducial phase space to be 2 . 45 ± 0 . 60 fb, consistent with the next-to-leading-order Standard Model prediction. The inclusive differential cross-section as a function of a CP-sensitive angular observable is also measured. The results are used to constrain anomalous CP-odd neutral triple gauge couplings.


Introduction
Measurements of  boson pair ( ) production in proton-proton ( ) collisions at the Large Hadron Collider (LHC) provide an important test of the Standard Model (SM) gauge structure in the electroweak (EW) sector.Representative leading-order Feynman diagrams of   production at the LHC are shown in Figure 1.With increased luminosity, experimental data enables studies beyond precision measurements of integrated and differential cross-sections, such as the weak boson polarisation and charge conjugation (C) and parity (P) properties.The polarisation measurement of massive weak bosons is a direct probe of the Electroweak Symmetry Breaking mechanism, through which the  and  bosons obtain their longitudinally polarised states.Diboson polarisation measurements, especially those probing longitudinally polarised vector bosons, provide unique sensitivity to physics beyond the SM (BSM) [1].
Polarisation measurements performed with LHC data have mainly focused on single bosons such as in -boson production [2, 3], -boson production [4,5], and  bosons from top-quark decays [6][7][8].Single-boson polarisation states have also been measured in   production by the ATLAS [9] and CMS [10] Collaborations.Recently first measurements on joint-polarisation states of weak bosons have also been reported.The production of two transversely polarised  bosons has been measured by the CMS Collaboration in EW production of the same-sign  ±  ± boson pairs [11], while the the ATLAS Collaboration has measured the joint-polarisation states of longitudinally or transversely polarised bosons in inclusive   production [12].In the latter, the fraction of diboson events with a simultaneous longitudinal polarisation (LL) was observed with a significance of 7.1 standard deviations.This paper presents a measurement of the production of two longitudinally polarised  bosons ( L  L ) in the decay channel   → ℓ + ℓ − ℓ ′+ ℓ ′− , where ℓ and ℓ ′ can be an electron or a muon.The -boson candidates are reconstructed with same-flavour, opposite-charge (SFOC) electron or muon pairs, and they are required to be on-shell with | ℓℓ −   | < 10 GeV, where  ℓℓ is the invariant mass of the lepton pair and   is the -boson pole mass [13].
The violation of CP symmetries is required to explain the matter-antimatter asymmetry in the universe, and it is well known that there is insufficient CP violation in the SM [14 -16].The measurement of CP-sensitive observables in diboson production can be utilised to explore new sources of CP violation in the gauge-boson sector.CP-violating effects in weak-boson self-interactions were studied in various measurements of diboson production at the LHC by constraining the CP-odd anomalous neutral triple gauge couplings (aNTGC), including those entering the    and   vertexes, using   production in ATLAS and CMS [17][18][19][20][21][22][23].Such experimental searches primarily derive constraints on anomalous triple gauge boson couplings (aTGC) using event rates or cross-section measurements without employing dedicated CP-sensitive observables.This paper presents the differential cross-section for a dedicated CP-odd angular observable, referred to as the Optimal Observable (OO).The OO is defined in Section 6 using the decay products of weak bosons in   production, in such a way as to be sensitive to BSM amplitudes through the interference to the SM [24,25].The results are then reinterpreted to constrain aNTGC using an effective vertex function approach [26].The ATLAS Collaboration has previously used such type of dedicated CP-sensitive observables in the EW    production to test CP violation in the weak-boson self-interactions [27].
The CP property is studied using an aNTGC vertex that can be parameterised with two coupling parameters  4  and  4  that violate the CP symmetry.By using such parameters, the cross-section in any given bin of the CP-sensitive observable can be parameterised as where the superscript  is the bin index of the CP-sensitive observable,  is the CP-odd aNTGC,   SM is the prediction from the SM,   interference is the linear interference between the SM and the aNTGC, and   quadratic is the quadratic contribution of the aNTGC.As pointed out in Ref. [28], for the aNTGC, the quadratic term dominates over the linear interference term.Existing constraints on a CP-odd aNTGC stem primarily from their effect on the cross-section in the high- T regime, derived using kinematic observables such as the leading   T .Such high- T sensitive kinematic observables cannot distinguish CP-even and CP-odd effects.This paper presents a search for CP violation using the unfolded OO which is sensitive to the interference terms.

ATLAS detector
The ATLAS experiment [29] at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and a near 4 coverage in solid angle. 1 It consists of an inner tracking detector (ID) surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic and hadron calorimeters, and a muon spectrometer (MS).The inner tracking detector covers the pseudorapidity range || < 2.5.It consists of silicon pixel, silicon microstrip, and transition radiation tracking detectors.Lead/liquid-argon (LAr) sampling calorimeters provide electromagnetic (EM) energy measurements with high granularity.A steel/scintillator-tile hadron calorimeter covers the central pseudorapidity range (|| < 1.7).The endcap and forward regions are instrumented with LAr calorimeters for both the EM and hadronic energy measurements up to || = 4.9.The muon spectrometer surrounds the calorimeters and is based on three large superconducting air-core toroidal magnets with eight coils each.The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector.The muon spectrometer includes a system of precision tracking chambers and fast detectors for triggering.A two-level trigger system is used to select events.The first-level trigger is implemented in hardware and uses a subset of the detector information to accept events at a rate below 100 kHz.This is followed by a software-based trigger that reduces the accepted event rate to 1 kHz on average depending on the data-taking conditions.An extensive software suite [30] is used in data simulation, in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.

Data and simulation
The data collected by the ATLAS experiment during the 2015-2018 data-taking period of the LHC at a centre-of-mass energy of √  = 13 TeV is analysed, and corresponds to an integrated luminosity of 140 fb −1 .The data events are triggered using single-lepton and dilepton triggers with the minimum trigger threshold, depending on data-taking periods, varying between 20 − 26 GeV for single-lepton triggers and 8 − 24 GeV for dilepton triggers [31,32].The trigger efficiency is nearly 100% for signal events after offline event selections.
Simulations of the signal SM processes of the on-shell production of two  bosons and background processes resulting in two SFOC lepton pairs are derived using Monte Carlo (MC) generators.The  q →   (Figure 1(a)) process was modelled using the Sherpa 2.2.2 generator [33].The matrix elements (ME) were calculated at next-to-leading-order (NLO) accuracy in QCD for up to one additional parton emission and leading-order (LO) accuracy for up to three additional parton emissions.The 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the -axis along the beam pipe.The -axis points from the IP to the centre of the LHC ring, and the -axis points upwards.Cylindrical coordinates (, ) are used in the transverse plane,  being the azimuthal angle around the -axis.
The pseudorapidity is defined in terms of the polar angle  as  = − ln tan(/2).Angular distance is measured in units of Δ ≡ √︁ (Δ) 2 + (Δ) 2 .matrix element calculations were matched and merged with the Sherpa parton shower based on the Catani-Seymour dipole factorisation [34,35], using the MEPS@NLO prescription [36][37][38].The virtual QCD corrections were provided by the OpenLoops library [39][40][41].The NNPDF3.0nnlo set of parton distribution functions (PDF) [42] and the dedicated set of tuned parton-shower parameters developed by the Sherpa authors were used for generating this sample.The higher-order corrections from NLO electroweak effects were taken into account for the  q →   sample by applying reweighting corrections as a function of the four-lepton invariant mass  4ℓ [43,44].
An alternative prediction of the  q →   process with Powheg Box v2 [45][46][47] is used for estimating the theory modelling uncertainty.This sample was generated at NLO accuracy in QCD and interfaced with Pythia 8.186 [48] for modelling the parton shower, hadronisation, and effect of underlying events using the AZNLO set of tuned parameters ('tune') [49].For the hard-scattering and parton showering, CT10 PDF [50] and CTEQ6L1 PDF [51] sets were used, respectively.A higher-order correction as a function of  4ℓ was obtained using a Matrix NNLO QCD prediction [52][53][54][55] and applied to these events.The correction was defined as the ratio of the cross-section at NNLO QCD accuracy to the one at NLO QCD accuracy.
The loop-induced  →   process (Figure 1(b)) was also modelled with the Sherpa 2.2.2 generator using ME calculated with LO accuracy for up to one additional parton emission.As detailed in Ref. [56], the higher-order QCD effects are accounted for by normalising the LO prediction to NNLO using NLO [57,58] and NNLO [59,60] predictions.The same parton showering, matching, and merging schemes as for the Sherpa  q →   simulations were used.The EW production of   in association with two jets,  →     (Figure 1(c)), was modelled by MadGraph5_aMC@NLO 2.6.7 [61].The matrix element was calculated at LO in QCD, and the NNPDF3.0nnloPDF set was used.Pythia 8.244 [62] was used to simulate parton showering, hadronisation, and underlying-event activity using the NNPDF2.3lo[63] PDF set and the A14 tune [64].The three processes,  q →  ,  →  , and  →    , collectively simulate the total signal events for this measurement.
In addition to the simulations for the inclusive   production, samples for different polarisation states were also produced.The polarised  q →   signal samples were simulated using Mad-Graph5_aMC@NLO 2.7.3 [61,65], with the ME calculated at LO in perturbative QCD and with the NNPDF3.0nlo[42] PDF set.The events were interfaced to Pythia 8.240 to model the parton shower, hadronisation, and underlying event, with parameters set according to the A14 tune and using the NNPDF2.3loset of PDFs.Polarised samples were simulated corresponding to the three helicity states,  T  T for two transversely polarised  bosons,  T  L for one transversely polarised  boson and one longitudinally polarised  boson, and  L  L for two longitudinally polarised  bosons, respectively, with each  boson decaying independently into an electron or a muon pair.The -boson polarisation is defined in the centre-of-mass (CM) frame of the two  bosons, which is a natural choice for diboson production [66].To account for the real part of the NLO corrections, the QCD-induced events were simulated with up to two jets in the matrix element at LO and merged with Pythia 8 parton showers using the CKKW-L scheme [67,68].The EW production of the polarised   samples was similarly generated with MadGraph5_aMC@NLO 2.7.3, and they are used together with the QCD-induced  q →   samples in the analysis.Generation of polarised MC events for the loop-induced  →   process is not possible, and a MC event reweighting procedure is applied for the modelling as described in detail in Section 6.1.
The off-shell / * contribution is included in the inclusive   samples but not in the polarised simulations, and the impact is found to be negligible after the event selection with a requirement of two on-shell  bosons, as described in Section 4.3.
Background events consist of four-lepton events originating from  t and   processes, as well as four-lepton events containing at least one non-prompt lepton.The  t events were modelled by the Sherpa 2.2.0 generator at LO accuracy with up to one additional parton emission using the MEPS@LO setup [37,38] and the NNPDF3.0nnloPDF set.The sample was scaled to reproduce the prior  t cross-section measurement from ATLAS [69].Events from the   processes were simulated by the Sherpa 2.2.2 generator with NLO accuracy in QCD for the inclusive process and LO accuracy in QCD for up to two additional parton emissions.The NNPDF3.0nnloPDF set and the default Sherpa parton showering scheme were used.The prediction was scaled to match the prior measurement of the triboson production from ATLAS [70].
The events with any of the four leptons originating from non-prompt sources are estimated by using a semi-data-driven method discussed in Section 5.The method uses information about the origin of the non-prompt leptons from simulations of the relevant processes.The primary source of this background is   production with the leptonic decays of vector bosons and was modelled with Sherpa 2.2.2.The same setup and parameters as with the  q →   modelling were used.The events with non-prompt leptons arising from  + jets processes were modelled using Sherpa 2.2.1 at NLO accuracy in QCD for up to two parton emission and LO accuracy for up to four parton emission.The ME were calculated using the Comix and OpenLoops libraries [34].Matching and merging were performed with the Sherpa parton shower scheme using the MEPS@NLO prescription.The NNPDF3.0nnloPDF set was used, and the events are normalised to a prediction accurate to NNLO in QCD [71] for the inclusive production.The non-prompt events originating from the  t processes were modelled using the Powheg Box v2 generator at NLO accuracy and the NNPDF3.0nloPDF set.The prediction was interfaced with Pythia 8.230 for modelling parton showering, hadronisation, and underlying event with parameters set according to the A14 tune.
The MC samples used in the BSM aNTGC interpretation of the CP study were generated at LO in QCD with MadGraph5_aMC@NLO 2.7.2 using a UFO model [72] as implemented in Ref. [28], and interfaced with Pythia 8.244 for the modelling of parton shower, hadronisation, and underlying event.The aNTGC predictions are reweighted by applying a per-bin -factor, derived by comparing the Sherpa prediction to the LO MadGraph5_aMC@NLO  q →   prediction in each bin of the OO to account for the missing higher-order effects in the BSM prediction.
The predictions from the MC simulations were passed through a detailed simulation of the ATLAS detector [73] based on Geant4 [74].The effect of multiple   interactions in the same bunch crossing, known as pile-up, was emulated by overlaying inelastic   collisions, simulated with Pythia 8.186 using the NNPDF2.3loset of PDFs and the A3 tune [75].The events were then reweighted to match the distribution of the average number of interactions per bunch crossing observed in the data for different data-taking periods.The simulated events were reconstructed using the same algorithms used for the data.Additionally, efficiencies for the trigger, lepton reconstruction, identification, and isolation in MC events were corrected to match those measured in the data.

Fiducial region, object and event selections 4.1 Fiducial region definition
The fiducial phase space is defined close to the detector's kinematic acceptance using particle-level prompt2 leptons produced by a MC event generator, without simulating the effects of the detector.The particle-level prompt leptons are dressed by adding the four-momenta of nearby prompt photons within a small cone of Δ < 0.1.The dressed electrons (muons) are required to be within the detector's acceptance such that they satisfy  T > 7 (5) GeV and || < 2.47 (2.7).The fiducial phase space of the analysis does not include the negligible contribution of prompt leptons from -lepton decays.Events must have a minimum of four prompt dressed leptons to be grouped into at least two SFOC lepton pairs, and the leading and sub-leading leptons must have  T > 20 GeV.The angular separation between any two leptons is required to satisfy Δ > 0.05 to reduce the double counting of detector signatures while keeping leptons from possible boosted production scenarios.The invariant mass of any SFOC lepton pair is required to satisfy  ℓℓ > 5 GeV.
The two SFOC lepton pairs whose invariant masses are closest and next closest to   are designated as an event quadruplet.In the on-shell   region, the resolution on   is comparable to the mass difference between the two  bosons, resulting in some events with inconsistent definitions of the leading and sub-leading pairs at particle and detector levels.This increases the resolution-induced bin migrations that need to be corrected for by the unfolding procedure.Therefore, once the quadruplet is formed, the leading (sub-leading)-lepton pair  1 ( 2 ) is identified as the one with the larger (smaller) value of absolute rapidity, i.e., | ℓℓ |.Based on these requirements, the events are divided into three categories, 4 events with two  +  − pairs, 4 events with two  +  − pairs and 22 events where one of the pairs is  +  − and the other is  +  − .The invariant mass of each SFOC lepton pair is required to be within | ℓℓ −   | < 10 GeV, and the invariant mass of the four leptons is required to be  4ℓ > 180 GeV, motivated to select only on-shell   events.

Object selection
Events are required to contain at least one reconstructed   collision vertex candidate with at least two associated ID tracks with  T > 0.5 GeV.The vertex with the largest sum of  2 T of tracks is considered to be the primary interaction vertex.
Electrons are reconstructed by matching the topological energy clusters deposited in the electromagnetic calorimeters to the tracks in the ID [76].The electron identification is based on a multivariate-likelihood technique that takes information about clusters' shower shapes in the electromagnetic calorimeters, ID track properties, and the quality of track-cluster matching.Baseline electrons, used for the non-prompt background estimate, are required to satisfy  T > 7 GeV, || < 2.47 and the 'VeryLoose' identification criteria [76] and loose association with the primary hard-scatter vertex by requiring | 0 sin | < 0.5 mm, where  0 is the longitudinal impact parameter and  is the polar angle of the track.Signal electrons that define signal events are required to satisfy all of the baseline electron criteria and the stricter 'LooseBLayer' identification criteria and transverse impact parameter significance of | 0 |/  0 < 5, where  0 is the transverse impact parameter relative to the beamline, and   0 is its uncertainty.Prompt leptons originating from hard scattering are characterised by low activity around them in the  −  plane.Therefore, a signal electron must be isolated from other particles by applying criteria for the  T -dependent isolation variable, which is defined using the electron's track and calorimeter energy deposits [76].
Muons are identified using information from various parts of the detector, the ID, the MS, and the calorimeters.A 'Loose' identification working point [77] is used.ID tracks identified as muons based on their calorimetric energy deposits or the presence of individual muon segments are included in the region || < 0.1, and the stand-alone MS tracks are added in the region 2.5 < || < 2.7.Baseline muons used in the non-prompt background estimate are required to satisfy the 'Loose' identification criteria,  T > 5 GeV, || < 2.7 and loose track-to-vertex association of | 0 sin | < 0.5 mm.Signal muons are required to satisfy all criteria for baseline muons and transverse impact parameter significance of | 0 |/  0 < 3. Similarly to the signal electrons, signal muons are required to satisfy additional isolation criteria on the  T -dependent isolation variable, defined using the tracks and particle-flow objects used in muon reconstruction [77].
Jets are reconstructed from particle-flow objects [78] using the anti-  algorithm with a radius parameter of  = 0.4 [79,80].The jet-energy scale is calibrated using simulation and further corrected with in situ methods [81].Candidate jets are required to satisfy  T > 15 GeV and || < 4.5, and are used in the non-prompt background estimation studies.A jet-vertex tagger [82] is applied to jets with  T < 60 GeV and || < 2.4 to suppress jets that originate from pile-up.In addition, -jets are identified using a multivariate -tagging algorithm [83].The chosen -tagging algorithm has an efficiency of 85% for -jets and a rejection factor of 33 against light-flavour jets, measured in  t events [84].
An overlap removal is applied to avoid double counting of the detector signal, favouring leptons with higher  T and preferring non-calorimeter tagged muons over electrons if they share an ID track.The overlap removal rejects any jets within Δ < 0.2 of an electron or jets associated with less than three ID tracks if they overlap with a muon.

Event selection
The reconstructed events must have at least four baseline leptons, and the leading and sub-leading leptons must satisfy  T > 20 GeV.In each event, all possible SFOC lepton pairs are formed by requiring  ℓℓ > 5 GeV and Δ ℓℓ > 0.05 to suppress contributions from the leptonic decays of resonance hadrons.Like the fiducial region, a quadruplet is formed from two SFOC pairs with the lowest values of | ℓℓ −   |.The SFOC lepton pair with the largest value of | ℓℓ | is defined as the leading  boson candidate.To ensure the selection of on-shell   events, the invariant masses of both SFOC lepton pairs are required to be | ℓℓ −   | < 10 GeV and the invariant mass of the quadruplet is required to satisfy  4ℓ > 180 GeV.Each lepton of the quadruplet is required to satisfy the signal lepton definition for the signal region selection.Events with either one or more lepton in a quadruplet failing to meet the signal lepton requirement are used in non-prompt background estimate.Due to the on-shell requirement on both  bosons, in 4 and 4 final states there are less than 2% of events with mispaired leptons, which has no significant impact for both the CP and polarisation measurements.

Background estimation
After the event selection, the background consists of events with one or more of the reconstructed leptons in the quadruplet not originating from a  boson decay.The background processes with prompt leptons from  t and fully leptonic decays of triboson processes are estimated using MC simulations.The measurement also accounts for an additional source of background from non-prompt leptons originating from hadron decays or misidentification of jets.A data-driven fake-factor method is used to estimate the non-prompt background.
A fake-factor quantity, defined as the ratio of signal leptons to the number of baseline leptons failing to meet the signal lepton criteria, is measured from data in dedicated control regions (CR) enriched with non-prompt leptons.The CR where the fake factors are evaluated consists of events with two prompt leptons from a physics process such as  + jets or  t and additional leptons from other sources such as jets.The same triggers and object-level kinematic selections as in the signal region are also applied to select the CR events.The  + jets CR is selected by requiring an SFOC signal lepton pair from a -boson decay with an invariant mass of 76 <  ℓℓ < 106 GeV and additional baseline leptons.Similar to the signal region, SFOC lepton pairs are formed with leptons having Δ > 0.05 and  ℓℓ > 5 GeV.Additionally, events are required to have a missing transverse energy less than 50 GeV to suppress the contamination from the   process.The  t CR is defined by requiring an SFOC signal lepton pair, at least one -tagged jet and at least one additional baseline lepton.
The additional leptons in the  + jets and  t CRs originate predominantly from non-prompt sources such as heavy-or light-flavour jets and are used to derive the fake factors.The additional non-prompt leptons in  + jets events arise dominantly from light-flavour jets, whereas they arise dominantly from heavy-flavour jets in  t events.To match the heavy-flavour composition in the signal region, the two types of events are first weighted and combined.The combination weight is evaluated by comparing the fraction of the non-prompt leptons from heavy-flavour decays in the  + jets and  t CRs to that of the signal region using the simulated samples discussed in Section 3. Since the compositions for the non-prompt electrons and muons are different in the signal region, the weights are evaluated separately to combine the control region events that have additional baseline electrons and muons.When evaluating the fake factors, an estimate of the genuine baseline prompt leptons that fail the signal requirements from the   MC simulation is subtracted for both the  + jets and  t events.The fake factor is measured as a function of  T and  of the non-prompt leptons and the number of jets in an event.
A fake-lepton enriched region is selected with a minimum of four baseline leptons passing the same kinematic selection as in the signal region, but with at least one baseline-not-signal lepton, which is defined as leptons passing the baseline selection but failing the signal-lepton requirement.The non-prompt event yield in the signal region is then estimated by applying a fake-factor weight to each baseline-not-signal lepton.Four-lepton events containing prompt baseline-not-signal leptons are removed from the estimate using simulations.
The background estimate is validated by comparing the prediction of the fake-factor method to the data in two dedicated validation regions where the non-prompt background events are expected to be dominant.The first one is the different-flavour validation region which has the same event selection requirements as the signal region but requires the leptons forming one of the pairs to have different flavours.The second one is the same-charge validation region, which requires the leptons in one of the pairs to have the same charge.In both of the validation regions, the data agree closely with the sum of the estimated non-prompt background yield and the MC prediction within the statistical uncertainties of the data.
Three sources of uncertainties in the background yield are considered.The first uncertainty is related to the statistical uncertainty of the non-prompt leptons in the control region used to calculate the fake factors and ranges from 10 − 80%.The second set of uncertainty, ranging from 2 − 40%, is related to the theory uncertainties of the subtracted prompt-lepton component in the control region, which are dominated by the QCD scale variations.The third and dominant uncertainty, ranging from 20 − 100%, is the statistical precision on the number of events with at least four baseline leptons and at least one failing to satisfy the signal-lepton requirement.

Polarisation measurements
The  boson can be either transversely polarised or longitudinally polarised, and the polarisation fractions depend on the transverse momentum of the  boson [85].These effects lead to different kinematic properties of the production and the final decay states of the -boson pair.To extract the fraction of  L  L events from the reconstructed   candidate events, a multivariate technique based on a boosted decision tree (BDT) [86] is used to enhance the separation between  L  L and  T  X ( T  T or  T  L ) events.After a dedicated optimisation study to maximise the  L  L signal sensitivity, the input variables used in the BDT are the following: cos  1 (cos  3 ), where  1 ( 3 ) is the angle between the negatively charged final-state lepton in the  1 ( 2 ) rest frame and the direction of flight of the  1 ( 2 ) boson in the four-lepton rest frame; cos  *  1 , where  *  1 is the production angle of the  1 defined in the four-lepton rest frame; and Δ ℓ 1 ℓ 2 (Δ ℓ 3 ℓ 4 ), the azimuthal separation of the two leptons from  1 ( 2 ) defined in the four-lepton rest frame.The angles are illustrated in Figure 2. Other kinematic variables, such as the  T and rapidity of  1 and  2 also have substantial separation power, but they are not included in the BDT training to reduce theoretical modelling uncertainties.The MC templates for different   polarisation states were generated at LO in QCD.Higher-order corrections, in both the QCD and EW, on MC templates of the different   polarisation states have to be taken into account when extracting the  L  L fraction from data.Recently the combined NLO QCD and EW corrections and the loop-induced  →   corrections to the polarisation structure of   production were calculated at fixed-order with the MoCaNLO program [66].Differential cross-sections for different polarisation states of the  q →   and  →   processes were provided for several kinematic observables.As in the simulated MC samples, the polarisation definition in the MoCaNLO program is also based on the CM frame of the two  bosons.In order to incorporate the fixed-order higher-order corrections into the polarisation measurement, a three-step reweighting method is established, using either a one-dimensional (1D) observable or two-dimensional (2D) observables.
• Step 1: 1D reweighting for each individual polarisation state.In this step, the reweighting is done separately for the  q →   and  →   processes.For  q →  , the combined NLO QCD and EW corrections, in a multiplicative approach, as a function of the cos  1 variable are applied by taking the ratio of the differential cross-sections calculated from MoCaNLO at NLO and the ones from the MadGraph5_aMC@NLO MC samples at the particle level in the fiducial phase space of the measurement, for  T  T ,  T  L and  L  L events, respectively.An additional reweighting as a function of cos  1 is applied to account for the missing higher-order QCD and parton shower effects by taking the ratio of inclusive Sherpa  q →   predictions at the particle level to the inclusive MoCaNLO calculations.The impact of these two 1D reweighting corrections on the BDT discriminant is mostly of about 20% on the normalisation, and the shape variation is only of 2 − 4%.For the  →   process, which contributes to about 15% of the total signal yield, only the inclusive Sherpa MC sample is available, and the MoCaNLO program provides polarised differential cross-sections at LO. Thus the inclusive Sherpa MC sample is reweighted to obtain polarised templates of  T  T ,  T  L and  L  L , by taking the fraction of polarised and inclusive cross-sections calculated by MoCaNLO as a function of cos  1 .No reweighting is applied to the EW  →     process and the original MC simulation of polarised samples is used.
• Step 2: 1D reweighting for the interference effect.The simulated polarised samples do not consider the interference effects among different polarisation states, while such interference effects are found to be non-negligible in some kinematic regions where the contribution could reach up to 5% [66].A dedicated template for the interference term is therefore constructed by reweighting the inclusive Sherpa  q →   events with MoCaNLO calculations that include interference contributions, by taking the difference between the inclusive cross-sections and the sum of the three polarised cross-sections as a function of cos  1 .For  →   events, the interference effect is found to be negligible and thus ignored.For the subleading EW  →     process, the interference effect is not included either.
• Step 3: 2D reweighting for the residual higher-order corrections.Four templates, including three polarisation states and the interference term, are obtained after the two reweighting steps described above.A closure test is performed by comparing the sum of the four templates and the prediction given by the inclusive Sherpa MC events.Residual discrepancies are observed, which could be due to the non-closure of the 1D reweighting method, resulting from missing higher-order QCD and parton shower effects.An additional 2D reweighting is applied to each of the three polarisation templates to correct the mismodelling by taking the ratio of the non-closure effect and the sum of the three polarisation templates as a function of cos  *  1 and Δ ℓ 1 ℓ 2 .The impact of this 2D reweighting on the BDT discriminant is mostly on the shape with a maximum variation of about 10%.
Figure 3 shows the BDT distribution of the three polarisation templates for the  q →   process before and after the reweighting procedure.To extract the fraction of  L  L events, a profile binned maximumlikelihood fit [87][88][89] to the BDT distribution is performed using the final templates for the three polarisation states and the interference term, and other non-  background contributions.The normalisation factors of the  L  L template,   , referred as the signal strength, and the combined  T  T +  T  L templates,    , are allowed to float in the fit.The signal strength is the ratio of the measured signal contribution relative to the SM expectation.A validation study is performed to show the robustness of the templates, by applying a similar fit to the inclusive Powheg Box MC events.The extracted  L  L fraction is compared with the predicted  L  L fraction, and they are found to be consistent within uncertainties.  Figure 3: BDT distributions of the three polarisation templates for the  q →   process, before (dashed lines) and after (solid lines) the reweighting procedure to account for higher-order corrections.All distributions are normalised to the same area.The lower panel shows the ratio of the templates after the corrections to those before the corrections.

CP-odd Optimal Observable
The OO defined for the CP study combines the CP-sensitive polar and azimuthal angles of both -boson systems, providing additional CP sensitivity from shape differences between the SM and aNTGC predictions.The CP-sensitive polar angles  1 ( 3 ) for the  1 ( 2 ) boson are already defined in Section 6.1 and illustrated in Figure 2. The CP-sensitive azimuthal angles  1 and  3 are reconstructed in a reference frame illustrated in Figure 2 that allows a direct measure of the -boson spin as discussed in Ref. [24,90].The CP-sensitive azimuthal angle  1 ( 3 ) is the azimuthal angle of the negative lepton in the  1 ( 2 ) rest frame in this new axis system.The differential cross-sections for  1 ( 3 ) and  1 ( 3 ) are symmetric in the SM but asymmetric in the presence of a CP-odd aNTGC.
To improve the sensitivity, the two CP-sensitive angles  1 ( 3 ) and  1 ( 3 ) are combined to form an angular observable  ,1(3) = sin  1(3) × cos  1(3) that enhances the asymmetry for each -boson system.sensitive and non-sensitive bins to enhance the sensitivity for the four-lepton system.Each bin of the O  ,1  ,3 observable represents approximately an L-shaped grouping of the bins around the  ,3 =  ,1 line as shown in Figure 5(a).The small fraction of events with mis-paired leptons in the   → 4 (4) final states was studied and found to have negligible impact on the CP-sensitivity of the OO.The data agree closely with the prediction within the measurement's statistical precision and systematic uncertainties.Figure 5(b) also shows an asymmetric prediction in the presence of a CP odd BSM coupling when  4  = 1.

Detector corrections
Particle-level differential cross-sections for the on-shell   production are obtained by correcting the detector effects such as inefficiency and resolution.The background-subtracted event yields are corrected using an iterative Bayesian unfolding method [91].
The first step of the correction multiplies each bin yield by a fiducial correction factor obtained from the Sherpa SM prediction, which accounts for the events that satisfy the detector level but fail to satisfy the fiducial-level event selections.This correction accounts for the 5 − 20% of the fake fiducial events in various bins caused by the resolution effects.Then, the detector resolution-induced bin migrations are corrected iteratively using the SM particle-level distribution as the initial prior.With an increasing number of iterations, the statistical uncertainty increases, and the residual bias relative to the prior decreases due to the improvement of its knowledge.Two iterations were deemed optimal as a compromise between the increasing statistical uncertainty and decreasing bias.The final step in the unfolding procedure is to correct for the detector inefficiency by dividing the per-bin yield by the ratio of the number of events satisfying both the particle-and reconstruction-level selections to the number of events passing the particle-level selections.
A data-driven closure test is performed to evaluate the model dependence of the unfolding method.This test first simulates a pseudo-data sample by reweighting the SM prediction to the shape observed in the data.The pseudo-data sample is then unfolded using the nominal SM prediction.The comparison of the unfolded pseudo-data with the reweighted particle-level prediction gives the intrinsic bias of the unfolding method, which was found to be less than 1% in each bin of the unfolded O  ,1  ,3 observable in the case of two iterations of unfolding.This resulting bias is taken as a systematic uncertainty of the final result.Moreover, the uncertainty related to the choice of the generator in the unfolding is studied using the alternative Powheg prediction of the  q →   process, which is reweighted to match the nominal Sherpa lineshape to avoid double counting of the data-driven bias.The generator bias estimated by comparing the difference between the unfolded results is negligible.
Additionally, an injection test is performed to evaluate the robustness of the unfolding algorithm in the presence of BSM physics in the data.A detector-level distribution for the BSM aNTGC parameter   4 = 1 is injected into the SM detector-level prediction.The BSM-injected detector level distribution is then unfolded using the inputs from the nominal SM prediction.When compared, the unfolded distribution agrees closely with the corresponding particle-level distribution within uncertainties.

Systematic uncertainties
Both the polarisation and CP property studies presented here are affected by some common sources of theoretical, experimental, and background-related uncertainties.
Theoretical systematic uncertainties affecting the results of this analysis arise from two sources: uncertainties related to the QCD scale dependence and uncertainties related to the choice of the PDF set.The impact of each is propagated to the final results.The QCD scale dependency is evaluated individually for three physics processes  q →  ,  →  , and the EW  →     by varying the default choice of renormalisation and factorisation scales independently by factors of two and one-half, and then removing combinations where the variations differ by a factor of four.For each process, the envelope of the effects from these variations is taken as the scale uncertainty.
The PDF-related uncertainty for each of the signal samples is estimated using the PDF4LHC recommendation [92].The PDF variations include a set of 100 replica variations on the nominal NNPDF set, two additional variations from the alternative PDF sets of MMHT2014nnlo [93] and CT14nnlo [94], and variations of the strong coupling constant by ±0.001 around the nominal value of   = 0.118.The total PDF uncertainty is the absolute envelope of standard deviations of 100 internal variations and of the two alternate PDF variations, added in quadrature with the envelope of the   variations.
For the polarisation measurements, dedicated systematic uncertainties are considered in the modelling of the polarisation templates, including theoretical uncertainties from higher-order corrections, PDFs and   described above, and the uncertainties associated with the reweighting methods as described in Section 6.1.For both the  q →   and  →   polarisation templates, the theoretical uncertainties from QCD scales are taken from the MoCaNLO calculations by varying the QCD scales as described above, while for the template of the interference term, which is reweighted from the inclusive Sherpa  q →   samples, these theoretical uncertainties are estimated from the Sherpa sample.The parton showering and hadronisation uncertainty is estimated for the signal by comparing the nominal Pythia 8 parton showering with the alternative Herwig 7 [95,96] algorithm.Uncertainties from the NLO EW corrections are estimated by taking the difference between the additive and the multiplicative prescription in MoCaNLO [66].For the 1D reweighting method in Section 6.1, the uncertainty is estimated by comparing the reweighted templates using the nominal observable and an alternative observable: the rapidity difference of the two  bosons.For the additional 2D reweighting, the residual difference between the sum of the four templates and the inclusive prediction from the Sherpa sample is taken as an uncertainty.
For the CP study, additional uncertainties associated with the higher-order corrections applied to the inclusive  q →   and  →   samples are considered.For the virtual NLO electroweak effects on the  q →   process, a 100% uncertainty is assigned to the reweighting function to account for non-factorisable effects in events with high QCD activity [97].The uncertainty in the NLO QCD -factors for the  →   process is evaluated differentially as a function of the  4ℓ as discussed in Ref. [58].
The uncertainties associated with reconstruction, identification, isolation and track-to-vertex matching efficiencies, and momentum resolution and scale of the leptons are the dominant experimental uncertainties and originate from imperfect modelling in the simulation and uncertainties in the determination of the correction factors.The uncertainties associated with each scale factor that is applied in the simulation are estimated by modifying the nominal values by their associated uncertainties [77,98].
An uncertainty of ±0.83% from the measurement of Run 2 data sample luminosity [99] is propagated to the final results.A systematic uncertainty related to the pile-up reweighting is included to cover the uncertainty in the ratio of the predicted and measured   inelastic cross-sections [100].
The background-related uncertainties are from two distinct sources, uncertainties related to the nonprompt background estimate, discussed in Section 5, and the uncertainty related to the background containing four prompt leptons simulated from MC.The simulation of  t and triboson background processes are normalised to ATLAS measurements, as outlined in Section 3, and the uncertainty related is estimated by varying the normalisation of the simulated samples by the experimental precision of the ATLAS measurements.
For the CP study, the systematic uncertainties are propagated to the particle level through the unfolding method.The effect of each systematic source is evaluated using simulation.Since the theory variation applied affects both the detector and the particle-level yields, the resulting theory uncertainties on the unfolded cross-sections are minor.Another source of systematic uncertainty related to the intrinsic unfolding bias discussed in Section 6.2.2 is also included in the differential cross-sections and the BSM interpretation.

Polarisation measurements
The profile likelihood fit procedure described in Section 6.1 is performed on the BDT distribution of the data.Systematic uncertainties are modelled as Gaussian-constrained nuisance parameters in the likelihood.Figure 6 shows the post-fit BDT distribution in the signal region.The corresponding pre-fit and post-fit yields in the signal region are detailed in Table 1.The decrease of the yield uncertainty for  T  T is mainly due to the constraint on the normalisation from the fit and the increase of the yield uncertainty for  L  L is mainly due to the statistical uncertainty of the signal strength after the fit.The  L  L signal strength is measured to be   = 1.15 ± 0.27(stat.)± 0.11(syst.)= 1.15 ± 0.29, where the uncertainties are either statistical (stat.)or of systematic (syst.)nature.This corresponds to a significance of 4.3.Figure 7 shows the profile likelihood ratio [101] as a function of   .The measured results are consistent with the SM expectation of   = 1.00 ± 0.27 with an expected significance of 3.8.The normalisation factor for the  T  X template is measured to be    = 1.00 ± 0.05, also consistent with the SM prediction.
Table 1: Expected and observed numbers of events in the signal region.Numbers are presented before and after the fit to the BDT distribution.The 'Others' category represents the contribution from  t and  .For the  L  L signal, the pre-fit yield values correspond to the theoretical prediction and corresponding uncertainties.The uncertainties include both of the statistical and systematic contributions.The uncertainty in the total yield can be smaller than the quadrature sum of the contributions because of correlations resulting from the fit.

Pre-fit
Post-fit   fiducial cross-section of the  L  L production is measured to be  obs.
L  L = 2.10 ± 0.09 fb.The SM prediction includes NLO QCD and EW corrections for the  q →   process, the LO prediction for the  →   process, both of which are calculated from MoCaNLO, and the LO prediction for the EW  →     process.The uncertainty on the prediction is dominated by QCD scale and PDF uncertainties.The measurement is limited by data statistical uncertainty and the impact of uncertainties in the measured  L  L fiducial cross-section is shown in Table 2, with the leading contribution from the theoretical modelling of the polarisation templates.

Unfolded differential measurement
The differential cross-section for the   → 4ℓ production as a function of the OO O  ,1  ,3 is shown in Figure 8 and compared with two different SM predictions, where the  q →   contribution is predicted either using the Sherpa or Powheg generators.The band represents the total uncertainties, which are small and flat, related to the theoretical, experimental, and unfolding procedure-related uncertainties on the unfolded differential cross-section.Similarly, for the SM predictions, the uncertainty bands represent the total theoretical uncertainties on the total prediction.For each bin, the measured cross-section from the data agrees closely with both sets of predicted cross-sections.  .The grey, red and blue uncertainty bands represent the respective systematic uncertainties in the measured unfolded cross-section and predicted particle-level cross-sections, where the  q →   contribution is predicted using either the Sherpa or Powheg generators.The vertical error bars represent the total uncertainties in the measured differential cross-section.

BSM interpretation
The differential cross-section as a function of O  ,1  ,3 and the corresponding SM + aNTGC prediction are used to define a likelihood function, where ì  meas.and ì  pred.are the measured and predicted differential cross-sections as a function of O  ,1  ,3 , respectively, and  is the total covariance matrix defined by the sum of the statistical and systematic covariances in ì  meas.and ì  pred. .Each systematic uncertainty is treated as fully correlated across bins but uncorrelated with other uncertainty sources.Each source of theoretical uncertainty affecting the SM particle level prediction of the  q →  ,  →   and EW  →     processes is implemented as nuisance parameters, ì , with Gaussian constraints, G (  , 0, 1).
Limits on the two CP-odd aNTGC parameters can be constrained in two scenarios, either including or excluding the quadratic term of Equation 1.The main focus of the measurement is the linear interferenceonly terms, which results in an asymmetric differential cross-section.However, the difference in limits between the linear only and linear + quadratic terms indicates the role played by the neglected operators at higher dimension in the effective Lagrangian [28].Thus, limits are set in these two scenarios for each CP-odd aNTGC parameter under test when the other is set to zero.The test statistics are constructed based on the profile likelihood ratio [101], where  is the aNTGC parameter, ĉ and ì  are the unconditional maximum-likelihood estimators, and ì  () is the conditional maximum-likelihood estimator under the  hypothesis.
Using the test statistic in Equation (3), which is assumed to be distributed according to a  2 distribution with one degree of freedom following Wilks' theorem [102], a 95% confidence level for each aNTGC parameter is calculated, as shown in Table 3. Theoretical uncertainties related to the QCD scale, PDF,   , and parton showering for the  q →   predictions have the largest impact on the estimated confidence interval.
From Table 3, the constraints on the CP-odd  4   and  4  aNTGC parameters using the linear interference term when the quadratic term is set to zero are looser than those obtained using high- T kinematic observables with quadratic terms included [23], as expected.However, these are the first constraints on CP-odd aNTGC parameters using only linear interference terms with a dedicated CP-sensitive OO, based on angular observables.Additional likelihood fits are performed with the quadratic terms included, and the constraints on  4   and  4  are stronger by one order of magnitude, but since O  ,1  ,3 is not sensitive to the high- T regime, such constraints are still not tighter than those obtained with high- T kinematic observables.
Table 3: Expected and observed 95% confidence interval for the two aNTGC operators using the measured particlelevel differential cross-section as a function of O  ,1  ,3 .For the 'Full' case, both the interference and the quadratic terms are included in the aNTGC prediction.

Conclusion
Studies of the polarisation and CP properties in   → 4ℓ production are presented for the first time, using proton-proton collisions at the LHC collected with the ATLAS detector at √  = 13 TeV with an integrated luminosity of 140 fb −1 .The -boson candidates are reconstructed with same-flavour, opposite-charge electron or muon pairs, and they are required to be on-shell with | ℓℓ −   | < 10 GeV.The production of two simultaneously longitudinally polarised  bosons is measured with observed and expected significances of 4.3 and 3.8 standard deviations, respectively.The production cross-section of  L  L events is measured in a fiducial phase space as 2.45 ± 0.60 fb, consistent with the Standard Model prediction of 2.10 ± 0.09 fb with NLO QCD and EW corrections considered.The differential cross-section of the inclusive   boson production as a function of a CP-sensitive angular observable, O  ,1  ,3 , is also measured.Furthermore, the measured differential cross-section is used to constrain the CP-odd neutral triple gauge couplings      [101] G. Cowan, K. Cranmer

Figure 2 :
Figure2: Definition of the angles used for the polarisation measurement and the reference frame used to define the CP-sensitive angles.The -frame (dot-dashed) is the laboratory frame with the -axis along the beam direction.The  ′  ′  ′ -frame (solid) is a new frame used to define the CP-sensitive angles.The  ′ -axis is defined as the direction of motion of the  1 boson in the four-lepton rest frame.The  ′ -axis defines the reaction plane containing the laboratory -axis and the  ′ -axis.The right-hand rule gives the  ′ -axis.

Figures 4 (Figure 4 :
Figure 4: Particle level 2D differential cross-sections of   of the two  bosons for the  q →   → 4ℓ process as predicted by (a) the SM and (b) in the presence of the BSM aNTGC vertex.The BSM prediction shows the contribution of the interference effects only, excluding the quadratic term in Equation (1), when  4  = 1.

Figure 5 (
Figure 5(b) shows the data compared with the total SM signal and background MC prediction at the detector level of the OO O  ,1  ,3 .The bins 1 to 7 and 24 to 30 in Figure 5(b) represent the third quadrant and the first quadrant, respectively, of the 2D distribution of  ,1 vs  ,3 shown in Figure 4.In these two quadrants, the   observables for both  bosons have the same sign in the SM and are the most CP-sensitive region, along with the two central bins representing the bin number 15 and 16 of the OO.The data agree closely with the prediction within the measurement's statistical precision and systematic uncertainties.Figure 5(b) also shows an asymmetric prediction in the presence of a CP odd BSM coupling when  4  = 1.

Figure 5 :
Figure 5: (a) The 2 → 1 mapping and (b) the detector-level measurement of the OO O  ,1  ,3.The measured distribution is compared with the SM signal prediction and the total background.The 'Others' category includes the contribution from  t and   processes.The non-prompt background is estimated by using the fake-factor method.The grey band represents the effect of the total theoretical and experimental uncertainties for the detector-level predictions, and the vertical error bars on data represent the statistical uncertainties.The effect of CP-odd BSM coupling is also represented by the dashed histogram.

Figure 6 :
Figure6: The BDT distribution of the data and the post-fit SM predictions.The 'Others' category represents the contribution from  t and  .The lower panel shows the ratio of the data points to the post-fit total prediction.The arrow indicates that the ratio lies outside the range covered by the vertical axis.The uncertainty bands include both of the statistical and systematic uncertainties as obtained by the fit.

Figure 7 :
Figure 7: Observed and expected profile likelihood ratio, −2lnΛ, as a function of the signal strength   , with (solid curves) and without (dashed curves) systematic uncertainties included.
An additional profile likelihood fit is performed to convert the measured   to the measured fiducial crosssection, with the normalisation effects from the theoretical uncertainties of the  L  L template removed.The

Table 2 :
Impact of uncertainties in the measured fiducial cross-section   L  L from the fit.The impact from a group of nuisance parameters is defined via quadrature subtraction of the   L  L uncertainties: √︁ (Δ) 2 − (Δ ′ ) 2 , where Δ is the uncertainty of the   L  L from the nominal fit and the Δ ′ is the uncertainty of the   L  L when this group of nuisance parameters are fixed to their best-fit values from the nominal fit.Unfolded differential cross-section as a function of the Optimal Observable O  ,1  ,3