Measurement of the angular coefficients in $Z$-boson events using electron and muon pairs from data taken at $\sqrt{s}=8$ TeV with the ATLAS detector

The angular distributions of Drell-Yan charged lepton pairs around the $Z$-boson mass peak probe the underlying QCD dynamics of $Z$-boson production. This paper presents a measurement of the complete set of angular coefficients $A_{0-7}$ describing these distributions in the $Z$-boson Collins-Soper frame. The data analysed correspond to 20.3 fb$^{-1}$ of $pp$ collisions at $\sqrt{s} = 8$ TeV, collected by the ATLAS detector at the CERN LHC. The measurements are compared to the most precise fixed-order calculations currently available ($\mathcal{O}(\alpha^{2}_{s})$) and with theoretical predictions embedded in Monte Carlo generators. The measurements are precise enough to probe QCD corrections beyond the formal accuracy of these calculations and to provide discrimination between different parton-shower models. A significant deviation from the $\mathcal{O}(\alpha^{2}_{s})$ predictions is observed for $A_{0}-A_{2}$. Evidence is found for non-zero $A_{5,6,7}$, consistent with expectations.


Introduction
The angular distributions of charged lepton pairs produced in hadron-hadron collisions via the Drell-Yan neutral current process provide a portal to precise measurements of the production dynamics through spin correlation effects between the initial-state partons and the final-state leptons mediated by a spin-1 intermediate state, predominantly the Z boson. In the Z-boson rest frame, a plane spanned by the directions of the incoming protons can be defined, e.g. using the Collins-Soper (CS) reference frame [1]. The lepton polar and azimuthal angular variables, denoted by cos θ and φ in the following formalism, are defined in this reference frame. The spin correlations are described by a set of nine helicity density matrix elements, which can be calculated within the context of the parton model using perturbative quantum chromodynamics (QCD). The theoretical formalism is elaborated in Refs. [2][3][4][5].
The full five-dimensional differential cross-section describing the kinematics of the two Born-level leptons from the Z-boson decay can be decomposed as a sum of nine harmonic polynomials, which depend on cos θ and φ, multiplied by corresponding helicity cross-sections that depend on the Z-boson transverse momentum (p Z T ), rapidity (y Z ), and invariant mass (m Z ). It is a standard convention to factorise out the unpolarised cross-section, denoted in the literature by σ U+L , and to present the five-dimensional differential cross-section as an expansion into nine harmonic polynomials P i (cos θ, φ) and dimensionless angular coefficients A 0−7 (p Z T , y Z , m Z ), which represent ratios of helicity cross-sections with respect to the unpolarised one, σ U+L , as explained in detail in Appendix A: dσ dp Z T dy Z dm Z d cos θ dφ = 3 16π dσ U+L dp Z T dy Z dm Z (1 + cos 2 θ) + 1 2 A 0 (1 − 3 cos 2 θ) + A 1 sin 2θ cos φ (1) + 1 2 A 2 sin 2 θ cos 2φ + A 3 sin θ cos φ + A 4 cos θ +A 5 sin 2 θ sin 2φ + A 6 sin 2θ sin φ + A 7 sin θ sin φ .
The dependence of the differential cross-section on cos θ and φ is thus completely manifest analytically.
In contrast, the dependence on p Z T , y Z , and m Z is entirely contained in the A i coefficients and σ U+L . Therefore, all hadronic dynamics from the production mechanism are described implicitly within the structure of the A i coefficients, and are factorised from the decay kinematics in the Z-boson rest frame. This allows the measurement precision to be essentially insensitive to all uncertainties in QCD, quantum electrodynamics (QED), and electroweak (EW) effects related to Z-boson production and decay. In particular, EW corrections that couple the initial-state quarks to the final-state leptons have a negligible impact (below 0.05%) at the Z-boson pole. This has been shown for the LEP precision measurements [6,7], when calculating the interference between initial-state and final-state QED radiation.
At leading order (LO) in QCD, only the annihilation diagram qq → Z is present and only A 4 is non-zero. At next-to-leading order (NLO) in QCD (O(α s )), A 0−3 also become non-zero. The Lam-Tung relation [8][9][10], which predicts that A 0 − A 2 = 0 due to the spin-1 of the gluon in the qg → Zq and qq → Zg diagrams, is expected to hold up to O(α s ), but can be violated at higher orders. The coefficients A 5,6,7 are expected to become non-zero, while remaining small, only at next-to-next-to-leading order (NNLO) in QCD (O(α 2 s )), because they arise from gluon loops that are included in the calculations [11,12]. The coefficients A 3 and A 4 depend on the product of vector and axial couplings to quarks and leptons, and are sensitive to the Weinberg angle sin 2 θ W . The explicit formulae for these dependences can be found in Appendix A.
The full set of coefficients has been calculated for the first time at O(α 2 s ) in Refs. [2][3][4][5]. More recent discussions of these angular coefficients may be found in Ref. [13], where the predictions in the NNLOPS scheme of the Powheg [14][15][16][17] event generator are shown for Z-boson production, and in Ref. [18], where the coefficients are explored in the context of W-boson production, for which the same formalism holds.
The CDF Collaboration at the Tevatron published [19] a measurement of some of the angular coefficients of lepton pairs produced near the Z-boson mass pole, using 2.1 fb −1 of proton-anti-proton collision data at a centre-of-mass energy √ s = 1.96 TeV. Since the measurement was performed only in projections of cos θ and φ, the coefficients A 1 and A 6 were inaccessible. They further assumed A 5 and A 7 to be zero since the sensitivity to these coefficients was beyond the precision of the measurements; the coefficients A 0,2,3,4 were measured as a function of p Z T . These measurements were later used by CDF [20] to infer an indirect measurement of sin 2 θ W , or equivalently, the W-boson mass in the on-shell scheme, from the average A 4 coefficient. These first measurements of the angular coefficients demonstrated the potential of this not-yet-fully explored experimental avenue for investigating hard QCD and EW physics.
Measurements of the W-boson angular coefficients at the LHC were published by both ATLAS [21] and CMS [22]. More recently, a measurement of the Z-boson angular coefficients with Z → µµ decays was published by CMS [23], where the first five coefficients were measured with 19.7 fb −1 of protonproton (pp) collision data at √ s = 8 TeV. The measurement was performed in two y Z bins, 0 < |y Z | < 1 and 1 < |y Z | < 2.1, each with eight bins in p Z T up to 300 GeV. The violation of the Lam-Tung relation was observed, as predicted by QCD calculations beyond NLO. This paper presents an inclusive measurement of the full set of eight A i coefficients using charged lepton pairs (electrons or muons), denoted hereafter by . The measurement is performed in the Z-boson invariant mass window of 80-100 GeV, as a function of p Z T , and also in three bins of y Z . These results are based on 20.3 fb −1 of pp collision data collected at √ s = 8 TeV by the ATLAS experiment [24] at the LHC. With the measurement techniques developed for this analysis, the complete set of coefficients is extracted with fine granularity over 23 bins of p Z T up to 600 GeV. The measurements, performed in the CS reference frame [1], are first presented as a function of p Z T , integrating over y Z . Further measurements divided into three bins of y Z are also presented: 0 < |y Z | < 1, 1 < |y Z | < 2, and 2 < |y Z | < 3.5. The Z/γ * → e + e − and Z/γ * → µ + µ − channels where both leptons fall within the pseudorapidity range |η| < 2.4 (hereafter referred to as the central-central or ee CC and µµ CC channels) are used for the y Z -integrated measurement and the first two y Z bins. The Z/γ * → e + e − channel where one of the electrons instead falls in the region |η| > 2.5 (referred to hereafter as the central-forward or ee CF channel) is used to extend the measurement to the high-y Z region encompassed by the third y Z bin. In this case, however, because of the fewer events available for the measurement itself and to evaluate the backgrounds (see Section 4), the measurement is only performed for p Z T up to 100 GeV using projections of cos θ and φ, making A 1 and A 6 inaccessible in the 2 < |y Z | < 3.5 bin.
The high granularity and precision of the specific measurements presented in this paper provide a stringent test of the most precise perturbative QCD predictions for Z-boson production in pp collisions and of Monte Carlo (MC) event generators used to simulate Z-boson production. This paper is organised as follows. Section 2 summarises the theoretical formalism used to extract the angular coefficients and presents the fixed-order QCD predictions for their variations as a function of p Z T . Section 3 describes briefly the ATLAS detector and the data and MC samples used in the analysis, while Section 4 presents the data analysis and background estimates for each of the three channels considered. Section 5 describes the fit methodology used to extract the angular coefficients in the full phase space as a function of p Z T and Section 6 gives an overview of the statistical and systematic uncertainties of the measurements. Sections 7 and 8 present the results and compare them to various predictions from theoretical calculations and MC event generators, and Section 9 summarises and concludes the paper.

Theoretical predictions
The differential cross-section in Eq. (1) is written for pure Z bosons, although it also holds for the contribution from γ * and its interference with the Z boson. The tight invariant mass window of 80-100 GeV is chosen to minimise the γ * contribution, although the predicted A i coefficients presented in this paper are effective coefficients, containing this small contribution from γ * . This contribution is not accounted for explicitly in the detailed formalism described in Appendix A, which is presented for simplicity for pure Z-boson production. Throughout this paper, the leptons from Z-boson decays are defined at the Born level, i.e. before final-state QED radiation, when discussing theoretical calculations or predictions at the event-generator level.
The p Z T and y Z dependence of the coefficients varies strongly with the choice of spin quantisation axis in the Z-boson rest frame (z-axis). In the CS reference frame adopted for this paper, the z-axis is defined in the Z-boson rest frame as the external bisector of the angle between the momenta of the two protons, as depicted in Fig. 1. The positive direction of the z-axis is defined by the direction of positive longitudinal Z-boson momentum in the laboratory frame. To complete the coordinate system, the y-axis is defined as the normal vector to the plane spanned by the two incoming proton momenta and the x-axis is chosen to define a right-handed Cartesian coordinate system with the other two axes. Polar and azimuthal angles are calculated with respect to the negatively charged lepton and are labelled θ CS and φ CS , respectively. In the case where p Z T = 0, the direction of the y-axis and the definition of φ CS are arbitrary. Historically, there has been an ambiguity in the definition of the sign of the φ CS angle in the CS frame: this paper adopts the recent convention followed by Refs. [13,23], whereby the coefficients A 1 and A 3 are positive.
The coefficients are not explicitly used as input to the theoretical calculations nor in the MC event generators. They can, however, be extracted from the shapes of the angular distributions with the method proposed in Ref. [3], owing to the orthogonality of the P i polynomials. The weighted average of the angular distributions with respect to any specific polynomial isolates an average reference value or moment p l p x y ẑ^θ CS φ CS Figure 1: Sketch of the Collins-Soper reference frame, in which the angles θ CS and φ CS are defined with respect to the negatively charged lepton (see text). The notationsx,ŷ andẑ denote the unit vectors along the corresponding axes in this reference frame.
The moment of each harmonic polynomial can thus be expressed as (see Eq. (1)): One thus obtains a representation of the effective angular coefficients for Z/γ * production. These effective angular coefficients display in certain cases a dependence on y Z , which arises mostly from the fact that the interacting quark direction is unknown on an event-by-event basis. As the method of Ref. [3] relies on integration over the full phase space of the angular distributions, it cannot be applied directly to data, but is used to compute all the theoretical predictions shown in this paper.
The inclusive fixed-order perturbative QCD predictions for Z-boson production at NLO and NNLO were obtained with DYNNLO v1.3 [25]. These inclusive calculations are formally accurate to O(α 2 s ). The Z-boson is produced, however, at non-zero transverse momentum only at O(α s ), and therefore the calculation of the coefficients as a function of p Z T is only NLO. Even though the fixed-order calculations do not provide reliable absolute predictions for the p Z T spectrum at low values, they can be used for p Z T > 2.5 GeV for the angular coefficients. The results were cross-checked with NNLO predictions from FEWZ v3.1.b2 [26][27][28] and agreement between the two programs was found within uncertainties. The renormalisation and factorisation scales in the calculations were set to E Z T = (m Z ) 2 + (p Z T ) 2 [29] on an event-by-event basis. The calculations were done using the CT10 NLO or NNLO parton distribution functions (PDFs) [30], depending on the order of the prediction.
The NLO EW corrections affect mostly the leading-order QCD cross-section normalisation in the Z-pole region and have some impact on the p Z T distribution, but they do not affect the angular correlations at the Z-boson vertex. The DYNNLO calculation was done at leading order in EW, using the G µ scheme [31]. This choice determines the value of A 4 at low p Z T , and for the purpose of the comparisons presented in this paper, both A 3 and A 4 obtained from DYNNLO are rescaled to the values predicted when using the measured value of sin 2 θ eff W = 0.23113 [32]. The theoretical predictions are shown in Fig. 2 and tabulated in Table 1 for three illustrative p Z T bins. The binning in p Z T is chosen based on the experimental resolution at low p Z T and on the number of events at high p Z T and has the following boundaries (in GeV) used consistently throughout the measurement: The predictions show the following general features. The A 0 and A 2 coefficients increase as a function of p Z T and the deviations from lowest-order expectations are quite large, even at modest values of p Z T = 20-50 GeV. The A 1 and A 3 coefficients are relatively small even at large p Z T , with a maximum value of 0.08. In the limit where p Z T = 0, all coefficients except A 4 are expected to vanish at NLO. The NNLO corrections are typically small for all coefficients except A 2 , for which the largest correction has a value of −0.08, in agreement with the original theoretical studies [2]. The theoretical predictions for A 5,6,7 are not shown because these coefficients are expected to be very small at all values of p Z T : they are zero at NLO and the NNLO contribution is large enough to be observable, namely of the order of 0.005 for values of p Z T in the range 20-200 GeV.
The statistical uncertainties of the calculations, as well as the factorisation and renormalisation scale and PDF uncertainties, were all considered as sources of theoretical uncertainties. The statistical uncertainties of the NLO and NNLO predictions in absolute units are typically 0.0003 and 0.003, respectively. The larger statistical uncertainties of the NNLO predictions are due to the longer computational time required than for the NLO predictions. The scale uncertainties were estimated by varying the renormalisation and factorisation scales simultaneously up and down by a factor of two. As stated in Ref. [2], the theoretical uncertainties due to the choice of these scales are very small for the angular coefficients because they are ratios of cross-sections. The resulting variations of the coefficients at NNLO were found in most cases to be comparable to the statistical uncertainty. The PDF uncertainties were estimated using the CT10 NNLO eigenvector variations, as obtained from FEWZ and normalised to 68% confidence level. They were found to be small compared to the NNLO statistical uncertainty, namely of the order of 0.001 for A 0−3 and 0.002 for A 4 .  Figure 2: The angular coefficients A 0−4 and the difference A 0 − A 2 , shown as a function of p Z T , as predicted from DYNNLO calculations at NLO and NNLO in QCD. The NLO predictions for A 0 − A 2 are compatible with zero, as expected from the Lam-Tung relation [8][9][10]. The error bars show the total uncertainty of the predictions, including contributions from statistical uncertainties, QCD scale variations and PDFs. The statistical uncertainties of the NNLO predictions are dominant and an order of magnitude larger than those of the NLO predictions.  (5)(6)(7)(8), mid , and high (132-173 GeV) p Z T for the y Z -integrated configuration. The uncertainty represents the sum of statistical and systematic uncertainties.
−0.0003 detector, electromagnetic (EM) and hadronic calorimeters, and a muon spectrometer. The inner tracker provides precision tracking of charged particles in the pseudorapidity range |η| < 2.5. This region is matched to a high-granularity EM sampling calorimeter covering the pseudorapidity range |η| < 3.2 and a coarser granularity calorimeter up to |η| = 4.9. A hadronic calorimeter system covers the entire pseudorapidity range up to |η| = 4.9. The muon spectrometer provides triggering and tracking capabilities in the range |η| < 2.4 and |η| < 2.7, respectively. A first-level trigger is implemented in hardware, followed by two software-based trigger levels that together reduce the accepted event rate to 400 Hz on average. For this paper, a central lepton is one found in the region |η| < 2.4 (excluding, for electrons, the electromagnetic calorimeter barrel/end-cap transition region 1.37 < |η| < 1.52), while a forward electron is one found in the region 2.5 < |η| < 4.9 (excluding the transition region 3.16 < |η| < 3.35 between the electromagnetic end-cap and forward calorimeters). 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 z-axis along the beam pipe. The x-axis points from the IP to the centre of the LHC ring, and the y-axis points upwards. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). Angular distance is measured in units of ∆R ≡ (∆η) 2 + (∆φ) 2 .

Data and Monte Carlo samples
The data were collected by the ATLAS detector in 2012 at a centre-of-mass energy of √ s = 8 TeV, and correspond to an integrated luminosity of 20.3 fb −1 . The mean number of additional pp interactions per bunch crossing (pile-up events) in the data set is approximately 20.
The simulation samples used in the analysis are shown in Table 2. The four event generators used to produce the Z/γ * → signal events are listed in Table 2. The baseline PowhegBox (v1/r2129) sample [14][15][16][17], which uses the CT10 NLO set of PDFs [33], is interfaced to Pythia 8 (v.8.170) [34] with the AU2 set of tuned parameters [35] to simulate the parton shower, hadronisation and underlying event, and to Photos (v2.154) [36] to simulate QED final-state radiation (FSR) in the Z-boson decay. The alternative signal samples are from PowhegBox interfaced to Herwig (v.6.520.2) [37] for the parton shower and hadronisation, Jimmy (v4.31) [38] for the underlying event, and Photos for FSR. The Sherpa (v.1.4.1) [39][40][41][42] generator is also used, and has its own implementation of the parton shower, hadronisation, underlying event and FSR, and uses the CT10 NLO PDF set. These alternative samples are used to test the dependence of the analysis on different matrix-element calculations and parton-shower models, as discussed in Section 6. The Powheg (v2.1) + MiNLO event generator [43] was used for the Z+jet process at NLO to normalise certain reference coefficients for the ee CF analysis, as described in Section 5. The number of events available in the baseline PowhegBox + Pythia 8 signal sample corresponds to approximately 4 (25) times that in the data below (above) p Z T = 105 GeV. Backgrounds from EW (diboson and γγ → production) and top-quark (production of top-quark pairs and of single top quarks) processes are evaluated from the MC samples listed in Table 2. The W + jets contribution to the background is instead included in the data-driven multijet background estimate, as described in Section 4; W-boson samples listed in Table 2 are thus only used for studies of the background composition.
All of the samples are processed with the Geant4-based simulation [44] of the ATLAS detector [45]. The effects of additional pp collisions in the same or nearby bunch crossings are simulated by the addition of so-called minimum-bias events generated with Pythia 8.

Data analysis 4.1. Event selection
As mentioned in Sections 1 and 3, the data are split into three orthogonal channels, namely the ee CC channel with two central electrons, the µµ CC channel with two central muons, and the ee CF channel with one central electron and one forward electron. Selected events are required to be in a data-taking period in which the beams were stable and the detector was functioning well, and to contain a reconstructed primary vertex with at least three tracks with p T > 0.4 GeV.
Candidate ee CC events are obtained using a dielectron trigger requiring two electron candidates with p T > 12 GeV, combined with high-p T single-electron triggers. Electron candidates are required to have p T > 25 GeV and are reconstructed from clusters of energy in the electromagnetic calorimeter matched to inner detector tracks. The electron candidates must satisfy a set of "medium" selection criteria [50,51], which have been optimised for the level of pile-up present in the 2012 data. Events are required to contain exactly two electron candidates of opposite charge satisfying the above criteria. Candidate µµ CC events are retained for analysis using a dimuon trigger requiring two muon candidates with p T > 18 GeV and 8 GeV, respectively, combined with single high-p T muon triggers. Muon candidates are required to have p T > 25 GeV and are identified as tracks in the inner detector which are matched and combined with track segments in the muon spectrometer [52]. Track-quality and longitudinal and transverse impact-parameter requirements are imposed for muon identification to suppress backgrounds, and to ensure that the muon candidates originate from a common primary pp interaction vertex. Events are required to contain exactly two muon candidates of opposite charge satisfying the above criteria.
Candidate ee CF events are obtained using a single-electron trigger, requiring an isolated central electron candidate with p T > 24 GeV, combined with a looser high-p T single-electron trigger. The central electron candidate is required to have p T > 25 GeV. Because the expected background from multijet events is larger in this channel than in the ee CC channel, the central electron candidate is required to satisfy a set of "tight" selection criteria [50], which are optimised for the level of pile-up observed in the 2012 data. The forward electron candidate is required to have p T > 20 GeV and to satisfy a set of "medium" selection criteria, based only on the shower shapes in the electromagnetic calorimeter [50] since this region is outside the acceptance of the inner tracker. Events are required to contain exactly two electron candidates satisfying the above criteria.
Since this analysis is focused on the Z-boson pole region, the lepton pair is required to have an invariant mass (m ) within a narrow window around the Z-boson mass, 80 < m < 100 GeV. Events are selected for y Z -integrated measurements without any requirements on the rapidity of the lepton pair (y ). For the y Z -binned measurements, events are selected in three bins of rapidity: |y | < 1.0, 1.0 < |y | < 2.0, and 2.0 < |y | < 3.5. Events are also required to have a dilepton transverse momentum (p T ) less than the value of 600 (100) GeV used for the highest bin in the ee CC and µµ CC (ee CF ) channels. The variables m , y , and p T , which are defined using reconstructed lepton pairs, are to be distinguished from the variables m Z , y Z , and p Z T , which are defined using lepton pairs at the Born level, as described in Section 2. The simulated events are required to satisfy the same selection criteria, after applying small corrections to account for the differences between data and simulation in terms of reconstruction, identification and trigger efficiencies and of energy scale and resolution for electrons and muons [50][51][52][53]. All simulated  Figure 3: Comparison of the expected yields (left) and acceptance times efficiency of selected events (right) as a function of y Z (top) and p Z T (bottom), for the ee CC , µµ CC , and ee CF events. Also shown are the expected yields at the event generator level over the full phase space considered for the measurement, which corresponds to all events with a dilepton mass in the chosen window, 80 < m Z < 100 GeV. events are reweighted to match the distributions observed in data for the level of pile-up and for the primary vertex longitudinal position. Figure 3 illustrates the different ranges in p Z T and y Z expected to be covered by the three channels along with their acceptance times selection efficiencies, which is defined as the ratio of the number of selected events to the number in the full phase space. The difference in shape between the ee CC and µµ CC channels arises from the lower reconstruction and identification efficiency for central electrons at high values of |η| and from the lower trigger and reconstruction efficiency for muons at low values of |η|. The central-central and central-forward channels overlap in the region 1.5 < |y Z | < 2.5.

Backgrounds
In the Z-boson pole region, the backgrounds from other processes are small, below the half-percent level for the ee CC and µµ CC channels and at the level of 2% for the ee CF channel. The backgrounds from prompt isolated lepton pairs are estimated using simulated samples, as described in Section 3, and consist predominantly of lepton pairs from top-quark processes and from diboson production with a smaller contribution from Z → ττ decays. The other background source arises from events in which at least one of the lepton candidates is not a prompt isolated lepton but rather a lepton from heavy-flavour hadron Table 3: For each of the three channels, yield of events observed in data and expected background yields (multijets, top+electroweak, and total) corresponding to the 2012 data set and an integrated luminosity of 20.3 fb −1 . The uncertainties quoted include both the statistical and systematic components (see text).

Channel Observed
Expected background Multijets (from data) Top+electroweak (from MC) Total ee CC 5.5 × 10 6 6000 ± 3000 13000 ± 3000 19000 ± 4000 µµ CC 7.0 × 10 6 9000 ± 4000 19000 ± 4000 28000 ± 6000 ee CF 1.5 × 10 6 28000 ± 14000 1000 ± 200 29000 ± 14000 decay (beauty or charm) or a fake lepton in the case of electron candidates (these may arise from charged hadrons or from photon conversions within a hadronic jet). This background consists of events containing two such leptons (multijets) or one such lepton (W + jets or top-quark pairs) and is estimated from data using the lepton isolation as a discriminating variable, a procedure described for example in Ref. [50] for electrons. For the central-central channels, the background determination is carried out in the full twodimensional space of (cos θ CS , φ CS ) and in each bin of p T . In the case of the central-forward channel, the multijet background, which is by far the dominant one, is estimated separately for each projection in cos θ CS and φ CS because of the limited amount of data. This is the main reason why the angular coefficients in the central-forward channel are extracted only in projections, as described in Section 1. Figure 4 shows the angular distributions, cos θ CS and φ CS , for the three channels for the data, the Zboson signal MC sample, and the main sources of background discussed above. The total background in the central-central events is below 0.5% and its uncertainty is dominated by the large uncertainty in the multijet background of approximately 50%. The uncertainty in the top+electroweak background is taken conservatively to be 20%. In the case of the central-forward electron pairs, the top+electroweak background is so small compared to the much larger multijet background that it is neglected for simplicity in the fit procedure described in Section 5. Table 3 summarises the observed yields of events in data for each channel, integrated over all values of p T , together with the expected background yields with their total uncertainties from multijet events and from top+electroweak sources. More details of the treatment of the background uncertainties are discussed in Section 6.
There are also signal events that are considered as background to the measurement because they are present in the data only due to the finite resolution of the measurements, which leads to migrations in mass and rapidity. These are denoted "Non-fiducial Z" events and can be divided into four categories: the dominant fraction consists of events that have m Z at the generator level outside the chosen m mass window but pass event selection, while another contribution arises from events that do not belong to the y Z bin considered for the measurement at generator level. The latter contribution is sizeable only in the ee CF channel. Other negligible sources of this type of background arise from events for which the central electron has the wrong assigned charge in the ee CF channel or both central electrons have the wrong assigned charge in the ee CC channel, or for which p Z T at the generator level is larger than 600 GeV. These backgrounds are all included as a small component of the signal MC sample in Fig. 4. Their contributions amount to one percent or less for the ee CC and µµ CC channels, increasing to almost 8% for the ee CF channel because of the much larger migrations in energy measurements in the case of forward electrons. For the 2 < |y Z | < 3.5 bin in the ee CF channel, the y Z migration contributes 2% to the non-fiducial Z background. The fractional contribution of all backgrounds to the total sample is shown explicitly for each channel as a function of p T in Fig. 5 together with the respective contributions of the multijet and CS θ cos   top+electroweak backgrounds. The sum of all these backgrounds is also shown and templates of their angular distributions are used in the fit to extract the angular coefficients, as described in Section 5.

Angular distributions
The measurement of the angular coefficients is performed in fine bins of p Z T and for a fixed dilepton mass window on the same sample as that used to extract from data the small corrections applied to the lepton efficiencies and calibration. The analysis is thus largely insensitive to the shape of the distribution of p Z T , and also to any residual differences in the modelling of the shape of the dilepton mass distribution. It is, however, important to verify qualitatively the level of agreement between data and MC simulation for the cos θ CS and φ CS angular distributions before extracting the results of the measurement. This is shown for the three channels separately in Fig. 6, together with the ratio of the observed data to the sum of predicted events. The data and MC distributions are not normalised to each other, resulting in normalisation differences at the level of a few percent. The measurement of the angular coefficients is, however, independent of the normalisation between data and simulation in each bin of p Z T . The differences in shape in the angular distributions reflect the mismodelling of the angular coefficients in the simulation (see Section 7).

Coefficient measurement methodology
The coefficients are extracted from the data by fitting templates of the P i polynomial terms, defined in Eq. (1), to the reconstructed angular distributions. Each template is normalised by free parameters for its corresponding coefficient A i , as well as an additional common parameter representing the unpolarised cross-section. All of these parameters are defined independently in each bin of p Z T . The polynomial P 8 = 1 + cos 2 θ CS in Eq. (1) is only normalised by the parameter for the unpolarised cross-section.
In the absence of selections for the final-state leptons, the angular distributions in the gauge-boson rest frame are determined by the gauge-boson polarisation. In the presence of selection criteria for the leptons, the distributions are sculpted by kinematic effects, and can no longer be described by the sum of the nine P i polynomials as in Eq. (1). Templates of the P i terms are constructed in a way to account for this, which requires fully simulated signal MC to model the acceptance, efficiency, and migration of events. This process is described in Section 5.1. Section 5.2 then describes the likelihood that is built out of the templates and maximised to obtain the measured coefficients. The methodology for obtaining uncertainties in the measured parameters is also covered there. The procedure for combining multiple channels is covered in Section 5.3, along with alternative coefficient parameterisations used in various tests of measurement results from different channels.

Templates
To build the templates of the P i polynomials, the reference coefficients A ref i for the signal MC sample are first calculated with the moments method, as described in Section 2 and Eq. (5). These are obtained in each of the 23 p Z T bins in Eq. (6), and also in each of the three y Z bins for the y Z -binned measurements. The information about the angular coefficients in the simulation is then available through the corresponding functional form of Eq. (1). Next, the MC event weights are divided by the value of this function on an event-by-event basis. When the MC events are weighted in this way, the angular distributions in the full phase space at the event generator level are flat. Effectively, all information about the Z-boson polarisation is removed from the MC sample, so that further weighting the events by any of the P i terms yields the shape of the polynomial itself, and if selection requirements are applied, this yields the shape of the selection efficiency. The selection requirements, corrections, and event weights mentioned in Section 4 are then applied. Nine separate template histograms for each p Z T and y Z bin j at generator level are finally obtained after weighting by each of the P i terms. The templates t i j are thus three-dimensional distributions in the measured cos θ CS , φ CS , and p T variables, and are constructed for each p Z T and y Z bin. Eight bins in cos θ CS and φ CS are used, while the binning for the reconstructed p T is the same as for the 23 bins defined in Eq. (6). By construction, the sum of all signal templates normalised by their reference coefficients and unpolarised cross-sections agrees exactly with the three-dimensional reconstructed distribution expected for signal MC events. Examples of templates projected onto each of the dimensions cos θ CS and φ CS for the y Z -integrated ee CC channel in three illustrative p Z T ranges, along with their corresponding polynomial shapes, are shown in Fig. 7. The polynomials P 1 and P 6 are not shown as they integrate to zero in the full phase space in either projection (see Section 5.2). The effect of the acceptance on the polynomial shape depends on p Z T because of the event selection, as can be seen from the difference between the template polynomial shapes in each corresponding p Z T bin. This is particularly visible in the P 8 polynomial, which is uniform in φ CS , and therefore reflects exactly the acceptance shape in the templated polynomials. In Appendix B, two-dimensional versions of Fig. 7 are given for all nine polynomials in Figs. 21 -23. These two-dimensional views are required for P 1 and P 6 , as discussed above.
Templates T B are also built for each of the multijet, top+electroweak, and non-fiducial Z-boson backgrounds discussed in Section 4.2. These are normalised by their respective cross-sections times luminosity, or data-driven estimates in the case of the multijet background. The templates for the projection measurements in the ee CF channel are integrated over either the cos θ CS or φ CS axis at the end of the process.
Templates corresponding to variations of the systematic uncertainties in the detector response as well as in the theoretical modelling are built in the same way, after varying the relevant source of systematic uncertainty by ±1 standard deviation (σ). If such a variation changes the A ref i coefficients in the MC prediction, for example in the case of PDF or parton shower uncertainties, the varied A ref i coefficients are used as such in the weighting procedure. In this way, the theoretical uncertainties on the predictions are not directly propagated to the uncertainties on the measured A i coefficients. However, they may affect indirectly the measurements through their impact on the acceptance, selection efficiency, and migration modelling.

Likelihood
A likelihood is built from the nominal templates and the varied templates reflecting the systematic uncertainties. A set of nuisance parameters (NPs) θ = {β, γ} is used to interpolate between them. These are constrained by auxiliary probability density functions and come in two categories: β and γ. The first category β are the NPs representing experimental and theoretical uncertainties. Each β m in the set β = β 1 , ..., β M are constrained by unit Gaussian probability density functions G(0|β m , 1) and linearly interpolate between the nominal and varied templates. These are defined to have a nominal value of zero, with β m = ±1 corresponding to ±1σ for the systematic uncertainty under consideration. The total number of β m is M = 171 for the ee CC + µµ CC channel and M = 105 for the ee CF channel. The second category γ are NPs that handle systematic uncertainties from the limited size of the MC samples. For The p T dimension that normally enters through migrations is also integrated over. The differences between the polynomials and the templates reflect the acceptance shape after event selection. each bin n in the reconstructed cos θ CS , φ CS , and p T distribution, γ n in the set γ = γ 1 , ..., γ N bins , where N bins = 8 × 8 × 23 is the total number of bins in the reconstructed distribution, has a nominal value of one and normalises the expected events in bin n of the templates. They are constrained by Poisson probability density functions P(N n eff |γ n N n eff ), where N n eff is the effective number of MC events in bin n. The meaning of "effective" here refers to corrections applied for non-uniform event weights. When all signal and background templates are summed over with their respective normalisations, the expected events N n exp in each bin n can be written as: where: The summation over the index j takes into account the contribution of all p Z T bins at generator level in each reconstructed p T bin. This is necessary to account for migrations in p T . The likelihood is the product of Poisson probabilities across all N bins bins and of auxiliary constraints for each nuisance parameter β m : Unlike in the ee CC and µµ CC channels that use both angular variables simultaneously, the ee CF measurements are performed in projections (see Eq. (2) and Eq. (3)), and therefore the A 1 and A 6 coefficients are not measured in this channel. The P i polynomials that normally integrate to zero when projecting onto one angular variable in full phase space may, however, not integrate to zero if their shape is distorted by the event selection. The residual shape is not sufficient to properly constrain their corresponding A i , and therefore an external constraint is applied to them. For the A i that are largely independent of y Z (A 0 and A 2 ), the constraints are taken from the independent y Z -integrated measurements in the combined ee CC + µµ CC channel. For the y Z -dependent coefficients A 1 , A 3 , and A 4 , which are inaccessible to the ee CC + µµ CC channels in the y Z bin in which ee CF is used, predictions from Powheg + MiNLO [43] are used.
The migration of events between p T bins leads to anti-correlations between A i in neighbouring bins which enhance the effects of statistical fluctuations. To mitigate this effect and aid in resolving underlying structure in the A i spectra, the A i spectra are regularised by multiplying the unregularised likelihood by a Gaussian penalty term, which is a function of the significance of higher-order derivatives of the A i with respect to p Z T . The covariance terms between the A i j coefficients are taken into account and computed first with the unregularised likelihood. This has parallels with, for example, regularised Bayesian unfolding, where additional information is added through the prior probability of unfolded parameter values [54,55]. As is the case there, the choice of penalty term (or in the Bayesian case, the choice of added information) must be one that leads to a sound result with minimal bias. See Appendix C for more details.
The uncertainties in the parameters are obtained through a likelihood scan. For each parameter of interest A i j , a likelihood ratio is constructed as In the denominator, the likelihood is maximised unconditionally across all parameters of interest and NPs.
In the numerator, the likelihood is maximised for a specific value of a single A i j . The maximum likelihood estimators for the other parameters of interestÂ and NPsθ are in general a function of A i j , hence the explicit dependence is shown in the numerator. The Minuit package is used to perform numerical minimisation [56] of −2 log Λ(A i j ), and a two-sided test statistic is built from the likelihood ratio: This is asymptotically distributed as a χ 2 with one degree of freedom [57]. In this case, the ±1σ confidence interval of A i j is defined by the condition

Combinations and alternative parameterisations
When applicable, multiple channels are combined through a simple likelihood multiplication. Each likelihood can be decomposed into three types of terms: those that contain the observed data in each channel, denoted L i (A, σ, θ), the auxiliary terms that constrain the nuisance parameters θ, denoted A i (θ i ), and the auxiliary term that imposes the regularisation, A reg (A). There are a total of M cb NPs, corresponding to the total number of unique NPs, including the total number of bins, across all combined channels. With this notation the combined likelihood can be written as: There are several instances in which a combination of two channels is performed. Within these combinations, the compatibility of the channels is assessed. The measurements in the first two y Z bins and the y Z -integrated configuration are obtained from a combination of the ee CC and µµ CC channels. The y Zintegrated µµ CC and ee CF channels are also combined in order to assess the compatibility of the high y Z region probed by the ee CF channel and the lower rapidity region probed by the central-central channels.
The compatibility of channels is assessed through a reparameterisation of the likelihood into parameters that represent the difference between the coefficients in two different channels. For coefficients A a i j and A b i j in respective channels a and b, difference parameters ∆A i j ≡ A a i j − A b i j are defined that effectively represent the difference between the measured coefficients in the two channels. Substitutions are made in the form of A a i j → ∆A i j + A b i j . These new parameters are measured with the same methodology as described in Section 5.2. Similar reparameterisations are also done to measure the difference between the A 0 and A 2 coefficients. These reparameterisations have the advantage that the correlations between the new parameters are automatically taken into account.

Measurement uncertainties
Several sources of statistical and systematic uncertainty play a role in the precision of the measurements presented in this paper. In particular, some of the systematic uncertainties impact the template-building procedure described in Section 5.1. For this reason, templates are rebuilt after each variation accounting for a systematic uncertainty, and the difference in shape between the varied and nominal templates is used to evaluate the resulting uncertainty.
A description of the expected statistical uncertainties (both in data in Section 6.1 and in simulation in Section 6.2) and systematic uncertainties (experimental in Section 6.3, theoretical in Section 6.4, and those related to the methodology in Section 6.5) associated with the measurement of the A i coefficients is given in this section. These uncertainties are summarised in Section 6.6 in three illustrative p Z T bins for the ee CC , µµ CC (and their combination), and ee CF channels. The evolution of the uncertainty breakdown as a function of p Z T is illustrated there as well.

Uncertainties from data sample size
Although the harmonic polynomials are completely orthogonal in the full phase space, resolution and acceptance effects lead to some non-zero correlation between them. Furthermore, the angular distributions in a bin of reconstructed p T have contributions spanning several generator-level p Z T bins. This leads to correlations between the measured coefficients which increase their statistical uncertainties. The amount of available data is the largest source of uncertainty, although the resolution and binning in the angular variables also play a role. A discussion of the categorisation of this uncertainty may be found in Appendix D.

Uncertainties from Monte Carlo sample size
Statistical uncertainties from the simulated MC samples are treated as uncorrelated between each bin of the three-dimensional (p T , cos θ CS , φ CS ) distribution. Although the events used to build each template are the same, they receive a different weight from the different polynomials, and are therefore only partially correlated. It was verified that assuming that the templates are fully correlated yields slightly more conservative uncertainties, but central values identical to those obtained using the fully correct treatment. For simplicity, this assumption is used for this uncertainty.

Experimental systematic uncertainties
Experimental systematic uncertainties affect the migration and efficiency model of the detector simulation, impacting the variables used to construct the templates and the event weights applied to the simulation.
Lepton-related systematic uncertainties: Scale factors correcting the lepton reconstruction, identification, and trigger efficiencies to those observed in data [50][51][52] are applied to the simulation as event weights. The statistical uncertainties of the scale factors tend to be naturally uncorrelated in the kinematic bins in which they are measured, while the systematic uncertainties tend to be correlated across these bins. Lepton calibration (electron energy scale and resolution as well as muon momentum scale and resolution) [52,53] and their associated uncertainties are derived from data-driven methods and applied to the simulation. Whereas the charge misidentification rate for muons is negligible, the probability for the electron charge to be misidentified can be significant for central electrons at high |η|, due to bremsstrahlung in the inner detector and the subsequent conversion of the photon. This uncertainty is a potential issue in particular for the ee CF channel, where the measured charge of the central electron sets the charge of the forward electron (where no charge determination is possible). Measurements of the per-electron charge misidentification rate using same-charge electron pairs have been done in data and compared to that in simulation; the systematic uncertainty coming from this correction has a negligible impact on the measurement.
Background-related systematic uncertainties: Uncertainties in the multijet background estimate come from two sources. The first source is the statistical uncertainty in the normalisation of the background in each bin of reconstructed p T . The second is the systematic uncertainty of the overall background normalisation, which is evaluated using alternative criteria to define the multijet background templates. These uncertainties are applied to all three channels and treated as uncorrelated amongst them. In addition, a 20% systematic uncertainty uncorrelated across p T bins but correlated across the ee CC and µµ CC channels is applied to the estimation of the top+electroweak background.
Other experimental systematic uncertainties: Several other potential sources of experimental systematic uncertainty are considered, such as event pileup or possible detector misalignments which might affect the muon momentum measurement, and are found to contribute negligibly to the overall measurement uncertainties. The uncertainty in the integrated luminosity is ± 2.8%. It is derived following the same methodology as that detailed in Ref. [58]. It only enters (negligibly) in the scaling of the background contributions evaluated from the Monte Carlo samples.

Theoretical systematic uncertainties
Theoretical systematic uncertainties due to QCD renormalisation/factorisation scale, PDFs, parton-shower modelling, generator modelling, and QED/EW corrections are considered. They are evaluated using either event weights, for example through PDF reweighting, or templates built from alternative Monte Carlo samples. The templates built after each variation accounting for a systematic uncertainty have their own set of reference coefficients so that each variation starts from isotropic angular distributions. This procedure is done so that uncertainties in the simulation predictions for the coefficients propagate minimally to the uncertainties in the measurement; rather, uncertainties in the measurement are due to the theoretical uncertainty of the migration and acceptance modelling. A small fraction of the acceptance can, however, be attributed to the behaviour of coefficients outside the accessible y Z range. In this specific case, the theoretical predictions used for the coefficients can have a small influence on the uncertainties in the measured coefficients.
QCD scale: These systematic uncertainties only affect the predictions over the small region of phase space where no measurements are available. They are evaluated by varying the factorisation and renormalisation scale of the predicted coefficients in the region |y Z | > 3.5 (see Fig. 3). The changes induced in the templates due to the variation in acceptance are used to assess the impact of this uncertainty, which is found to be negligible.
PDF: These systematic uncertainties are computed with the 52 CT10 eigenvectors representing 26 independent sources. The CT10 uncertainties are provided at 90% CL, and are therefore rescaled by a factor of 1.64 to bring them to 68% CL variations. Events are also reweighted using the NNPDF2.3 [59] and MSTW [60] PDFs and are treated as independent systematics. These two-point variations are symmetrised in the procedure.
Parton showers: The Powheg + Herwig samples are used to compute an alternative set of templates. The shape difference between these and the templates obtained from the baseline Powheg + Pythia 8 samples are used to evaluate the underlying event and parton shower uncertainty.
Event generator: These systematic uncertainties are evaluated through the reweighting of the rapidity distribution of the nominal Powheg + Pythia 8 MC sample to that from Sherpa, which corresponds approximately to a 5% slope per unit of |y Z |. An alternative set of signal templates is built from this variation, using the new set of reference coefficients averaged over rapidity after the reweighting.
QED/EW corrections: The impact of the QED FSR corrections on the measurements is accounted for by the uncertainties in the lepton efficiencies and scales. The contribution of the EW corrections to the calculation of the Z-boson decay angular distributions in Eq. (1) is estimated to be negligible around the Z-pole mass, based on the extensive and detailed work done at LEP in this area [6,7,61], as discussed in Section 1.
The PDF uncertainties were found to be the only non-negligible source of theoretical systematic uncertainty in the measured A i coefficients, and are in particular the dominant source of uncertainty in the measurement of A 0 at low p Z T .

Systematic uncertainties related to the methodology
Systematic uncertainties related to the template building, fitting, and regularisation methodology are considered. These could manifest through sensitivity to the p Z T shape in the simulation, the particular shape of the A i coefficient being fitted, or possible biases caused by the regularisation scheme. p T shape: The sensitivity to the shape of the p T spectrum is tested in two different ways. First, the shape of the p T spectrum in simulation is reweighted with a polynomial function so that it approximately reproduces the reconstructed spectrum in data. The impact of this procedure is expected to be small, since the signal is normalised to the data in fine bins of p Z T . Second, the p T shape within each p T bin is reweighted to that of the data. Since the binning is fine enough that the p T shape does not vary too rapidly within one bin, the impact of this is also small. Regularisation: The potential bias induced by the regularisation is evaluated with pseudo-experiments. A sixth-order polynomial is fit to the PowhegBox + Pythia 8 set of A ref i coefficients to obtain a continuous reference spectrum y i j . Pseudo-data are randomised around the expected distribution obtained from this fit using a Poisson distribution for each bin. The difference between y i j and the expectation value E[A i j ] of the distribution of fitted and regularised A i j is an estimate of the potential bias in the regularised A i j . The envelope of |E[A i j ] − y i j | is symmetrised and taken to be the bias uncertainty. (See Appendix C for more details.) The effect of p Z T reweighting and closure of A i spectra were found to be negligible. The only nonnegligible source of uncertainty in the methodology was found to be the regularisation bias, which can have a size approaching the statistical uncertainty of A 0 and A 2 . Tables 4-7 show the uncertainties in each measured coefficient in three representative p Z T bins, along with the impact of each category of systematics. The theoretical uncertainties are dominated by the PDF uncertainties, which in a few cases are larger than the statistical uncertainties. The experimental uncertainties are dominated by the lepton uncertainties and are the leading source of systematic uncertainty for low values of p Z T for the A 2 coefficient. The large uncertainties assigned to the multijet background estimates and their shape have a negligible impact on this measurement.

Summary of uncertainties
The dominant uncertainty in the measurements of the A i coefficients is in most cases the statistical uncertainty, even in the most populated bins at low p Z T , which contain hundreds of thousands of events. The exception is the A 0 coefficient where PDF and electron efficiency uncertainties dominate for p Z T values below 80 GeV. The next largest uncertainty is due to the signal MC statistical uncertainty. This is reflected in Fig. 8, which shows the uncertainty evolution versus p Z T for A 0 − A 2 , including a breakdown of the systematic uncertainties for both the unregularised and regularised measurements. The evolution versus p Z T of the total, statistical, and systematic uncertainties is shown for all other coefficients in Fig. 9 for the regularised measurement.

Results
This section presents the full set of experimental results. The compatibility between channels is assessed in Section 7.1. The measured A i coefficients are then shown in Section 7.2. A test is also performed to check for non-zero values of the A 5,6,7 coefficients. Several cross-checks are presented in Section 7.3, including a test of the validity of the nine-term decomposition, probing for the presence of higher-order P i polynomial terms.

Compatibility between channels
Given that a complex fitting procedure based on reconstructed observables is used, the compatibility between different channels is assessed with a strict quantitative test. The likelihood is parameterised in terms of ∆A i j ≡ A µµ i j − A ee i j for coefficient index i and p Z T bin j, as described in Section 5.3. The compatibility of the ∆A i j with zero can be quantified via a χ 2 test taking into account all systematic Table 4: Summary of regularised uncertainties expected for A 0 , A 2 , and A 0 − A 2 at low (5-8 GeV), mid , and high (132-173 GeV) p Z T for the y Z -integrated configuration. The total systematic uncertainty is shown with the breakdown into its underlying components. Entries marked with "-" indicate that the uncertainty is below 0.0001.    Table 5: Summary of regularised uncertainties expected for A 1 , A 3 and A 4 at low (5-8 GeV), mid (22-25.5 GeV), and high (132-173 GeV) p Z T for the y Z -integrated configuration. The total systematic uncertainty is shown with the breakdown into its underlying components. Entries marked with "-" indicate that the uncertainty is below 0.0001.   Table 6: Summary of regularised uncertainties expected for A 5 , A 6 and A 7 at low (5-8 GeV), mid (22-25.5 GeV), and high (132-173 GeV) p Z T for the y Z -integrated configuration. The total systematic uncertainty is shown with the breakdown into its underlying components. Entries marked with "-" indicate that the uncertainty is below 0.0001.   Table 7: Summary of regularised uncertainties expected for the coefficients at low (5-8 GeV), mid (22-25.5 GeV), and high (73.4-85.4 GeV) p Z T for the 2 < |y Z | < 3.5 configuration. The total systematic uncertainty is shown with the breakdown into its underlying components. Entries marked with "-" indicate that the uncertainty is below 0.0001.    uncertainty correlations. The χ 2 values are first computed for each coefficient i and across all p Z T bins j, then for all coefficients and p Z T bins simultaneously. This test is done in the y Z -integrated case for the differences between the measurements extracted from the µµ CC and ee CC events and from the µµ CC and ee CF events, as well as in the first two y Z bins for the µµ CC and ee CC events. The χ 2 values are tabulated in Table 8 and indicate almost all the differences are compatible with zero.
The ∆A i j spectra are shown in Fig. 10 for the y Z -integrated ee CC and µµ CC channels. The regularised and unregularised spectra are overlayed. Visually, it appears that these results are compatible with zero. In some cases, the unregularised ∆A i j show alternating fluctuations above and below zero due to anticorrelations between neighbouring p Z T bins. These are smoothed out in the regularised results, which come at the expense of larger bin-to-bin correlations.

Results in the individual and combined channels
The measurements represent the full set of y Z -integrated coefficients, including the difference A 0 − A 2 , as a function of p Z T , as well as the y Z -dependent coefficients as a function of p Z T in the available y Z bins. The combination of the ee CC and µµ CC channels is used for the y Z -integrated measurements and the measurements in the first two y Z bins, while the ee CF channel is used for the measurements in the  Figure 9: The total uncertainty as a function of p Z T along with a breakdown into statistical and systematic components for all coefficients in the regularised y Z -integrated ee CC + µµ CC measurement. The full (open) circles represent the measured differences before (after) regularisation. The error bars represent the total uncertainty in the measurements. Table 8: Tabulation of the compatibility of the measured ∆A i with zero reported as χ 2 per degree of freedom (N DoF ), where ∆A i represents the difference between the A i coefficient extracted from the µµ CC and ee CC events (left) and from the µµ CC and ee CF events (right). For the ee CC versus µµ CC tests, the number of degrees of freedom is 23 for the tests of the individual coefficients and 184 for the tests of all coefficients simultaneously. Likewise, for the ee CF versus µµ CC tests, there are 19 degrees of freedom for the tests of the individual coefficients, 38 for the simultaneous test in the cos θ CS projection, and 76 for the simultaneous test in the φ CS projection. The comparisons are not performed for the A 1 and A 6 coefficients between the µµ CC and ee CF channels (see Section 5.2). last y Z bin. A summary of these measurements is tabulated in Tables 9-10 for three representative p Z T bins. Figure 11 shows the y Z -integrated measurements for all A i and overlays of the y Z -dependent A i in each accessible y Z bin. The A 1 and A 6 measurements are missing from the third y Z bin since they are inaccessible in the projections used in the ee CF channel (see Section 5.2). Also, a measurement of A 0 − A 2 is missing in this bin since A 0 and A 2 are accessible in different projections. Complete tables can be found in Appendix F along with additional figures in y Z bins. Similarly to the regularised ∆A i j measurements, there is a large degree of correlation from bin to bin. This, coupled with statistical fluctuations, can lead to correlated deviations in the spectra, for example near p Z T = 40 GeV for A 4 in the 2 < |y Z | < 3.5 bin, and for A 1 in the 0 < |y Z | < 1 bin. Visually, the coefficients A 5,6,7 all show a trend towards non-zero positive values in the region with p Z T around 100 GeV.

Cross-checks
Several cross-checks are performed to ensure that the fit is of good quality and that the underlying theoretical assumptions are valid to the extent of the precision of the analysis.
The signal MC distributions are reweighted to the full set of measured parameters. An event-by-event weight is calculated as a ratio using the right-hand side of Eq. (1): the numerator uses the measured parameters, and the denominator uses the reference values in the MC simulation. Distributions are obtained after applying this reweighting; the cos θ CS and φ CS distributions integrated in y Z are shown in Fig. 12, along with their bin-by-bin pulls, obtained by combining the statistical and systematic uncertainties. Overall, the data and MC simulation agree well. One observes significant pulls near cos θ CS = 0 for the ee CF channel, but the number of events in this region is very small and its impact on the coefficient measurements is negligible.
Finally, the degree to which the data follow the nine-P i polynomial decomposition is tested by checking for the presence of higher-order P i in the data. The original nine P i are up to second order in spherical harmonics. The template-building methodology described in Section 5.1 is extended to have more than nine P i by using third-and fourth-order spherical harmonics, corresponding to 16 additional P i . One additional P i template is fitted at a time. The higher-order coefficients are found to be compatible with zero using a χ 2 test as in Section 7.1, leading to the conclusion that any possible breaking of the nine P i polynomial decomposition is beyond the sensitivity of the analysis.

Comparisons with theory
In this section, the measurements are compared to the most precise fixed-order calculations currently available. They probe the dynamics of perturbative QCD, including the presence of higher-order corrections, and explore the effects from the V − A structure of Z-boson couplings. These comparisons are made with both the y Z -integrated and y Z -binned measurements. For the y Z -integrated measurements and for the 0 < |y Z | < 1 and 1 < |y Z | < 2 bins, the combined ee CC and µµ CC measurements are used, while the ee CF measurements are used for the 2 < |y Z | < 3.5 bin. In all cases, the regularised uncertainties described Table 10: Summary of the measured coefficients in the ee CC + µµ CC channel for the two bins 0 < |y Z | < 1 and 1 < |y Z | < 2 and in the ee CF channel for the 2 < |y Z | < 3.5 bin at low (5-8 GeV), mid (22-25.5 GeV), and high (132-173 GeV) p Z T . The uncertainties are given as ±(stat.) ± (syst.). 0.022 ± 0.010 ± 0.006 0.071 ± 0.013 ± 0.007  in Section 5 are used for the data. The measurements are also compared to various event generators, in particular to probe different parton-shower models and event-generator implementations.
The overlays of the y Z -integrated measurements are shown in Figs. 13-15 for all coefficients. The calculations from DYNNLO are shown at NNLO for p Z T > 2.5 GeV with their uncertainties computed as a sum in quadrature of statistical, QCD scale, and PDF uncertainties, as described in Section 2. The Powheg + MiNLO predictions, which are shown only including statistical uncertainties, were obtained using the Z + jet process at NLO [43] over the full p Z T range. Owing to numerical issues in the phase-space integration, the Powheg + MiNLO results show fluctuations beyond their statistical uncertainties. The formal accuracy of both calculations is the same, namely O(α s ) for the predictions of the A i coefficients as a function of p Z T . The left-hand plots in these figures illustrate the behaviour of each coefficient as a function of p Z T , while the right-hand plots, in which the data measurements are used as a reference, show to which extent the various theoretical predictions agree with the data. In the very first p Z T bin, φ CS has poor resolution and therefore suffers from larger measurement uncertainties. This is reflected in the deviation from the prediction in A 2 , for example, which is derived primarily from φ CS .
The predictions from the DYNNLO and Powheg + MiNLO calculations agree with the data within uncertainties for most coefficients. The striking exception is the A 2 coefficient, which rises more slowly as p Z T increases in the data than in the calculations. The data confirm that the Lam-Tung relation (A 0 − A 2 = 0) does not hold at O(α 2 s ). For p Z T > 50 GeV, significant deviations from zero, almost a factor of two larger than those predicted by the calculations, are observed. Since the impact of the PDF uncertainties on the calculations is very small, these deviations must be due to higher-order QCD effects.
In the case of the A 5,6,7 coefficients, the trend towards non-zero values at high p Z T discussed in Section 7 is also compatible with that from the predictions, although it is at the limit of the sensitivity for both the data and the calculations. As shown in Fig. 15 and also in Table 1, the predictions from DYNNLO suggest that the values of the A 5,6,7 coefficients should be at the level of 0.005 at high values of p Z T . A test is performed to quantify the significance of the deviation from zero (see Appendix E). A signed χ 2 test statistic is defined based on the tail probability of each individual measurement, taking into account the correlations between the parameters in bins of p Z T . An ensemble test is performed to compute the observed and expected significance of all three coefficients together, where pseudo-data from DYNNLO is used for the expected value. This test gives an observed (expected) significance of 3.0 (3.2) standard deviations.
The measurements of the A 1 , A 3 , and A 4 coefficients in the three y Z bins (only the first two bins are available for the A 1 coefficient) are compared to the predictions in Figs. 16-18. Overall, the predictions and the data agree for all three y Z bins. These coefficients are the only ones that display any significant y Z dependence and it is interesting to note that, for high values of p Z T , the A 1 and A 3 coefficients increase as y Z increases. As explained in Section 1 and detailed in Appendix A, at low values of p Z T , the measured value of the A 4 coefficient can be directly related to the Weinberg angle sin 2 θ W [62]. The strong dependence of the value of the A 4 coefficient on |y Z | is, however, mostly a consequence of the approximation made for the interacting quark direction in the CS reference frame on an event-by-event basis. The impact of this approximation decreases at higher values of |y Z |, and, as a result, the measured and expected values of the A 4 coefficient increase, as can be seen in Figs. 16 -18. The effect of the parton-shower modelling and matching scheme on the reference angular coefficients is explored in Fig. 19, which shows a comparison of the measurements of A 0 , A 1 , A 2 , and A 0 − A 2 with DYNNLO at NLO and NNLO, PowhegBox (without parton shower), and with the same process in PowhegBox with the parton shower simulated with Pythia 8 (PowhegBox + Pythia 8) and Herwig (PowhegBox + Herwig). The predictions from DYNNLO at NLO and PowhegBox without parton shower, which are formally at the same level of accuracy, agree for A 1 and A 2 . For the A 2 coefficient, which is the most sensitive one to higher-order corrections, adding the parton-shower simulation to the Powheg-Box Z-boson production process brings the predictions closer to DYNNLO at NNLO. This is consistent with the assumption that the parton-shower model emulates higher-order effects, although the discrepancy between the measurements and the parton-shower models is larger than that with DYNNLO at NNLO. The A 0 coefficient has an unexpected offset of −0.025 at low values of p Z T in the PowhegBox implementation. This effect is also reflected in the predictions for A 0 − A 2 and has been corrected in the more recent version of PowhegBox (v2.1) used in this paper for the Z + jet predictions with Powheg + MiNLO [14][15][16][17]. The predictions from DYNNLO at NLO and NNLO agree well with the data measurements for the A 0 coefficient, but overestimate the rise of the A 2 coefficient at higher values of p Z T , as discussed above. Finally, it is interesting to note that, whereas the agreement between Pythia 8 and Herwig is good for most of the coefficients, the A 1 coefficient displays significant differences between the two predictions over most of the p Z T range. Although this might be ascribed to differences between the parton-shower model and matching schemes at intermediate values of p Z T , it is somewhat surprising to observe large differences for the highest values of p Z T . Figure 20 shows a comparison of the measurements of A 0 , A 2 and A 0 − A 2 with Sherpa 1.4 (up to five jets at LO) and Sherpa 2.1 [39][40][41][42]. The effect of simulating Sherpa 2.1 events (up to two jets at NLO and up to five jets at LO for higher jet multiplicities) is explicitly shown. None of the configurations correctly predict the behaviour of A 0 or A 2 . The Sherpa 2.1 version follows the data more closely than the Sherpa 1.4 version. In addition, in all versions except Sherpa 2.1 with two jets, significant higher-order polynomial behaviour was found to be present. This is probably due to the matrix-element matching scheme used in the event generator for the calculation of the Z + n-jet process for n > 2.            Figure 19: Distributions of the angular coefficients A 0 , A 2 , A 0 − A 2 and A 1 (from top to bottom) as a function of p Z T . The results from the y Z -integrated measurements are compared to the DYNNLO predictions at NLO and at NNLO, as well as to those from PowhegBox + Pythia8 and PowhegBox + Herwig (left). The differences between the calculations and the data are also shown (right), with the shaded band around zero representing the total uncertainty in the measurements. The error bars for the calculations show the total uncertainty for DYNNLO, but only the statistical uncertainties for PowhegBox.

Summary
This paper presents a precise set of measurements of the Z-boson production dynamics in the Z-boson pole region, through the angular distributions of the leptons. The data analysed correspond to 20.3 fb −1 of pp collisions at √ s = 8 TeV, collected by the ATLAS detector at the CERN LHC. The measurements are obtained as a function of p Z T , integrated over y Z and in bins of y Z , covering almost the full range of y Z spanned by Z-boson production at √ s = 8 TeV. This is made possible by exploiting the decomposition of the production cross-section into nine terms, where in each term the angular coefficients that encapsulate the production dynamics are factorised from the decay dynamics described by angular polynomials. Templates of the nine polynomials folded to the detector level are fitted to the data to extract the angular coefficients in the full phase space of the Z boson.
Over most of the phase space, the measurements that are obtained from samples of electron and muon pairs covering respectively the ranges 0 < |y Z | < 3.5 and 0 < |y Z | < 2.5 are limited only by statistical uncertainties in the data. These uncertainties are small and range from 0.002 at low p Z T to 0.008 at p Z T = 150 GeV. The experimental systematic uncertainties are much smaller in almost all cases. The theory systematic uncertainties are minimised through the template-building procedure, such that the PDF uncertainties, which are the dominant source of theoretical uncertainties, are below 0.004 in all cases.
The measurements, when compared to theoretical calculations and to predictions from MC generators, are precise enough to probe QCD corrections beyond the formal accuracy of the calculations. A significant deviation from the O(α 2 s ) predictions from DYNNLO is observed for A 0 − A 2 , indicating that higherorder QCD corrections are required to describe the data. Evidence at the 3σ level is found for non-zero A 5,6,7 coefficients, consistent with expectations from DYNNLO at O(α 2 s ). The measurements also provide discrimination between various event generators, in particular in terms of the related implementation of different parton-shower models.
The measurements of the A i coefficients, in particular through the correlation of the angular distributions with the lepton transverse momentum distributions, are thus an important ingredient to the next steps in precision measurements of electroweak parameters at the LHC, both for the effective weak mixing angle sin 2 θ W and for the W-boson mass.

Acknowledgements
We thank CERN for the very successful operation of the LHC, as well as the support staff from our institutions without whom ATLAS could not be operated efficiently.

Appendix
A. Theoretical formalism Following the notation in Ref. [3], the lepton-hadron correlations in the pp → Z → process are described by the contraction of the lepton tensor L µν with the parton-level hadron tensor H µν . The tensor L µν acts as an analyser of the structure of H µν , which carries effective information about the polarisation of the Z boson mediating the interaction. The angular dependence can be elucidated by introducing nine helicity density matrix elements where m, m = +, 0, − and are the polarisation vectors for the Z boson, defined with respect to its rest frame. The angular dependence of the differential cross-section can be written as: dσ α dp Z T dy Z dm Z , M = {U + L, L, T, I, P, A, 7, 8, 9}, (14) where the g α (θ, φ) are second-order harmonic polynomials, multiplied by normalisation constants. The helicity cross-sections σ α are linear combinations of the helicity density matrix elements H mm : The unpolarised cross-section is denoted historically by σ U+L , whereas σ L,T,I,P,A,7,8,9 characterise the Zboson polarisation. Respectively, these are the contributions to the Z-boson cross-section from longitudinally and transversely polarised states, transverse-longitudinal interference, etc., as described in Ref. [2].
The individual helicity cross-sections depend on the coupling coefficients of the Z boson as follows: where v q (v ) and a q (a ) denote the vector and axial-vector coupling of the Z boson to the quarks (leptons). The cross-sections σ U+L,L,T,I,9 receive contributions from the parity-conserving component of the hadron tensor, while the others, σ P,A,7,8 , are proportional to the parity-violating component of H µν . However, the angular polynomials g P,A,9 (θ, φ) are parity-violating as well, so contributions to the angular distributions from σ U+L,L,T,I,P,A are parity-conserving.
It is standard notation to factorise out the unpolarised cross-section, and to present the five-dimensional differential cross-section as an expansion in harmonic polynomials P i (cos θ, φ) and dimensionless angular coefficients A 0−7 , which represent ratios of helicity cross-sections with respect to the unpolarised one, as follows: This leads to Eq. (1), as discussed in Section 1.

B. Additional Templates
To expand upon the one-dimensional templates shown in Section 5, the two-dimensional versions are shown here. The dimension corresponding to migrations in p T is integrated over. Figures 21-23 show each analytical polynomial for each corresponding coefficient along with the templated versions in three representative p Z T bins after acceptance and selection requirements. The differences between the analytical polynomials and their templates reflect primarily the effect of the acceptance shape in the angular variables, and to a lesser extent resolution effects.

C. Regularisation
The migration of events between p Z T bins leads to anti-correlations between the measured A i in neighbouring p Z T bins which enhance the effects of statistical fluctuations. To mitigate this effect and aid in resolving the underlying structure of the A i spectra, the A i coefficients are regularised by imposing a Gaussian penalty term on the significance of their higher-order derivatives with respect to p Z T . This penalty multiplies the likelihood in Eq. (8).
The exact derivative order is chosen based on the expected reduction in statistical uncertainties of the measurement and the potential bias that the regularisation scheme may introduce. The smaller statistical uncertainties come with increased positive correlation between neighbouring coefficients. The n th derivative of A i j is defined as: where A (0) i j ≡ A i j . The derivatives are staggered between even and odd orders in order to create a derivative definition more symmetric around each p Z T bin. Since the measurement is determined from the exact likelihood expression for the coefficients, the covariance matrix Σ of the coefficients can be derived based on the second-order partial derivatives of the likelihood [57]. Along with the uncertainties in A i and A j , namely σ(A i ) and σ(A j ), their correlation can be computed as ρ i j = Σ i j /σ(A i )σ(A j ). This is done based on pseudo-data taken from Powheg + Pythia 8 and is shown in Fig. 24 before and after regularisation.
A Jacobian matrix J (and its transpose J T ) is used to transform the covariance matrix of the coefficients to the covariance matrix of their derivatives. A regularisation strength γ is introduced to control the amount by which the derivatives are penalised. The penalty term applied to the likelihood that controls the regularisation is therefore defined as In all channels, a regularisation scheme using sixth-order derivatives is used.
In the limit that the regularisation procedure described above has infinite strength, an n th -order derivative regularisation fixes the measured spectrum to be an (n − 1) th -order polynomial. This can be seen in Fig. 25, which shows the residual of a fifth-order polynomial fit to the A 0 spectrum regularised with sixth-order derivatives and strength γ = 100; the fit is nearly perfect (the regularisation strength used in this case is large but not infinite and so there are some small non-zero residuals). Also shown is the residual of a fourth-order polynomial fit to this same spectrum; a fifth-order term can be clearly observed in the residual. The regularisation bias B[A i j ] in the coefficients is evaluated using pseudo-experiments based on the difference between the expectation value of the best-fit coefficient E[A i j ] and the value of the coefficient y i j used to randomise the data: The choice of y i j is derived from a sixth-order polynomial fit to the Powheg + Pythia 8 reference coefficients.
The derived uncertainty due to the regularisation bias in the y Z -integrated A 0 coefficient in the ee CC +µµ CC channel is shown in Fig. 26 for four different regularisation strengths, along with the corresponding statistical uncertainty of the coefficient for each strength. As can be seen, the regularisation uncertainty increases with increasing regularisation strength, while the corresponding statistical uncertainty decreases, as expected. In the limit that the regularisation strength goes to zero, the statistical uncertainty approaches the unregularised one. Along with the decrease in statistical uncertainty comes an increase in correlation among the measured coefficients of neighbouring p Z T bins. The regularisation bias uncertainty appears to plateau between γ = 10 and γ = 100, which corresponds to the limit that the spectrum is fixed to a sixth-order polynomial, as described above. Based on these studies, a strength of γ = 100 is chosen for the ee CC and µµ CC channels, while the scheme in the ee CF channel is based on a strength of γ = 5.
Overlays of the regularised measurements for the ee CC + µµ CC channel in the y Z -integrated configuration are shown in Fig. 27 for A 0−7 and in Fig. 28 for A 0 − A 2 . In the unregularised results, there are many binto-bin fluctuations that enter primarily through anti-correlations between neighbouring measurements. In contrast, the regularised results are largely correlated from bin-to-bin and are much smoother.   ← Parameter Figure 29: Categorisation of parameters leading to the data statistical uncertainty in the measured coefficients illustrated by the uncertainty categorisation for A 0 in p Z T bin 0.

D. Categorisation of statistical uncertainties
The categorisation of statistical uncertainties is illustrated in Fig. 29 for A 0 in p Z T bin 0. Uncertainties due to the parameter of interest alone are labelled as "Uncorr.-stat" in the solid red box. Boxes directly below the solid red box represent parameters common to the same p Z T bin as the parameter of interest, and are therefore non-migration parameters. The other boxes represent parameters in different p Z T bins, and are categorised as migration parameters. The categorisation can be broken down as follows: • Parameters in the dashed green box are from different coefficient numbers but the same p Z T bin and are labelled as "Shape" parameters.
• Parameters in the dotted blue box are from the same coefficient number (A 0 ) but in a different p Z T bin and are labelled as "Self-migration" parameters.
• The complements to these two categories are the parameters in the single-lined orange box and are labelled as "Shape-migration"; they are outside of both the chosen p Z T bin and coefficient number. These separations are done as well for the cross-section parameters, and are labelled as "Norm" and "Norm-migration" in the dot-dashed blue and double-lined purple boxes, respectively.
An illustration of this categorisation of the various components of the statistical uncertainty is shown for the ee CC + µµ CC , y Z -integrated measurement of the coefficient A 0 in Fig. 30.   Figure 30: Statistical uncertainty decomposition for the unregularised (left) and regularised (right) measurement of the A 0 coefficient in the ee CC + µµ CC channel for the integrated y Z configuration.

E. Quantifying A 5,6,7
The coefficients A 5,6,7 are expected to be zero at NLO, but are expected to receive NNLO contributions as large as 0.005 at high p Z T . The data measurements appear to be consistent with this, although the level at which the data measurements are non-zero should be quantified. A simple method to quantify this would be a standard χ 2 test of the measured spectra with respect to the null hypothesis of zero, but this has several disadvantages. First, the coefficients are expected to be non-zero only at high p Z T , and therefore a χ 2 test across the entire spectrum would be diluted by the low p Z T bins. Performing the test only for high p Z T could improve this locally, although this introduces some model dependence due to the choice of p Z T cutoff, as well as introducing a look-elsewhere-effect. Second, a χ 2 test is insensitive to the sign of the measured coefficients in each bin. Finally, it does not optimally account for positive or negative trends in the observation.
A signed covariant test statistic Q cov signed based on pseudo-experiments was developed for the purpose of quantifying the observed spectra. This takes into account pair-wise correlations between coefficients in neighbouring p Z T bins, as well as correlations between the different coefficients in the same p Z T bin. The contribution of each coefficient measurement to the test statistic is signed. Measurements below zero have a negative contribution, while measurements above zero have a positive contribution. Q cov signed is computed both on observed data and simulated data based on DYNNLO predictions at NNLO. The distribution of Q cov signed is obtained from ensemble tests under the null hypothesis. A p-value is obtained by integrating this distribution from the observed and simulated values to positive infinity, and converted to a one-sided statistical significance.
To compute Q cov signed (for any observed, simulated, or pseudo data), an initial set of pseudo-experiments are used to obtain the distribution ofÂ pseudo 5,6,7 . A fit is first performed to the data under the null hypothesis A 5,6,7 = 0 to obtainÂ null 0−4 ,σ null , andθ null . Pseudo-data is then generated in each likelihood bin around the expected events N n exp (Â null 0−4 ,σ null ,θ null ) (see Eq. (7)). A fit is performed to the pseudo-data to obtain A pseudo 5,6,7 .
Q cov signed is computed based onÂ 5,6,7 in any particular dataset in conjunction with the distribution ofÂ pseudo 5,6,7 from pseudo-data. It includes several components, which are described here. The significance Z i of the deviation from zero of each of the three A 5,6,7 coefficients in every p Z T bin is computed as depicted in Fig. 31. A weight, w i j = (sign(Z i )Z 2 i + sign(Z j )Z 2 j )/(Z 2 i + Z 2 j ), is computed based on the individual Z i values for every coefficient pair, both in coefficient number and bin in p Z T . A weight w i j has the property that w i j = +1 or −1 if Z i and Z j are both above or below zero, respectively, while it is a weighted difference between them otherwise. The correlation coefficient ρ i j between these pairs is extracted from their two-dimensional distributions. The pair-wise significance Z i j is computed from the tail probability of the measurement being more outward in it's quadrant than is observed. An example of this for each quadrant is shown in Fig. 31. Q cov signed is defined using these components as follows: The distribution of Q cov signed is finally obtained from a second set of pseudo-experiments. The observed value is computed along with the value from the DYNNLO expectation. The distribution is shown in Fig. 32, with vertical bars representing the observed and expected values. A total of 7800 pseudoexperiments were used in this computation. Integrating from the observed value to the right, the fraction of events in the tail is 0.14%, corresponding to a significance of 3.0σ. Similarly, the expected significance is 3.2σ.