First measurement of the $Z\rightarrow \mu^+ \mu^-$ angular coefficients in the forward region of $pp$ collisions at $\sqrt{s}=13$ TeV

The first study of the angular distribution of $\mu^+ \mu^-$ pairs produced in the forward rapidity region via the Drell-Yan reaction $pp \rightarrow \gamma^{*}/Z +X \rightarrow l^+ l^- + X$ is presented, using data collected with the LHCb detector at a centre-of-mass energy of 13TeV, corresponding to an integrated luminosity of 5.1 $\rm{fb}^{-1}$. The coefficients of the five leading terms in the angular distribution are determined as a function of the dimuon transverse momentum and rapidity. The results are compared to various theoretical predictions of the $Z$-boson production mechanism and can also be used to probe transverse-momentum-dependent parton distributions within the proton.

Drell-Yan muon pairs (µ + µ − ) produced near the Z-boson mass pole provide an excellent source of information about electroweak (EW) parameters, probe the proton structure and test the underlying strong-interaction dynamics of Z-boson production in proton-proton (pp) collisions, as described by quantum chromodynamics (QCD). With a dedicated detector instrumented in the forward region, the LHCb experiment plays a unique role in the study of Z → µ + µ − processes in the forward rapidity (y Z ) region [1], especially for y Z > 3.5. Therefore, the LHCb experiment is complementary to the ATLAS and CMS experiments, which are more efficient in the central region. The distributions of the Z-boson transverse momentum, p T , and the scattering angle, φ * η , of the muons with respect to the proton beam direction in the rest frame of the dimuon system [2] have already been measured by the LHCb collaboration [1,[3][4][5]. Key information is also encoded in the Z → µ + µ − angular distribution in the forward region, which has not been fully explored.
In the Drell-Yan neutral-current process, the lepton pair is produced via a spin-1 gauge boson. The kinematic distribution of the final-state leptons ( ± ) provides a direct probe of the polarization of the intermediate gauge boson, which is in turn sensitive to the underlying QCD production mechanisms. The QCD mechanisms of the process pp → γ * /Z + X → + − + X, where X represents any other particles, can be expressed using eight frame-dependent angular coefficients, A i (i = 0, ..., 7). The coefficients depend on the invariant mass, p T and rapidity of the lepton pair. At Born level, the angular distribution of the leptons in the boson rest frame is given by [6,7] dσ d cos θdφ ∝ (1 + cos 2 θ) + 1 2 A 0 (1 − 3 cos 2 θ) + A 1 sin 2θ cos φ + 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 φ, where θ and φ are the polar and azimuthal angles of the µ + lepton in the Collins-Soper frame [8]. At leading-order (LO) approximation within the framework of QCD, all angular coefficients vanish as the dilepton transverse momentum approaches zero, with the exception of A 4 , which is nonzero due to parity violation in the weak interaction. At next-to-leading-order (NLO), A 0−3 become nonzero, while A 5−7 have small deviations from zero only at next-to-next-to-leading-order (NNLO). The equality A 0 = A 2 is known as the Lam-Tung relation [9], which is valid for both qq and gq processes at LO [10]. However, small violations of the Lam-Tung relation occur in higher-order perturbative QCD calculations [11,12], as well as in QCD resummation calculations to all orders [13]. Therefore, precision measurements of A 0 − A 2 can test the Lam-Tung relation. Nonperturbative effects can lead to violation of the Lam-Tung relation via a correlation between the transverse spin and transverse momentum of the initial-state quark or antiquark, represented by the Boer-Mulders transverse-momentumdependent (TMD) parton distribution function (PDF) [14,15]. Significant violation of the Lam-Tung relation, with large cos 2φ modulations, has been observed in pion-induced Drell-Yan measurements [16][17][18]. The cos 2φ coefficient (A 2 ) in Eq. (1) is proportional to the convolution of two Boer-Mulders functions, of the quark and the antiquark [15]. Therefore, the angular measurements of the Drell-Yan process can improve constraints on nonperturbative partonic spin-momentum correlations within unpolarised protons via phenomenological extractions of the Boer-Mulders TMD PDF from these data in conjunction with data from other experiments. The Z-boson cross-section measurements at low p T ( 20 GeV/c) by the LHCb  This Letter reports the first measurement of the Z → µ + µ − angular coefficients in the forward rapidity region (2 < y Z < 5) by the LHCb experiment, using Z-boson candidates selected in the mass region 50 < M µµ < 120 GeV/c 2 , where M µµ is the invariant mass of the µ + µ − pair. The parameters A 0 , A 1 , A 2 and A 3 are determined as a function of the Z-boson p T and y Z , with the candidates selected in the mass region 75 < M µµ < 105 GeV/c 2 . Events in the low mass range, 50 < M µµ < 75 GeV/c 2 , and high mass range, 105 < M µµ < 120 GeV/c 2 , are used to probe the Boer-Mulders function; these regions are dominated by contributions from virtual photons and their interference with the Z-boson amplitude, with measurements in multiple mass regions adding sensitivity to the evolution of the TMD PDF with the hard scale. With the higher-order QCD effects related to the Z-boson p T being integrated out, the measurements of the angular coefficients as a function of y Z can be used to test the resummation calculations in the forward region. Due to the relatively small data sample containing high-p T (> 50 GeV/c) Z → µ + µ − events, the coefficients A 5 to A 7 are fixed to zero. Since this study focuses on probing the QCD dynamics of the Z-boson production, the A 4 coefficient, sensitive to the weak mixing angle [26], is not reported. Nevertheless, in order to investigate its variation across the kinematic range, the difference with respect to its mean value, ∆A 4 , is measured.
The LHCb detector [27, 28] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. A silicon-strip vertex detector [29] surrounding the pp interaction region allows c and b hadrons to be identified from their relatively long flight distance. A tracking system [30] provides a measurement of momentum, p, of charged particles, and two imaging Cherenkov detectors [31] are able to discriminate between different species of charged hadrons. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers [32]. The online event selection is performed by a trigger system [33], which consists of a hardware stage, based on information from the calorimeter and muon systems [34,35], followed by a software stage, which applies a full event reconstruction incorporating near-real-time alignment and calibration of the detector [36].
Simulated pp collisions are generated using Pythia8 [37] with a specific LHCb configuration [38]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [39] as described in Ref. [40]. Various theoretical predictions at the Born-level are compared with the measured results, including analytic resummed calculations, ResBos [41] and DYTurbo [42], and a fixed-order calculation with matching algorithms to veto the double counting of quantum electrodynamic (QED) final-state radiation (FSR) and parton shower, Powheg-BOX [43][44][45][46]. Predictions from the ResBos [41] generator are produced using the PDF known as CT14HERA2 [47]. The ResBos calculation combines NLO fixed-order perturbative calculations at high boson p T , with the soft-gluon Collins-Soper-Sterman resummation formalism [48][49][50] at low boson p T , with consideration of a nonperturbative contribution from the parton intrinsic transverse momentum, which is an all-orders summation of large terms from gluon emission. The DYTurbo [42] program is an NNLO generator for the calculation of the QCD transverse-momentum resummation of Drell-Yan cross sections up to next-tonext-to-leading logarithmic accuracy combined with the fixed-order results. DYTurbo includes the full kinematic dependence of the lepton pair with the corresponding spin correlations and finite-width effects. Powheg-BOX is an NLO generator and can be interfaced with Pythia8 for QCD and EW showering, but without the angular damping factors.
Measurements are performed using 5.1 fb −1 of data at a center-of-mass energy of 13 TeV collected by the LHCb detector during the years 2016-2018. For each Z → µ + µ − candidate, at least one of the muons is required to pass hardware and software singlemuon trigger decision stages. All muon candidates are required to have transverse momentum p µ T > 20 GeV/c and be in the range 2.0 < η < 4.5. The relative uncertainty in the momentum of the muon track must be less than 10%. To further suppress background sources, muon candidates are required to pass an isolation requirement which is parameterised as is the sum of the muon p T and that of the tracks within δR < 0.5, and δR is δη 2 + δΦ 2 . The quantities δη and δΦ give the separation between the muon and neighbouring tracks in pseudorapidity η and the laboratory azimuthal angle Φ, respectively. Candidate Z-bosons are composed of two muon candidates with opposite charge that originate from a common primary pp interaction vertex. In total, 818 074 (745 343) Z → µ + µ − candidates are selected in the mass region 50 < M µµ < 120 GeV/c 2 (75 < M µµ < 105 GeV/c 2 ). Background contributions from heavy-flavour processes and misidentified hadrons are estimated using data control samples [1], while background sources from W + W − , W ± Z, ZZ, Z → τ + τ − and tt processes are estimated from simulation. The total background contribution is determined to be 0.2%.
Detector misalignment effects are studied and suppressed using the Z → µ + µ − data events, where the mass peak position of the selected Z-boson candidates is calibrated in different kinematic and geometric regions to the world average value [51]. To improve the agreement between data and simulation, selection efficiencies determined from simulation are corrected to corresponding values measured in the data. The selection efficiencies include the muon trigger, muon track reconstruction, muon isolation and identification. These efficiencies are estimated using Z → µ + µ − data and simulation events, with the tag-and-probe method [1], as functions of muon η and Φ. Additionally, corrections to the momentum scale and resolution are applied to the simulation.
The angular coefficients are determined in six intervals of the Z-boson p T and five intervals of y Z , by fitting the two decay angles cos θ and φ of the selected Z → µ + µ − candidates, using an unbinned maximum-likelihood method in each interval. The probability density function used in the fit consists of a signal-only function, obtained from Eq. (1), convoluted with the detector resolution and acceptance effects. The angular coefficients (A 0−4 ) are the only free parameters in the fit, while the background contributions are effectively subtracted from the data sample by including background events with negative weights derived from simulation and data control samples. Simulated Pythia8 events are used to study the background modeling along with the detector resolution and acceptance effects, with FSR turned on at the generator level, thereby including the Born/Bare corrections [52]. This way, the measured A i coefficients are corrected to the Born level. The normalization weights method is used to avoid computing-intensive numerical integration of the probability density function in the fit. This method was introduced and used by the BaBar collaboration [53,54], and has been extensively used in LHCb heavy-flavour studies [55]. To avoid potential bias from the input angular coefficients in the simulation, an iterative method is used. The measured results are found to be stable after four iterations.
Several sources of systematic uncertainty are considered. Simulation is used to determine the detector acceptance effects in the fit process. To estimate the uncertainty from the finite number of simulated events, the bootstrap method [56] is used. This uncertainty varies from 0.0016 to 0.0402 in different kinematic regions, and is the dominant systematic uncertainty. Another significant uncertainty arises from the selection efficiencies. In the determination of the uncertainty from the efficiency corrections, two independent sources are considered: an uncertainty from the size of the control sample and one from the kinematic dependence of the corrections. For the second uncertainty, the default two-dimensional (muon η and φ) efficiency corrections are replaced with one-dimensional corrections, and variations of the results are taken as uncertainties. The fit method is validated with simulation samples. The angular coefficients of the simulation are changed to values far from the default (increased by a factor of 1.5, or set to 0). In the validation, the new angular coefficients are compared with the values from the fit, and differences between them are taken as uncertainties. The uncertainty from the muon momentum scaling and smearing are calculated by varying the corrections by their uncertainties, and the changes in the measured results are taken as the associated systematic uncertainties. These are found to be negligible.
The systematic uncertainties associated to the background are estimated by varying the background contribution by ±σ, where σ includes the uncertainties from the theoretical predictions of the cross-sections and the statistical uncertainties from the simulation. The systematic uncertainty from the FSR is estimated by comparing the measured results using the simulation sample weighted to the Powheg prediction with different FSR algorithms of Photos and Pythia8. Differences from different FSR approaches are taken as uncertainties. The systematic uncertainty from the PDF is estimated using 58 sub-PDF sets of CT18NNLO [57], with the prescription recommended by the PDF sets groups [58]. All systematic uncertainties are considered as uncorrelated, and the total systematic uncertainty is taken as their sum in quadrature.
The measured angular coefficients A 0−3 as well as ∆A 4 and the difference A 0 − A 2 , are presented as a function of the Z-boson p T in Fig. 1. The measurements are compared with the predictions of ResBos, DYTurbo, Powheg-BOX+Pythia8, and Pythia8 with LHCb tuning. The measurement of A 0 is in good agreement with DYTurbo, ResBos, and Powheg-BOX predictions, but disfavors the Pythia8 results in the high p T region (> 20 GeV/c). However, the Pythia8 calculations show similar increasing trends as the measurements in the high p T region. The measurement of A 1 disagrees with the Pythia8 result, but is in good agreement with other calculations in all p T regions. Within large statistical uncertainties, the measurements of A 2 are mostly in reasonable agreement with predictions, while for A 3 there is a noticeable underestimation of most predictions up to 55 GeV/c 2 . The measured ∆A 4 is in good agreement with theoretical predictions in all p T regions. Finally, for the A 0 − A 2 measurement, the LHCb results are compatible with zero only in the highest p T interval, but there the uncertainty is relatively large. Also, most predictions are not compatible with zero in the high-p T intervals. This violation of the Lam-Tung relations measured by the LHCb collaboration is consistent with previous results from the CMS [24] and ATLAS collaborations [25]. The measured A 0 − A 2 results are significantly smaller than the predictions at low p T (< 20 GeV/c) and slightly smaller at high p T (> 80 GeV/c), while larger than the predictions in other p T regions. A summary of the full results, including measurements in bins of y Z , is provided in the Supplemental Material [59]. With the exception of A 0 , the measured angular coefficients do not vary significantly as a function of y Z , when the measurement is performed for the full p T range. There is reasonable agreement between the measurements and ResBos calculations for A 1−3 and ∆A 4 , while for A 0 there is tension between the measurements and predictions, especially in the highest y Z region. The p-value between measurements and theoretical predictions is calculated to be 0.06. The differences between A 0 and the ResBos predictions in the y Z distribution indicate that there is a y Z -dependence in the QCD resummation or additional higher-order effects. More detailed studies on the predictions as a function of y Z are necessary.
The nonperturbative Boer-Mulders TMD PDF can be probed with the measurement of A 2 in the lower p T region. Measurements of A 2 in the low, middle and high M µµ regions are shown in Fig. 2, where the dimuon p T region is divided into four intervals from 0-20 GeV/c. Despite the limited sample size with the finer p T intervals, several observations can be made. In the low M µµ , the measured A 2 value in the lowest dimuon p T region (p T < 3 GeV/c) deviates significantly from all predictions. It is unclear if nonperturbative spin-momentum correlations in the proton, described by the Boer-Mulders distribution, could lead to such variations as no phenomenological calculations are available. None of the calculations used in the comparison include this type of nonperturbative effect. In other p T regions, reasonable agreement between measurements and predictions is seen.
In summary, the first measurements of the angular coefficients of Drell-Yan µ + µ − pairs in the forward rapidity region of pp collisions are presented. These quantities provide more information on the Z-boson production compared to other traditional observables. The measurements, which are given as a function of the Z-boson p T and y Z , in both the Z peak and lower and higher M µµ regions, are compared with various predictions. Some tension between the measurements and the theoretical predictions appears in the lower M µµ and low-p T region, and in the y Z dependence (See the Supplemental Material [59] for further details). More studies of the theoretical models are needed to shed light on the apparent discrepancies. This analysis provides important and unique inputs for future phenomenological extractions of spin-momentum correlations in the proton in terms of the Boer-Mulders TMD PDF and its evolution with the hard scale.

Supplemental material
The number of data candidate events in different Z-boson p T and rapidity intervals are summarized in Table 1. Summaries of measured angular coefficients A i are shown in Table 2 in intervals of Z-boson p T , and in Table 3 and Fig. 3 in intervals of Z-boson rapidity. Summaries of measured angular coefficients and regularised uncertainties of A 2 in the Z-boson low-p T region, with events in different mass regions, are given in Table 4.  Table 2: Summary of measured angular coefficients and regularised uncertainties for A i and A 0 − A 2 , in intervals of Z-boson p T , in the region of Z-boson y > 2 and 75 < M µµ < 105 GeV/c 2 . The total systematic uncertainty is shown with the breakdown into its underlying components. Entries marked with '-' indicate that the uncertainty is below 0.0001.  p T (Z) ∈ [35,55][ GeV]   Table 3: Summary of measured angular coefficients and regularised uncertainties for A i and A 0 − A 2 , in intervals of Z-boson rapidity, in the region of Z-boson p T < 100 GeV/c and 75 < M µµ < 105 GeV/c 2 . 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 4: Summary of measured angular coefficients and regularised uncertainties for A 2 , in intervals of Z-boson p T , in the region of Z-boson y > 2 and different mass ranges. The total systematic uncertainty is shown with the breakdown into its underlying components. Entries marked with '-' indicate that the uncertainty is below 0.0001. LHCb collaboration