Z-boson production in p-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=8.16$ TeV and Pb-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV

Measurement of Z-boson production in p-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=8.16$ TeV and Pb-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV is reported. It is performed in the dimuon decay channel, through the detection of muons with pseudorapidity $-4<\eta_{\mu}<-2.5$ and transverse momentum $p_{\rm T}^{\mu}>20$ GeV/$c$ in the laboratory frame. The invariant yield and nuclear modification factor are measured for opposite-sign dimuons with invariant mass $60<m^{\mu\mu}<120$ GeV$c^2$ and rapidity $2.5<y_{cms}^{\mu\mu}<4$. They are presented as a function of rapidity and, for the Pb-Pb collisions, of centrality as well. The results are compared with theoretical calculations, both with and without nuclear modifications to the Parton Distribution Functions (PDFs). In p-Pb collisions the center-of-mass frame is boosted with respect to the laboratory frame, and the measurements cover the backward ($-4.46<y_{cms}^{\mu\mu}<-2.96$) and forward ($2.03<y_{cms}^{\mu\mu}<3.53$) rapidity regions. For the p-Pb collisions, the results are consistent within experimental and theoretical uncertainties with calculations that include both free-nucleon and nuclear-modified PDFs. For the Pb-Pb collisions, a $3.4\sigma$ deviation is seen in the integrated yield between the data and calculations based on the free-nucleon PDFs, while good agreement is found once nuclear modifications are considered.


Introduction
Measurements of W and Z electroweak vector boson production are useful probes for studying the initial conditions of heavy-ion collisions. Their production occurs predominantly via the Drell-Yan process, in which a quark-antiquark pair annihilates into a lepton pair [1,2]. At leading order, this is an electroweak process although at higher orders gluon radiation must be accounted for. Due to the large masses of these resonances, vector boson production occurs in the early stages of the collisions and their cross section in elementary parton-parton interactions can be calculated within perturbative Quantum Chromodynamics (pQCD). The large masses also allow for high-precision theoretical calculations, currently reaching up to Next-to-Next-to-Leading Order (NNLO) accuracy [3,4].
Since neither the vector bosons nor their leptonic decay products carry color charge, they do not interact strongly with the dense QCD medium formed in heavy-ion collisions.
There are hints that the muons undergo electromagnetic interactions with the dense QCD medium, seen by p T broadening of dimuon spectra [5]. However, this broadening can also be described by photoproduction [6]. Regardless of the physical origin, the scale of these p T -modifications is negligible compared to the average momentum of the muons, which is about half of the Z-boson mass. Thus, the information carried by the muons is not diluted due to final state interactions, allowing to probe the initial state directly. At the Large Hadron Collider (LHC), the center-of-mass energies and luminosities are large enough to allow the production of these bosons to be measured in heavy-ion collisions.
Since the production of electroweak bosons occurs predominantly through quark-antiquark annihilation, it is dependent on the longitudinal momentum distributions of the quarks in the initial state of the collision. These distributions are parametrized by the Parton Distribution Functions (PDFs) f i (x, Q 2 ) [7]. Here, f i gives the probability of finding a parton of type i (this could be either a gluon or a quark with a given flavor) with momentum fraction x of the parent nucleon (also known as Bjorken-x) and squared 4-momentum transfer vector Q 2 . In general, PDFs are obtained through global fits to data, combining information from multiple experiments. Most of these data come from Deep Inelastic Scattering (DIS) experiments, although data from the Tevatron and the LHC have recently been included as well [8][9][10][11].
It has been observed that in a nucleus with mass number A, the distributions of partons f A i (x, Q 2 ) differ from the free-nucleon PDFs scaled by the number of protons and neutrons A f nucleon i (x, Q 2 ) [12]. The modified distributions can be described by means of so-called nuclear Parton Distribution Functions (nPDFs). The Bjorken-x domain can be divided into four main regions, displaying various nuclear effects [13,14]. It should be noted that the precise values of the boundaries of the x-regions depend on the parton flavor, nPDF parametrization and Q 2 . The following values assume Q 2 = M 2 Z [13]. At low x, up to x ∼ 0.05, a depletion of partons is present in nPDFs compared to free-nucleon PDFs. This depletion is referred to as shadowing. Then, in x ∼ 0.05 − 0.3 an enhancement is seen, called antishadowing. Following this, another depletion region, the so-called EMC region, is present from x ∼ 0.3 to x ∼ 0.9. Lastly, x> ∼ 0.9 to unity is the so-called Fermi region, where again an enhancement is present. These nuclear modifications to the PDFs influence the production of electroweak bosons [15], but suffer from large uncertainties. Since the parametrizations are obtained through global fits to data, experimental uncertainties enter the nPDFs as well. Accurate measurements of W and Z bosons at the LHC can therefore help to constrain the nPDFs [13,16,17], which are fundamental ingredients to properly describe the initial state of heavy-ion collisions. An in-depth overview of nPDFs can be found in Ref. [18].
The production of electroweak bosons at the LHC has already been studied in several collision systems, at various energies and rapidities [19][20][21][22][23][24][25][26][27][28][29][30][31][32]. At midrapidity, the data in different collision systems are generally well described by theoretical calculations both with and without nuclear modification of the PDFs [21-31]. In fact, the covered Bjorken-x range at midrapidity extends over the antishadowing and shadowing region (depending on the transition value, even into the EMC region). As a result, their competing effects reduce the final effect of nuclear modifications. However, at larger rapidities and Z-boson production in p-Pb (Pb-Pb) collisions at √ s NN = 8.16 (5.02) TeV ALICE Collaboration therefore lower x, there are increasingly stronger deviations between calculations with models that either do or do not include nuclear modification of the PDFs [20, 31]. The data taken at the LHC are in a kinematic regime in (x, Q 2 ) which is also sensitive to contributions coming from quark flavors such as charm and strange, present as sea quarks in the nucleons [33]. The uncertainties in both the nPDFs and the free-nucleon PDFs (where nuclear fixed-target data were also used for the fits) for these flavors are large and a combination of heavy-ion and proton data can help in reducing them [33].
This paper presents the measurement of the Z-boson production at forward rapidity through the dimuon decay channel in p-Pb collisions at √ s NN = 8. 16 TeV as well as in Pb-Pb collisions at √ s NN = 5.02 TeV. The p-Pb measurement is the first at this energy, following an earlier ALICE publication with data taken at √ s NN = 5.02 TeV [19]. The Z-boson production in Pb-Pb collisions at √ s NN = 5.02 TeV has already been published by the ALICE Collaboration using data collected in 2015 [20]. In 2018, new Pb-Pb data were collected at the same collision energy, corresponding to an integrated luminosity twice that of 2015. The dataset used in this paper includes both the 2015 and 2018 samples and therefore supersedes the previous Pb-Pb results. The larger dataset allows for a more differential analysis as well as increased precision on the integrated cross section measurement.
The paper is organized as follows: the ALICE detector and data samples are detailed in Sec. 2, followed by the analysis procedure in Sec. 3. The main results are then given in Sec. 4 and the conclusions are drawn in Sec. 5.

ALICE detector and data samples
Z bosons are reconstructed via their dimuon decay channel using data from the ALICE muon spectrometer, which selects, identifies and reconstructs muons in the pseudorapidity range −4 < η µ < −2.5 [34]. The tracking system consists of five stations, each containing two multi-wire proportional cathode pad chambers. The third station is located inside a dipole magnet that provides an integrated magnetic field of 3 T × m. A conical absorber of 10 interaction lengths (λ i ) made of carbon, concrete and steel, is located in front of the tracking system to filter out the hadrons and low-momentum muons from the decay of light particles (such as pions or kaons). The muon trigger system consists of four resistive plate chamber planes arranged in two stations placed downstream of an iron wall of ∼7.2 λ i that reduces the contamination of residual hadrons leaking from the front absorber. Finally, a small-angle beam shield made of dense materials protects the whole spectrometer from secondary particles coming from beam-gas interactions and from interactions of large rapidity particles with the beam pipe.
Primary vertex reconstruction is performed by the Silicon Pixel Detector (SPD), the two innermost cylindrical layer of the Inner Tracking System (ITS) [35]. The first and second layer cover the pseudorapidity regions |η| < 2.0 and |η| < 1.4, respectively. Two arrays of scintillator counters (V0A and V0C [36]) are used to trigger events and to reject events from beam-gas interactions. The V0A and V0C detectors are located on both sides of the interaction point at z = 3.4 m and z = −0.9 m and cover the pseudorapidity regions 2.8 < η < 5.1 and −3.7 < η < −1.7, respectively. The V0 detectors are also used to estimate the centrality in Pb-Pb collisions by using a Glauber model fit to the sum of their signal amplitudes [37]. The events are then distributed in classes corresponding to a percentile of the total inelastic hadronic cross section. Finally, the Zero Degree Calorimeters (ZDC) [38], placed on both sides of the interaction point at z = ±112.5 m, are used to reject electromagnetic background. A complete description of the ALICE detector and its performance can be found in Refs. [39,40].
The analysis in p-Pb collisions is performed on data collected in 2016 at a center-of-mass energy √ s NN = 8.16 TeV. The data were taken in two beam configurations, with either the proton (p-going) or lead ion (Pb-going) moving towards the muon spectrometer. By convention, the proton moves toward positive rapidities. Because of the asymmetry in the proton and lead beam energies (E p = 6.5 TeV and E Pb = 2.56 TeV per nucleon), the resulting nucleon-nucleon center-of-mass system is boosted with respect to the √ s NN = 8.16 (5.02) TeV ALICE Collaboration laboratory frame by ∆y cms = ±0.465. Therefore, the rapidity acceptance of the muon spectrometer in the center-of-mass system is 2.03 < y µ µ cms < 3.53 for the p-going configuration and −4.46 < y µ µ cms < −2.96 for the Pb-going configuration. The data used in the Pb-Pb analysis were taken in 2015 and 2018 at √ s NN = 5.02 TeV and cover the rapidity 1 interval 2.5 < y µ µ cms < 4.
The events selected for the analyses require two opposite-sign muon candidates in the muon trigger system, each with a transverse momentum above a configurable threshold, in coincidence with a minimum bias (MB) trigger. The latter was defined by the coincidence of signals in the two arrays of the V0 detector. In the Pb-Pb analysis, only the events corresponding to the most central 90% of the total inelastic cross section (0âȂŞ90%) are used. For these events the MB trigger is fully efficient and the contamination by electromagnetic interactions is negligible. For p-Pb collisions, the Z-boson cross section is calculated using a luminosity normalization factor obtained via a reference process corresponding to the MB trigger condition itself. Therefore the MB trigger efficiency does not affect the cross section evaluation. Finally, the muon trigger threshold was p µ T 0.5 GeV/c for p-Pb and p µ T 1 GeV/c for Pb-Pb collisions. After the event selection, the integrated luminosity in Pb-Pb collisions is about 750 µb −1 . In the p-Pb analysis, where a precise value of the luminosity is needed to compute the Z-boson cross section, dedicated Van der Meer scans were performed [41]. The values of the luminosity amount to 8.40 ± 0.16 nb −1 and 12.74 ± 0.24 nb −1 for the p-going and Pb-going configuration, respectively.

Analysis procedure
The Z-boson signal extraction is performed by combining muons of high transverse momentum in pairs with opposite charge. Muon track candidates are reconstructed in the tracking system of the spectrometer using the algorithm described in Ref. [42]. In order to ensure a clean data sample, a selection is performed on the single muon tracks reconstructed in the muon spectrometer, requiring them to have a pseudorapidity −4 < η µ < −2.5 and a polar angle measured at the end of the front absorber of 170 o < θ abs < 178 o . This procedure removes tracks at the edge of the spectrometer acceptance, and rejects tracks crossing the high-density section of the absorber, which experience significant multiple scattering. The background from tracks not pointing to the nominal interaction vertex, mostly coming from beam-gas interactions and muons produced in the front absorber, is removed by applying a selection on the product of the track momentum and its distance of closest approach to the primary vertex (i.e. the distance to the primary vertex of the track trajectory projected in the plane transverse to the beam axis). Finally, a track is identified as a muon if the track reconstructed in the tracking system matches a track segment in the triggering stations.
Only muons with p µ T > 20 GeV/c are used, to reduce the contribution from low-mass resonances and semileptonic decays of charm and beauty hadrons. The µ + µ − pairs are counted in the invariant mass range 60 < m µ µ < 120 GeV/c 2 , where the Z-boson contribution is dominant with respect to the Drell-Yan process. The invariant mass distributions of the Z-boson candidates are shown in Fig. 1 for minimum bias p-Pb collisions in the p-going and Pb-going configurations, and Pb-Pb collisions in the centrality range 0-90%. Several background sources can contribute to the invariant mass distributions of oppositecharge dimuons. The combinatorial background from random pairing of muons in an event is evaluated by looking at the like-sign pairs (µ ± µ ± ), applying the same selection criteria as for the signal extraction. In the Pb-Pb sample, one pair is found in the invariant mass range considered, which is subtracted from the signal distribution. In p-Pb collisions, no such pairs are found in the region of interest. An upper limit for this background contribution is evaluated by releasing the p µ T selection, fitting the resulting invariant mass distribution between 2 and 50 GeV/c 2 and extrapolating the fit to the 60-120 GeV/c 2 range. Various functions of exponential and power law forms were tried. With this procedure, the number of same-charge events in the region of interest is much smaller than 1% of the opposite-charge √ s NN = 8.16 (5.02) TeV ALICE Collaboration one, and is therefore neglected.
Contributions from cc, bb, tt and the muon decay of τ pairs in the process Z → τ + τ − → µ + µ − + X were estimated with Monte Carlo (MC) simulations using the POWHEG event generator [43]. In p-Pb collisions, the sum of these contributions amounts to 1% of the signal in the p-going configuration, which is taken as a systematic uncertainty from this background source. This contribution is negligible for the Pb-going configuration. In Pb-Pb collisions, a value of 1% is estimated as described in the previous publication [20].
The low amount of background allows the signal to be extracted by counting the candidates in the interval 60 < m µ µ < 120 GeV/c 2 in the distributions shown in Fig. 1. In the p-Pb data sample, 64 ± 8 (34 ± 6) good µ + µ − pairs are counted in the forward (backward) rapidity region. In Pb-Pb collisions, 208 ± 14 Z bosons are counted after merging the 2015 and 2018 data samples. All quoted uncertainties are statistical.  The dimuon invariant mass distributions are compared with the mass shapes obtained by detectorlevel simulations of the Z → µ + µ − process, generated using the POWHEG generator [43] paired with PYTHIA 6.425 [44] for the parton shower. The CT10 [45] free-nucleon PDFs are used, with EPS09NLO [46] √ s NN = 8.16 (5.02) TeV ALICE Collaboration for nuclear modifications. The propagation of the particles through the detector is simulated with the GEANT3 transport code [47]. To account for the modification of the production due to the light-quark flavor content of the nucleus (isospin effect), the simulated distributions are obtained with a weighted average of all possible binary collisions: proton-proton, proton-neutron and for Pb-Pb collisions also neutron-neutron. At high p µ T , tracks are nearly straight so a small misalignment of the detector elements will generate large changes in the track parameters. Therefore, a detailed study of the alignment of the tracking chambers is of utmost importance in order to correctly reproduce the track reconstruction in the simulations. The absolute position of the detector elements was measured by photogrammetry before data taking. The relative position of the elements is then estimated using the MILLEPEDE [48] package, combining data taken with and without magnetic field, with a precision of about 100 µm. This residual misalignment is then taken into account in the simulations of the Z production and the efficiency computation. This method accounts for the relative misalignment of the detector elements but it is not sensitive to a global displacement of the entire muon spectrometer. The simulation of the response of the muon tracking system is based on a data-driven parametrization of the measured resolution of the clusters associated to a track [40], using extended Crystal-Ball (CB) functions [49] to reproduce the distribution of the difference between the cluster and the track positions in each chamber. The CB functions, having a Gaussian core and two power law tails, are first tuned to data and then used to simulate the smearing of the track parameters. The effect of a global misalignment of the spectrometer is implemented by applying a systematic shift, in opposite directions for positive and negative tracks, to the distribution of the angular deviation of the tracks in the magnetic field. This shift is tuned to reproduce the observed difference in the p µ T distributions of positive and negative tracks. In Pb-Pb collisions, the data were taken with two opposite magnetic field polarities of the muon spectrometer dipole magnet. In this case, the sign of the shift is inverted accordingly.
The Z-boson raw yields are corrected for the acceptance times efficiency (A × ε) of the detector. It is evaluated with the MC simulations of the Z → µ + µ − process with POWHEG described above. The A×ε is estimated as the ratio of reconstructed Z bosons with the same selections as for the data, to the number of generated ones with 2.5 < y µ µ lab < 4 for the dimuon pairs, and p µ T > 20 GeV/c and −4 < η µ < −2.5 for the muons. The dimuon invariant mass selection 60 < m µ µ < 120 GeV/c 2 is applied to both reconstructed and generated distributions. In p-Pb collisions, the efficiency is 74% (72%) for the p-going (Pb-going) sample. In Pb-Pb collisions, the efficiency depends on the detector occupancy and therefore on the centrality of the collision. To account for this effect, the generated signal is embedded in real Pb-Pb events. The efficiency is found to be stable from peripheral to semi-central collisions, with a value of about 77% (71%) in the 2015 (2018) data sample and decreases to 71% (66%) for the most central collisions. The centrality-integrated efficiency amounts to 73% in the 2015 dataset and 68% for the 2018 dataset. The Z-boson invariant yield is then computed by dividing the number of measured candidates, corrected for A × ε, by the corresponding number of minimum bias events. The latter is evaluated using the normalization factor F µ−trig/MB , corresponding to the inverse of the probability to observe an opposite-sign dimuon triggered event in a MB event. The value of F µ−trig/MB is evaluated with two methods: (i) by applying the opposite-sign dimuon trigger condition in the analysis of MB events, and (ii) by comparing the counting rate of the two triggers, both corrected for pile-up effects. The final value is the average over the two methods. In Pb-Pb collisions, the normalization factor is computed for all the centrality classes considered.
In the p-Pb analysis, the invariant yield is multiplied by the MB cross section to obtain the Z-production cross section [41]. In the Pb-Pb analysis, results are given both integrated and differential with respect to centrality and rapidity. The production is expressed as the invariant yield, normalized by the nuclear overlap function T AA . The centrality is expressed as N part , the average number of participant nucleons. The T AA and N part quantities are estimated via a Glauber model fit of the signal amplitude in the two arrays of the V0 detector [37]. The nuclear modification of the production of a hard process, such as those producing the Z boson, is measured by R AA , the ratio of the observed normalized yield in √ s NN = 8.16 (5.02) TeV ALICE Collaboration Pb-Pb collisions to that in pp collisions. Due to the insufficient integrated luminosity for pp collisions at √ s NN = 5 TeV, the pp reference is determined from pQCD theoretical calculations using the CT14 PDF set [50].
The relative systematic uncertainties for the p-Pb analysis are summarized in Table 1. The variation between the two methods for the computation of the normalization factor, which is less than 1%, is taken as its systematic uncertainty. The evaluation of A × ε is shown not to be affected by a change of PDF and nPDF set, or transport code in the MC simulations. The uncertainty on the Z-boson yield due to the tracking efficiency, evaluated to be 1% (2%) for the p-going (Pb-going) sample, is obtained by comparing the efficiency between data and MC, using the redundancy of the chambers of the tracking stations [40]. The systematic uncertainty due to the dimuon trigger efficiency is determined by propagating the uncertainty on the efficiency of the detector elements, estimated with a data-driven method based on the redundancy of the trigger chamber information. The matching condition between tracks reconstructed in the tracking and triggering systems introduces a 1% additional uncertainty. Finally, the systematic uncertainty associated to the alignment procedure is evaluated as the difference between the A × ε computed with the data-driven tuning of the cluster resolution, with a global shift, and a MC parametrization without shift. This uncertainty is 7.7% for the p-going dataset and 5.7% for the Pb-going dataset, the difference between the two originating from the difference in the signal shape, which depends on rapidity. The total systematic uncertainty is determined by summing in quadrature the uncertainty from each source.
The sources and values of systematic uncertainties for the Pb-Pb analysis are displayed in Tab. 2. The systematic uncertainties of the normalization factor, the tracking and trigger efficiencies, the trigger/tracker matching, and the alignment are evaluated in the same way as for the p-Pb analysis. The uncertainty of the centrality estimation and the average nuclear overlap function T AA are obtained by varying the centrality class limits by ± 0.5%, as detailed in Ref. [51]. The uncertainty of the theoretical pp cross section, which is used as a reference for the R AA computation, is obtained by varying the factorization and renormalization scales and accounting for the PDF uncertainty. This uncertainty is rapidity dependent and has values between 3.5% and 5.0%. The total systematic uncertainty is taken as the quadratic sum of all the sources.

Results
The production cross section for the Z → µ + µ − process in p-Pb collisions at √ s NN = 8.16 TeV with p µ T > 20 GeV/c and −4 < η µ < −2.5 is measured to be dσ Pb−going Z→µ + µ − /dy = 2.5 ± 0.4 (stat.) ± 0.2 (syst.) nb and dσ p−going Z→µ + µ − /dy = 6.8 ± 0.9 (stat.) ± 0.6 (syst.) nb. In Fig. 2   cms < 3.53) it is roughly between 10 −4 and 10 −3 . The former is expected to be mostly affected by antishadowing and EMC effects, while the latter is in the shadowing region. The observed difference between the backward and forward cross sections is mainly due to the asymmetry of the collision and is consistent with that predicted by theoretical calculations for nucleonnucleon collisions, as shown in the figure. The forward-y region is closer to midrapidity where production cross sections are known to be larger.
The measurements are compared with two model calculations based on pQCD at NLO. The first calculation utilizes the MCFM (Monte Carlo for FeMtobarn processes) code [52] using CT14 at NLO [50] as free-nucleon PDFs. The EPPS16 [53] parametrization of the nuclear modification to the PDFs is then considered to describe the lead environment. The second calculation uses the NNLO code FEWZ (Fully Exclusive W and Z Production) [54]. The lead nucleus is modelled with nCTEQ15 nuclear PDFs [33,55], while CT14 is used for the proton. Both EPPS16 and nCTEQ15 rely on NLO calculations. The latter is a full nPDF set while EPPS16 is anchored to the CT14 free-nucleon PDFs. More details on the approximations and experimental datasets included in the extraction of the nPDFs can be found in Ref. [18]. In all nuclear calculations, the proton and neutron contributions are weighted to reproduce the lead nucleus isospin. Figure 2 shows that the measurements reported here are consistent with pQCD calculations incorporating both free-nucleon and nuclear-modified PDFs, within experimental and theoretical uncertainties. In p-Pb collisions the nuclear effects modify the parton distributions of only one of the two colliding nucleons and the inclusion of the nuclear modification of the PDFs results in a small change if compared to theoretical uncertainties. Moreover, the backward-y region corresponds to a high Bjorken-x range where multiple nuclear effects are present. These lead to both enhancement and depletion compared to free-nucleon PDFs. Their resulting effect is expected to be less pronounced than the one at forward-y where only shadowing is present.
The Z-boson invariant yield normalized by the nuclear overlap function T AA measured in Pb-Pb collisions at √ s NN = 5.02 TeV is 6.1 ± 0.4 (stat.) ± 0.4 (syst.) pb. Because of the symmetry of the collision, the forward rapidity of this measurement probes simultaneously the high-x and low-x range of partons in the lead nucleus. As a result of the rapidity shift and the different nucleon-nucleon center-of-mass energy, these ranges are very close to those probed in p-Pb. In Fig. 3   ized yield is also compared with several pQCD calculations based on different codes (MCFM [52] or FEWZ [54]) and different parton distribution sets. Along with CT14, CT14+EPPS16 and nCTEQ15 calculations [50,53,55], the calculation with CT14 baseline PDF and EPS09s nPDFs is also included in the comparison [56]. Although as a whole EPS09 is superseded by the more recent EPPS16, EPS09s is used because it contains a centrality dependence of the parton distributions which is not provided in the EPPS16 nPDF set. Neutron and proton contributions are properly weighted according to the lead isospin. The uncertainties on the models include the uncertainty on the NLO calculations as well as the uncertainty on the parton distributions that are larger for those including nuclear effects. The large uncertainty of the EPPS16 calculation originates from the larger number of flavor degrees of freedom included in the parametrization [18].
The calculations using nuclear PDFs describe the yield measured in Pb-Pb collisions within uncertainties while the CT14-only calculation deviates from data by 3.4σ . This deviation is not observed in the p-Pb analysis for two main reasons. The first one is statistical. Although the Pb-Pb luminosity is smaller than the p-Pb one, the presence of more nuclear matter in Pb-Pb collisions makes the expected Z-boson yield greater than the one measured in the p-Pb samples, reducing the statistical uncertainty. Second, in Pb-Pb collisions, the distributions of both interacting partons experience nuclear modifications. In order to produce a Z boson at forward rapidity, a collision must occur between a low-x and high-x parton. This leads to a convolution of the shadowing effects at low x and the net nuclear effect observed in backward-y p-Pb collisions. Their combination enhances the suppression of the production with respect to what is separately measured in the two p-Pb rapidity regions.
At the moment most of the nPDF sets do not contain an explicit dependence on the position inside the nucleus, but they provide the average effect over all the nucleons in a given nucleus. Results which are fully inclusive in centrality can better accomodate in the global fitting procedure used to constrain the nPDFs. An estimation of the integrated normalized invariant yield in the 0-100% centrality interval is therefore important. Assuming that the yield scales with the number of nucleon-nucleon binary collisions N coll , and with a conservative estimation of the nuclear modification in the 90-100% centrality interval, the difference between the integrated normalized yields in 0-90% and 0-100% is found to be less than 1 per mille. This means that the present measurement can be regarded also as the normalized Z-boson production in p-Pb (Pb-Pb) collisions at √ s NN = 8.16 (5.02) TeV ALICE Collaboration invariant yield in the 0-100% centrality interval given the current uncertainties.  The Z-boson production in Pb-Pb is also studied as a function of rapidity and centrality. The left panel of Fig. 4 shows the normalized invariant yield in the rapidity intervals 2.5 < y µ µ cms < 2.75, 2.75 < y µ µ cms < 3, 3 < y µ µ cms < 3.25 and 3.25 < y µ µ cms < 4. The results are compared with CT14 predictions both with and without EPPS16 nuclear modification. A shadowing effect is foreseen in the full rapidity range. The right panel of Fig. 4 shows the rapidity dependence of the nuclear modification factor R AA , computed by dividing the yield normalized to T AA by the pp cross section at √ s = 5.02 TeV obtained with pQCD calculations with CT14 PDFs. For this observable, the uncertainties on the free-nucleon PDFs are factored out in the theoretical calculations, and the remaining uncorrelated uncertainties are related to the nuclear PDFs only. The measured R AA is in agreement within uncertainties with the EPPS16 calculations while, at large rapidity, it deviates from the free-nucleon calculations.
In Fig. 5, the normalized invariant yield is shown as a function of centrality. The CT14 calculations are based on free-nucleon PDFs and therefore, by construction, carry no centrality dependence. The data are also compared with calculations from EPS09s [56], which show a decrease in the invariant yield towards more central collisions, although the effect is very weak. Furthermore, in each centrality bin the EPS09s prediction is consistent with the more recent EPPS16 set, which does not implement a dependence on the impact parameter (the CT14+EPPS16 calculation is displayed in Fig. 3).
Within uncertainties, each data point is well described by models including nPDFs, while the CT14-only calculation overestimates the data, especially for the most central collisions where the difference is 3.9σ .

Conclusions
The Z-boson production has been studied at large rapidities in p-Pb collisions at √ s NN = 8.16 TeV and in Pb-Pb collisions at √ s NN = 5.02 TeV.
For the p-Pb collisions, the Z bosons were measured in the rapidity range −4.46 < y µ µ cms < −2.96 and 2.03 < y µ µ cms < 3.53. The production cross section at forward and backward rapidity has been compared with theoretical predictions, both with and without nuclear modifications. The data show little sensitivity Z-boson production in p-Pb (Pb-Pb) collisions at √ s NN = 8.16 (5.   The results are compared with the free-nucleon PDF prediction (CT14 [50]) and with calculations with the centrality-dependent EPS09s nPDFs [56].
to the presence of nuclear effects, partially because in p-Pb collisions, nuclear modifications to the PDFs affect only one of the two colliding particles. This is particularly true in the backward region, where enhancement and depletion effects on nPDFs tend to compensate. As a result, the calculations for the nuclear modification of the PDFs are very close to those without. In the forward region, low-x partons of the Pb nucleus are probed which are only sensitive to shadowing (which corresponds to a depletion in the nPDFs). Consequently, nuclear effects tend more clearly to induce a decrease in the cross section. Nonetheless it remains compatible within uncertainties with the one calculated neglecting such effects.
In the Pb-Pb data, the invariant yield normalized by the average nuclear overlap function has been Z-boson production in p-Pb (Pb-Pb) collisions at √ s NN = 8.16 (5.02) TeV ALICE Collaboration evaluated in the rapidity range from 2.5 < y µ µ cms < 4 and in the 0-90% centrality class. The results obtained in this paper supersede those from an earlier ALICE publication [20], where only part of the current dataset was used. The experimental data are, within uncertainties, in agreement with theoretical calculations that include various parametrizations of nuclear modification of the PDFs. On the contrary, the integrated yield deviates by 3.4σ from the prediction obtained using free-nucleon PDFs.
Comparisons with models of the measured differential yields versus centrality and rapidity were also carried out, generally showing agreement with nuclear modified PDFs. In contrast, a discrepancy with calculations based on free-nucleon PDFs was found. The differential measurements presented in this paper can provide additional constraints to the nPDFs.
The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The ALICE Collaboration gratefully acknowledges the resources and support provided by all Grid centres and the Worldwide LHC Computing Grid (WLCG) collaboration. The ALICE Collaboration acknowledges the following funding agencies for their support in building and running the ALICE detector:       [