Z boson production in Pb+Pb collisions at √ s NN = 5 . 02 TeV measured by the ATLAS experiment

The production yield of Z bosons is measured in the electron and muon decay channels in Pb+Pb collisions at √ s NN = 5 . 02 TeV with the ATLAS detector. Data from the 2015 LHC run corresponding to an integrated luminosity of 0.49 nb − 1 are used for the analysis. The Z boson yield, normalised by the total number of minimum-bias events and the mean nuclear thickness function, is measured as a function of dilepton rapidity and event centrality. The measurements in Pb+Pb collisions are compared with similar measurements made in proton–proton collisions at the same centre-of-mass energy. The nuclear modiﬁcation factor is found to be consistent with unity for all centrality intervals. The results are compared with theoretical predictions obtained at next-to-leading order using nucleon and nuclear parton distribution functions. The normalised Z boson yields in Pb+Pb collisions lie 1–3 σ above the predictions. The nuclear modiﬁcation factor measured as a function of rapidity agrees with unity and is consistent with a next-to-leading-order QCD calculation including the isospin effect. 3 .


Introduction
The measurement of electroweak (EW) boson production is a key part of the heavy-ion (HI) physics programme at the Large Hadron Collider (LHC). Isolated photons and heavy vector bosons, Z and W , are powerful tools to probe the initial stages of HI collisions. After being created at the initial stage of the collision in high-momentum exchange processes, Z and W bosons decay much faster than the timescale of the medium's evolution. Their leptonic decay products are generally understood to not be affected by the strong interaction; hence they carry information about the initial stage of the collision and the partonic structure of the nuclei.
Measurements performed by the ATLAS and CMS experiments with Z and W bosons decaying leptonically show that production rates of these non-strongly interacting particles are proportional to the amount of nuclear overlap, quantified by the mean nuclear thickness function, T AA [1-6]. Results obtained with isolated high-energy photons [7,8] are also consistent with this observation.
The transverse momentum and rapidity distribution of Z bosons and the pseudorapidity distribution of muons originating from W bosons measured in Pb+Pb collisions at √ s NN = 2.76 TeV have been found to be generally consistent with perturbative quantum chromodynamics (pQCD) calculations of nucleon-nucleon E-mail address: atlas .publications @cern .ch.
collisions scaled by T AA . Production of Z bosons in Pb+Pb collisions was found to be consistent with next-to-leading-order (NLO) pQCD calculations that do not include nuclear modifications in the treatment of parton distribution functions (PDFs). However, some nuclear modification of PDFs could not be excluded within the precision of the existing Pb+Pb measurements [1,2,7]. The recent ALICE result at √ s NN = 5.02 TeV shows better agreement with nPDF calculations at forward rapidities [9]. On the other hand, the study of asymmetric p+Pb collisions at √ s NN = 5.02 TeV shows that including nuclear modifications of PDFs gives a substantially better description of the data than using a free proton PDF. This is seen by comparing the Z boson cross section in p+Pb collisions with pQCD calculations [10-13] and recently also the W boson cross section at √ s NN = 5.02 TeV and √ s NN = 8.16 TeV [13,14]. In addition, studies of Z bosons differentially in p+Pb centrality demonstrate that they are a sensitive test of the Glauber model description of nuclear geometry [11]. This letter presents results on Z boson production yield measurement in the Z → μμ and Z → ee decay channels in Pb+Pb collisions at √ s NN = 5.02 TeV with the ATLAS detector at the LHC.
The data sample was collected in November 2015 and corresponds to an integrated luminosity of 0.49 nb −1 . The observables under study are the yield of produced Z bosons in the fiducial kinematic region defined by detector acceptance and lepton kinematics normalised to the number of minimum-bias events, measured differentially in rapidity and event centrality. The Pb+Pb data are compared with pQCD calculations, and the nuclear modification factor is measured relative to pp cross section previously measured by the ATLAS experiment [15].

ATLAS detector
The ATLAS detector [16] covers nearly the entire solid angle 1 around the collision point. It consists of an inner tracking detector surrounded by a thin superconducting solenoid, electromagnetic and hadronic calorimeters, and a muon spectrometer incorporating three large superconducting toroid magnets.
The inner-detector system is immersed in a 2 T axial magnetic field and provides charged-particle tracking in the range |η| < 2.5.
The high-granularity silicon pixel detector covers the vertex region and typically provides four measurements per track, the first hit being in the insertable B-layer [17,18] in operation since 2015. It is followed by the silicon microstrip tracker, which usually provides eight measurements per track. These silicon detectors are complemented by the transition-radiation tracker, which enables radially extended track reconstruction up to |η| = 2.0.
The calorimeter system covers the pseudorapidity range |η| < 4.9. Within the region |η| < 3.2, electromagnetic (EM) calorimetry is provided by barrel and endcap high-granularity liquid-argon (LAr) sampling calorimeters, with an additional thin LAr presampler covering |η| < 1.8, to correct for energy loss in material upstream of the calorimeters. Hadronic calorimetry is provided by the scintillator-tile calorimeter, segmented into three barrel structures within |η| < 1.7, and two LAr hadronic endcap calorimeters. The forward calorimeter (FCal) is a LAr sampling calorimeter located on either side of the interaction point. It covers 3.1 < |η| < 4.9 and each half is composed of one EM and two hadronic sections. The FCal is used to characterise the centrality of Pb+Pb collisions as described below. Finally, zero-degree calorimeters (ZDC) are situated at large pseudorapidity, |η| > 8.3, and are primarily sensitive to spectator neutrons.
The muon spectrometer comprises separate trigger and highprecision tracking chambers measuring the deflection of muons in a magnetic field generated by superconducting air-core toroids. The precision chamber system covers the region |η| < 2.7 with three layers of monitored drift tubes, complemented by cathode-strip chambers in the forward region, where the background is highest. The muon trigger system covers the range |η| < 2.4 with resistiveplate chambers in the barrel, and thin-gap chambers in the endcap regions.
A two-level trigger system is used to select events of interest for recording [19]. The level-1 (L1) trigger is implemented in hardware and uses a subset of the detector information to reduce the event rate. The subsequent, software-based high-level trigger (HLT) selects events for recording. Both the electron and muon event selection used in this analysis combine L1 and HLT decision algorithms.

Data sets and event selection
All of the analysed data were recorded in periods with stable beam, detector, and trigger operations. Candidate events are required to have at least one primary vertex reconstructed from the inner-detector tracks. In addition, a trigger selection is applied, requiring a muon or an electron candidate above a p T threshold of 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). 8 GeV or 15 GeV, respectively. The electron-trigger candidate is further required to satisfy a set of loose criteria for the electromagnetic shower shapes [20]. The trigger algorithm implements an event-by-event estimation and subtraction of the underlying-event contribution to the transverse energy deposited in each calorimeter cell [21]. For both the electron and muon candidates, further requirements are applied to suppress electromagnetic background contributions, as described in Section 4.2.
Muon candidates reconstructed offline must satisfy p T > 20 GeV and |η| < 2.5 and pass the requirements of 'medium' identification optimised for 2015 analysis conditions [22]. Offline selected electron candidates are required to have p T > 20 GeV and |η| < 2.47, although candidates within the transition region between barrel and endcap calorimeters (1.37 < |η| < 1.52) are rejected. In addition, likelihood-based identification is applied, developed for the Pb+Pb data conditions and based on a general strategy described in Ref. [23].
Events with a Z boson candidate are selected by requiring exactly two opposite-charge muons or electrons, at least one of which is matched to a lepton selected at trigger level. The dilepton invariant mass must satisfy the requirement 66 < m < 116 GeV consistent with previous ATLAS measurements. A total of 5347 Z boson candidates are found in the muon channel and 4047 in the electron channel.
In order to estimate the geometric characteristics of HI collisions, it is common to classify the events according to the amount of nuclear overlap in the collision. The quantity used to estimate the collision geometry is called the 'collision centrality'. The centrality determination is based on the total transverse energy measured by both FCal detectors in each event, E FCal T . This quantity is then mapped to geometric quantities, such as the average number of participating nucleons, N part , and the mean nuclear thickness function, T AA , which quantifies the amount of nuclear overlap in a centrality class and is evaluated using a Glauber calculation [24,25]. The mapping is based on specific studies of an event sample without additional Pb+Pb collisions within the same or neighbouring bunch crossings (pile-up) collected with minimum-bias (MB) triggers. A special treatment is employed for events in the 20% most peripheral interval, where diffractive and photonuclear processes contribute significantly to the MB event sample. This requires extrapolating from the total number of MB events in this region and employing a special requirement on the Z boson event topology, as described in Section 4.2. Table 1 summarises the relationship between centrality, N part , and T AA as calculated with Glauber MC v2.4 [6,26], which incorporates nuclear densities averaged over protons and neutrons. The total number of MB events in the 0-80% centrality interval is (2.99 ± 0.04) × 10 9 , which is then distributed in different centrality intervals according to their size. The quoted uncertainty on the number of MB events includes variations on the E FCal T value corresponding to the 0-80% centrality interval estimated with the Glauber model. This sample is obtained by selecting events passing MB triggers and excluding the events with a pile-up contribution, where the total sampled integrated luminosity corresponds to the signal selection [25].
Simulated samples of Monte Carlo (MC) events are used to evaluate the selection efficiency for signal events and the contribution of several background processes to the analysed data set. All of the samples were produced with the Geant4-based simulation [27,28] of the ATLAS detector. Dedicated efficiency and calibration studies with data are used to derive correction factors to account for residual differences between experiment and simulation.
The processes of interest containing Z bosons were generated with the Powheg-Box v1 MC program [29][30][31][32] interfaced to the Pythia 8.186 parton shower model [33]. The CT10 PDF set [34] was used in the matrix element, while the CTEQ6L1 PDF set [35]  A sample of top-quark pair (tt) production was generated with the Powheg-Box v2 generator, which uses NLO matrix-element calculations [38] together with the CT10f4 PDF set [39]. The parton shower, fragmentation and underlying event in nucleonnucleon collisions were simulated using Pythia 6.428 [40] with the CTEQ6L1 PDF set and the corresponding Perugia 2012 tune (P2012) [41]. The top-quark mass was set to 172.5 GeV. The Evt-Gen v1.2.0 program [42] was used to model bottom and charm hadron decays for all versions of Pythia. The total Z boson and top-quark yields in MC samples are normalised using the results of NLO QCD calculations.
The signal MC samples were produced with different nucleonnucleon combinations (pp, pn, nn) weighted to reflect the isospin composition of lead nuclei. For lead, A = 208 and Z = 82, so all samples with two neutrons have a weight of ([A − Z ]/Z ) 2 =36.7%, events with two protons have a weight of 15.5% and the unlike nucleon combinations (pn, np) each have weights of 23.9%.
Once produced, the simulated events were overlaid with MB events taken during the Pb+Pb run. The overlay of data events was done such that the MC simulation accurately reflects detector occupancy conditions present in the Pb+Pb run. The MB events used for the overlay were sampled such that the centrality distribution, based on the total transverse energy deposited in the forward calorimeters, approximates that of Z boson events, which are generally biased to more-central collisions. The simulated events were finally reconstructed by the standard ATLAS reconstruction software.

Measurement procedure
The differential Z boson production yield per MB event is measured within a fiducial phase space defined by p T > 20 GeV, |η | < 2.5 and 66 < m < 116 GeV. The yields in both the electron and muon channel are calculated using where N Z and B Z are the number of selected events in data and the expected number of background events, respectively, and Z trig is the trigger efficiency per Z boson candidate measured in data and described in Section 4.3.
A correction for the reconstruction efficiency, momentum resolution and the final-state radiation effects is applied with the bin-by-bin correction factor C Z which is obtained from MC simulation as is the number of events passing the signal selection at the detector level. The number of selected events is corrected for measured differences between data and simulation in lepton reconstruction and identification efficiencies. The denominator N MC,fid Z is computed by applying the fiducial phase-space requirements to the generator-level leptons originating from Z boson decays. The measurement is corrected for QED final-state radiation effects by using the generator-level lepton kinematics before photon radiation.
The value of C Z in the fiducial acceptance averaged over all centralities is determined from MC simulation after reweighting as explained in Section 3. It is 0.659 and 0.507 in the Z → μ + μ − and Z → e + e − decay channels, respectively. The uncertainty due to the size of the simulated sample is at the level of 0.1% for each decay channel and is not the dominant uncertainty.
The rapidity, momentum and centrality dependence of C Z is calculated from the simulation as where F (p T , y) is the centrality-averaged efficiency calculated per y and p T interval of the dilepton system and G( y, E FCal T ) is a parabolic parameterisation of a correction factor accounting for the centrality and rapidity dependences of the efficiency. In each rapidity bin, the factor G is obtained from a fit of the ratio of the efficiency in a particular centrality bin to the value averaged over all possible centrality values. Nuclear modification is quantified by measuring the ratio of the Z boson production rate, scaled by the mean nuclear thickness function, to the Z boson production cross section in pp collisions, a quantity known as the nuclear modification factor: where T AA is the nuclear thickness function in a given central- Pb+Pb /dy is the differential yield of Z bosons per inelastic MB event measured in Pb+Pb collisions and dσ Z pp /dy is the differential Z boson cross section measured in pp collisions [15]. A deviation from unity in R AA indicates the nuclear modification of the observable. The value of R AA is expected to be greater than unity by about 2.5%, based on MC simulation, due to the higher Z boson production cross section in proton-neutron and neutron-neutron interactions which are present in Pb+Pb collisions and amount to 84.5% of the total hadronic cross section. This is later referred to as the 'isospin effect' and is not accounted for in the definition of R AA .

Background determination
There are two background source categories studied in this analysis. The first includes the same background sources that are studied in pp collisions [15] and the second includes additional background sources specific to the Pb+Pb collision system.
Background contributions in the first category are expected from Z → τ + τ − , top-quark pair production and multi-jet events.
The first two contributions are evaluated from dedicated simulation samples, whereas the multi-jet background contribution is derived using a data-driven approach. The Z → τ + τ − background is found to be 0.05% of all signal candidates in the muon channel and 0.06% in the electron channel. The top-quark background amounts to 0.08% in the muon channel and 0.05% in the electron channel. The background contribution from W boson decays and W +jet production is found to be negligible. The multi-jet background originates from jets, misidentified hadrons and, in the electron channel, from converted photons. In the muon channel, its contribution is estimated from the distribution of the same-sign Z boson candidates in rapidity and centrality. Due to the low charge misidentification rate in the muon spectrometer, their invariant mass distribution does not exhibit a peak in the Z boson mass region. In the muon channel this background amounts to 0.5% of all signal candidates. In the electron channel, there is a significant contribution from charge misidentification, fakes and conversions. The electron same-sign pairs therefore cannot be used to estimate the multi-jet background. This contribution to the selected event sample in the electron channel is estimated using a background template obtained from the data in Z boson rapidity and event centrality. The template is derived from a subset of the signal sample that corresponds to electrons from jets, i.e. electrons with a very poor reconstruction quality [23] that is the total transverse momentum measured inside the cone of size R = 0.2 around the electron track. The template is normalised to the number of same-sign data candidates in the low-mass region of 60 < m ee < 70 GeV after the signal MC subtraction. Due to the small number of signal candidates satisfying this condition, samesign electron pairs with the same kinematic requirements are also added to this background template. The shape of the obtained multijet template is shown in Fig. 1. This background amounts to 2% of all signal candidates in the electron channel.
The background contributions specific to Pb+Pb come from two main sources. The first is due to pile-up when more than a single Pb+Pb collision is recorded simultaneously or in a nearby bunch crossing. The second is the production of additional Z boson can-didates by photon-induced reactions produced by the intense electromagnetic fields generated by the colliding ions (below referred to as 'electromagnetic background'). Pile-up distorts the transverse energy measured in the FCal and causes reconstructed Z bosons to be assigned to an incorrect centrality interval. Pile-up events from other collisions in the same bunch crossing (in-time pile-up) increase the E FCal T , shifting the Z boson candidate to a morecentral interval. Alternatively, if a pile-up collision precedes the trigger event (out-of-time pile-up), its contribution to the E FCal T can be negative, due to the time response of the electronic signal shapers used in the calorimeters [43]. In this case, the Z boson candidate is shifted to a more-peripheral interval. Both processes depend on the instantaneous luminosity during data taking. At any time during the HI run the number of hadronic interactions per bunch crossing was less than 0.01. To preserve the accuracy of the total yield measurement, no pile-up removal procedure was applied to the selected events. However, due to the fact that the Z boson production scales linearly with T AA [1], the increase in the FCal transverse energy in in-time pile-up events transfers candidates from less populated to more-populated centrality intervals, thus having a very small effect and changing the average number of counts in the most central collisions by an insignificant amount. Contrary to that, the reduction in the E FCal T transfers out-of-time pile-up events from more-populated to less-populated centrality intervals, thus making a larger relative contribution to the moreperipheral events. The effect has been studied using several independent data-driven and simulation-based approaches. The largest contribution to the most peripheral 80-100% centrality interval due to this type of the pile-up is less than 2%, i.e. less than one count, and is significantly less in any other centrality range.
A non-negligible relative contribution to the dileptons in the Z boson mass range in peripheral centrality intervals is expected from electromagnetic background sources. On the other hand, the expected rate of signal events in those peripheral centrality bins is low. Two photon-induced processes are expected to contribute to the background: photon-photon scattering, γ γ → + − [44][45][46] and photon-nucleus scattering γ + A → Z → + − [47]. Although measurements of exclusive high-mass dilepton production have been performed in pp and Pb+Pb collisions by ATLAS [48] and in pp collisions by CMS [49] a photo-nuclear Z boson production has not yet been observed in HI collisions. When the impact parameter of the photon-induced processes is larger than twice the nuclear radius such processes are referred to as ultra-peripheral collisions, and in these they are not obscured by hadronic interactions. Both physics processes are characterised by large rapidity gaps on one or both sides of the detector (regions with no particle production recorded in the detector), which are used in this analysis to measure and subtract these backgrounds. The rapidity gap estimation is implemented using a similar technique as developed in Ref. [50].
In the 50-100% (peripheral) centrality intervals, there is an additional requirement of a ZDC signal coincidence in order to suppress the electromagnetic background contributions. The energy measured in either detector is required to be at least 1 TeV, corresponding to 40% of the energy deposition of a single neutron. Without using any ZDC coincidence requirement in the event selection, 34 events with a rapidity gap greater than 2.5 units are found in the two decay channels. Since the estimated number of hadronic Z boson candidates with such a gap is below 0.05, all of these events are considered to be produced by photon-induced reactions and are removed from the sample. Events without gaps can have a photon-induced dilepton pair as well, if the rapidity gap is filled by particles originating from a simultaneous nucleonnucleon interaction. These events would appear in the centrality intervals defined by the E FCal T deposition from the hadronic interaction.
Following Ref.
[51], photon-induced reactions occurring simultaneously with hadronic collisions can be identified using both the angular and momentum correlations of final-state dilepton pairs. One variable used to quantify these correlations is the dilepton acoplanarity, defined as α ≡ 1 − |φ + − φ − |/π , where φ ± are the azimuthal angles of the two opposite-charge leptons. The same observable is used in this analysis to extract the contribution of photon-induced reactions to the measured Z boson production. Based on the MC signal simulation and measurements in the 0-50% centrality interval, (13.3 ± 0.4)% of Z boson candidates produced by hadronic collisions have acoplanarity below 0.01. On the other hand, among the 34 events with a rapidity gap greater than 2.5 units which were rejected from the sample as pure photoninduced events, 26 are found to have α < 0.01, corresponding to a fraction of 76%. This demonstrates that the acoplanarity is sensitive to photon-induced reactions in the Z boson invariant-mass region. This allows the photon-induced background to be estimated in all centrality intervals by comparing the number of Z boson candidates in a given centrality interval with the number of candidates with α < 0.01.
It is estimated that in the 80-100% centrality interval, besides the events with a large rapidity gap, 7 ± 3 out of 28 remaining candidates originate from photon-induced reactions. In the 60-80% interval this number is 15 ± 5 out of 182, and in the 50-60% interval it is 18 ± 8 out of 258, where the uncertainty includes statistical and systematic uncertainties in the background, but not in the number of event candidates. In more-central collisions, the method of estimating the photon-induced background is limited by the statistical precision of the sample and this contribution is consistent with zero. Rapidity distributions in centrality intervals are normalised to the estimated number of signal events and an additional systematic error is assigned in each interval to account for the EM background subtraction procedure. Fig. 1 shows the invariant mass distribution of backgrounds from all considered sources for both decay channels integrated in event centrality. The EM background shape shown in Fig. 1 is derived from the events containing large rapidity gaps and is normalised to the total estimated number of background candidates quoted above.

Detector performance corrections
After subtracting background contributions, the number of Z boson candidates is corrected for the trigger efficiency and detection efficiency, according to Eq. (1). All the correction factors are derived directly from the current data set used in the analysis, with the exception of the lepton momentum calibration corrections that are derived from pp collision data, and extrapolated to the Pb+Pb dataset conditions. The trigger efficiency per reconstructed Z candidate Z trig is derived from the efficiency of the single-lepton trigger via the relation Z where the indices refer to the two leptons forming the candidate pair. In order to obtain Z trig as a function of the dilepton p T and y which is further applied as a correction per dilepton candidate, kinematic distributions of the decay products are taken from MC simulation.
Muon and electron trigger efficiencies are measured using a tag-and-probe method [19,22,23] as a function of η and φ. The tag lepton is required to be reconstructed with high quality and very low probability of background contamination and to be matched to a lepton selected at trigger level. The probe lepton, satisfying the analysis reconstruction and identification requirements, is paired to it to give an invariant mass in the range 66 < m < 116 GeV. The background contribution to this measurement is estimated from the number of same-sign pairs and amounts to 0.8% and 3.5% in the muon and electron channels, respectively. The single-muon trigger efficiency in the endcap region of the detector (1.05 < |η| < 2.4) is measured to be around 85%, and in the barrel region (|η| < 1.05) it varies between 60% and 80%. A significant dependence of the efficiency on the muon azimuthal angle φ was measured and thus the trigger efficiency correction is derived as a function of both φ and η. The single-electron trigger efficiency is measured to be around 95% in the endcap region of the calorimeter (1.52 < |η| < 2.47) and it increases slightly to 97% in the barrel (|η| < 1.37). A significant dependence on the electron p T was measured, and the efficiency rises from 85% to 97% in the range from 20 to 100 GeV integrated over η. The single-electron trigger efficiency is thus derived as a function of p T and η. The Selection requirements including the muon reconstruction and identification are imposed on muon candidates used in the analysis. The efficiency of the selection criteria is measured using a tag-and-probe method in Z → μ + μ − events [22] and compared with simulation. Ratios of the efficiencies determined in data and simulation are applied as scale factors (SF) to correct the simulated events. Since the measured efficiencies are found to have negligible dependence on the muon momentum in the selected kinematic region and a very weak centrality dependence, the SF are evaluated only as a function of muon η. The centrality dependence of the SF is taken into account in the evaluation of systematic uncertainties. The combined reconstruction and identification efficiency for medium-quality muons typically exceeds 84% in both the data and simulation with good agreement between the two estimates.
The largest difference is observed in the endcap region (|η| > 1.8).
Electron candidates used for the analysis are required to satisfy selection criteria related to reconstruction and identification. The efficiency of the selection is measured using a tag-and-probe method in Z → e + e − events, as described in Ref. [23], and compared with simulation to derive electron scale factors. Measurements are performed as a function of the electron η, p T and event centrality. The electron reconstruction and identification efficiency is measured to be typically 70% in the endcap (|η| > 1.52) with good agreement between the data and the simulation. The SF is measured to be 1% away from unity with a precision of 3% in that region. In the barrel region (|η| < 1.37) the efficiency is measured to be around 80% while in the MC simulation the efficiency reaches 85%. Therefore, in this region a significant SF is applied, measured with a precision of 3-5%.
The lepton momentum scale and resolution corrections are derived using the pp signal MC samples and are applied to the simulation for both the electrons and muons. For the reconstructed muons, these corrections are derived as a function of the muon η and φ [22]. The correction factors are chosen such that they minimise the χ 2 between the muon-pair invariant mass distributions in data and simulation. The energy scale of reconstructed electrons is corrected by applying to the data a per-electron correction factor. The momentum scale correction factors are derived from a comparison of the electron-pair invariant mass between simulation and data.

Systematic uncertainties
In both channels, the trigger efficiency is derived from the tag-and-probe results in the data. The statistical limitation of the measured sample determines the uncertainty associated with the trigger correction. Although the uncertainties in each bin are relatively large in the muon trigger efficiency, after propagating this uncertainty to the dimuon efficiency, where only one of the muons is required to fire the trigger, the total uncertainty is quite small, between 0.2% and 0.5% when derived as a function of centrality and 0.1-0.2% when derived as a function of Z boson rapidity. The uncertainty is propagated using MC pseudo-experiments and the uncertainties in the linear fit coefficients of the trigger efficiency as a function of centrality. In the electron channel this uncertainty is at most 0.5%.
The NLO cross section of the background samples of Z → τ + τ − and top-quark production is varied by 10% to derive the corresponding uncertainties [15]. However, in both decay channels the multi-jet background dominates the uncertainty contribution. In the muon channel, the multijet background contribution is varied by 10% according to the level of statistical uncertainty in the number of same-sign candidates. This source produces 0.01-0.1% uncertainty in rapidity bins and up to 0.2% uncertainty in centrality bins. In the electron channel, the multi-jet template normalisation is varied by 20%, which corresponds to the level of statistical uncertainty in the number of same-sign candidates in the low-mass control region. The overall contribution to the systematic uncertainty is about 0.5% in rapidity bins and 0.5-2% in centrality bins.
In the 50-100% centrality interval, the uncertainty in subtracting the photon-induced background is evaluated by considering two sources. The first source is the compatibility between the acoplanarity cut efficiencies for hadronic Z boson production evaluated from data and simulation. An uncertainty of 0.4% accounts for this difference. For the second source, the uncertainty in the background rejection efficiency of the acoplanarity cut is evaluated from the candidates with large rapidity gaps. This uncertainty has two contributions. One is the statistical uncertainty of the event sample with large rapidity gaps, which amounts to 7%. The other contribution is due to the observed differences between the acoplanarity distributions for electrons and muons, which amounts to about 8%. In the 0-50% centrality interval, where the background subtraction is not performed, an uncertainty of 0.4%, evaluated from the difference between the data and simulation acoplanarity distributions, is introduced to account for a possible residual EM background contribution.
Uncertainties in the determination of lepton reconstruction and identification efficiency scale factors, as well as the parameterisation of the centrality dependence of the total correction affect the measurements through the correction factors C Z .
In the muon channel, the scale factors in the three η regions described in Section 4.3 are modified by their errors to derive the corresponding systematic uncertainty of C Z . In addition, the impact of the measured SF dependence of the final C Z value on the event centrality is also evaluated. The total relative uncertainty from these two sources ranges from 3.1% at midrapidity (| y| < 0.5) to 4.5% at forward Z boson rapidities and gives a contribution, constant as a function of the event centrality, of ∼3.4% for Z boson yields.
The main contribution to the systematic uncertainty in the electron channel comes from the uncertainties in measuring the reconstruction efficiency scale factor. Uncertainties related to this efficiency are classified as either correlated or uncorrelated, and are propagated accordingly to the final measurement uncertainty. The correlated uncertainty component of the SF is obtained by varying the requirements on the tag electron identification and isolation and on the invariant mass of the tag-and-probe pair. The statistical, uncorrelated, components of the scale factor uncertainties are propagated to the measurements via MC pseudo-experiments, while the systematic components are propagated as a single variation fully correlated across all electron η intervals. This source gives a 2.5-5% uncertainty as a function of rapidity and around 3% for all centrality intervals.
The effect of the calibration and energy scale correction uncertainty of electrons and muons is negligible.
An additional uncertainty of the bin-by-bin correction is due to the parameterisation of the rapidity and centrality dependence of C Z described in Eq. (2) and it stems primarily from the statistical uncertainty of the MC data set. To estimate uncertainties associated with these assumptions the parameters of the function G( y, E FCal T ) are varied by the errors of the parameters of the parabolic fit including covariance between the parameters. The data are corrected with these variations, and the difference between these results and the standard correction are taken as an estimate of the systematic uncertainty. In the muon channel the uncertainty associated with this source ranges from 0.4% to 1.4% in rapidity bins and is constant at ∼0.5% as a function of centrality, with the exception of the most peripheral bin where the uncertainty is ∼1%. In the electron channel, the uncertainty as a function of the rapidity ranges from 0.4% to 1.6% and is ∼1% for most centrality intervals, although in the most peripheral bin this contribution rises to ∼2%.
For both channels and the combined result, the uncertainties of the geometric parameters ( T AA and N part ) listed in Table 1 range from about 1% in central collisions to about 12% in peripheral collisions. These are treated as fully correlated between the channels. Finally, the total uncertainty for the pp measurement [15] used for the R AA calculation is 2.3%.

Results
The rapidity distributions of the Z boson yield for the muon and electron decay channels, normalised by the number of MB events and mean nuclear thickness function T AA , are shown in the upper-left panel of Fig. 2. The upper-right panel of the figure shows the N part dependence of the normalised Z boson yield in the fiducial acceptance where the systematic uncertainties of the N part values are scaled by a factor of three for clarity. The measurements performed in the two channels are combined using the Best Linear Unbiased Estimate (BLUE) method [52], accounting for the correlations of the systematic uncertainties across the channels and measurement bins. The combined result is shown in Fig. 2 together with the combined statistical and systematic uncertainties. The level of agreement between the channels shown in the lower  panels of the figure is quantified as χ 2 /N dof = 1.7/5 as a function of rapidity and χ 2 /N dof = 21.6/14 as a function of centrality.
The measured Z boson yields are compared with theoretical predictions obtained using a modified version of DYNNLO 1.5 [53,54] optimised for fast computations. The calculation is performed at O (α S ) in QCD and at leading order in the EW theory, with parameters set according to the G μ scheme [55]. The input parameters (the Fermi constant G F , the masses and widths of W and Z bosons, and the CKM matrix elements) are taken from Ref. [56]. The DYNNLO predictions are calculated using the free proton PDF set CT14 NLO [57] typically used to compare with the pp data and, additionally, the nuclear PDF sets nCTEQ15 NLO [58] and EPPS16 NLO [59], which are averaged over each Pb nucleus. In addition, the parton-level NLO prediction from the MCFM code [60], interfaced to the CT14 NLO PDF set, is calculated. This takes into account the isospin effect, due to different partonic compositions of protons and neutrons in the Pb nuclei, which is neglected in the DYNNLO calculations. The renormalisation and factorisation scales, respectively denoted by μ r and μ f , are set to the value of lepton pair invariant mass. The uncertainties of these predictions are derived as follows. The effects of PDF uncertainties are evaluated from the variations corresponding to each NLO PDF set. Uncertainties due to the scales are defined by the envelope of the variations obtained by changing μ r and μ f by a factor of two from their nominal values and imposing 0.5 ≤ μ r /μ f ≤ 2. The uncertainty induced by the strong coupling constant is estimated by varying α S by ±0.001 around the central value of α S (m Z ) = 0.118, following the prescription of Ref. [57]; the effect of these variations is estimated by comparing the CT14nlo_as_0117 and CT14nlo_as_0119 PDF sets [57] to CT14 NLO. Imperfect knowledge of the proton PDF and the scale variations are the main contributions to the total theory uncertainties. In calculating the R AA predictions, only the nuclear PDF uncertainties contribute since the CT14 NLO uncertainties cancel.
In Fig. 3 the normalised Z boson yield is compared between the combined measurement and the theoretical predictions calculated with the CT14, nCTEQ15 and EPPS16 NLO PDF sets, with uncertainties assigned as previously described. All calculations lie 1-3σ below the data in all rapidity intervals, integrated over event centrality. Calculations using nuclear PDF sets deviate from the data more strongly than calculations based only on the CT14 NLO PDF set. A similar observation for the CT14 PDF was made in the pp collision system [15] where systematic deviations from the measured values are observed for calculations made at the NNLO. When comparing the measured R AA with calculations, shown in the right panel of Fig. 3, residual deviations from the data are observed. The trend observed in data is consistent with the isospin effect only, expected from the different valence quark content of   [11,14] which favour nPDF sets to describe the data. Fig. 4 shows the normalised Z boson yield as a function of rapidity for three centrality intervals. The results are consistent with each other within their respective statistical uncertainties. The size of the current data sample precludes making a more definitive statement about any possible modification of the Z boson rapidity distribution with centrality. Fig. 5 shows the centrality dependence of the normalised Z boson yield and of R AA compared with results from pp collisions and Glauber MC. The point corresponding to the pp cross section [15] is shown in the plot at N part = 2. The results are derived from Glauber MC v2.4 and a newer version v3.2 following the same procedure as described in Ref. [6]. The results are found to be consistent with each other within experimental uncertainties. The new Glauber MC calculation [61] implements a more advanced treatment of the nuclear density profile and an updated experimental value of the nucleon-nucleon inelastic cross section. From these two updates, the former directly affects the normalised Z boson yield derived in this analysis while the updated cross section value has no appreciable effect. The new model also leads to a set of reduced systematic uncertainties compared with the previous version.
In the estimation of the isospin effect contribution shown with dashed and full lines in Fig. 5, the Glauber MC v3.2 model accounts for the slightly larger radius of the neutron distribution compared to the protons in the Pb nucleus, often called the neutron skin effect. However, due to the weak dependence of Z boson production on the isospin content of the colliding baryons, the predictions of the two Glauber MC versions give essentially the same result.
The normalised yields are consistent with the pp cross section at all measured centralities and show only a weak dependence on N part . The values of R AA , shown in the lower panel, are consistent with unity within the total uncertainty. When the isospin effect is taken into account as shown with the dashed line, the model seems to agree better with the data at low N part values rather than at high values. To quantify the dependence of R AA on N part , the data are fit to a linear function. Including statistical and systematic uncertainties the decrease in the R AA value between the most-peripheral (80-100%) and most-central (0-2%) centrality intervals is found to be (10 ± 7)% and (6 ± 6)% for Glauber MC v2.4 and v3.2, respectively.
The Z boson measurement is used to compare the N part dependence of W [6] and Z boson production by calculating their yield ratios as shown in Fig. 6, where the uncertainties of the two measurements are conservatively treated as uncorrelated. The data points are compared with the ratio measured in pp collisions [15] that is scaled by the isospin factors calculated using the CT14 NLO PDF set. The measurements for both channels are found to be consistent with the scaled pp measurement and show a constant behaviour as a function of centrality.
The trend of the points shown in Fig. 5 for Z bosons is different from the trend observed by the ALICE Collaboration in the measurements with charged hadrons with high transverse momentum [62]. It was recently shown in Ref. [63], that the R AA in peripheral nucleus-nucleus collisions can deviate from unity due to a biased classification of the event geometry for events containing a hard process. In that analysis, the value of R AA without any nuclear effects was determined by using the HG-Pythia model [63], which can create an ensemble of events where the Hijing [64] event generator is used to determine the number of hard sub-interactions for each event, and the particle production is determined solely by superimposing a corresponding number of Pythia 6.4 [40] events. The model is able to qualitatively explain the ALICE measurement in the peripheral region. Fig. 7 shows the R AA for W ± [6] and Z bosons compared with the HG-Pythia model. To compare the model with the measurements of massive electroweak bosons the results are corrected for the isospin effect using the CT14 NLO PDF. All three data sets show trends that are consistent between the species, but that differ from their corresponding model predictions. This suggests that the apparent suppression mechanism [63] that explains the ALICE data [62] for the yields of high-p T charged particles does not have the same effect on the yields of massive electroweak bosons.

Summary
The Z boson production yield per minimum-bias collision, scaled by the mean nuclear thickness function T AA is reported in Pb+Pb collisions at a nucleon-nucleon centre-of-mass energy √ s NN = 5.02 TeV. The measurement is based on data taken by the ATLAS detector at the LHC corresponding to an integrated luminosity of 0.49 nb −1 . Normalised yields are reported in the electron and muon decay channels, differentially in rapidity and collision centrality in the mass window 66 < m < 116 GeV. The fiducial region is defined using the lepton kinematics and detector acceptance. The electron channel and muon channel results are found to agree within the measurement precision and are combined for the final result. The normalised Z boson yields measured in Pb+Pb collisions are 1-3σ higher than NLO pQCD predictions with both free and nuclear PDF sets, where the difference increases towards forward Z boson rapidity. Calculations using nuclear PDF sets deviate from the data more strongly than calculations based only on the CT14 NLO PDF set which is in contrast with recent EW boson measurements performed in the p+Pb system. The nuclear modification factor is measured differentially as a function of Z boson rapidity and event centrality. It is found to be consistent with unity in centrality and to agree with the prediction based on the CT14 PDF set that takes isospin into account. This behaviour is also consistent with the ATLAS measurement performed with W bosons. The yield ratios W /Z are found to be constant as a function of centrality within the uncertainties of the measurements. Unlike high-p T charged hadrons measured by the ALICE Collaboration, W and Z bosons show no indication of yield suppression in peripheral collisions.

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. We