Measurement of the $Z$ boson production cross-section in proton-lead collisions at $\sqrt{s_\mathrm{NN}}=8.16\,\mathrm{TeV}$

This article presents the first measurement of the differential $Z$-boson production cross-section in the forward region using proton-lead collisions with the LHCb detector. The dataset was collected at a nucleon-nucleon centre-of-mass energy of $\sqrt{s_\mathrm{NN}}=8.16\,\mathrm{TeV}$ in 2016, corresponding to an integrated luminosity of $30.8\,\mathrm{nb}^{-1}$. The forward-backward ratio and the nuclear modification factors are measured together with the differential cross-section as functions of the $Z$ boson rapidity in the centre-of-mass frame, the transverse momentum of the $Z$ boson and a geometric variable $\phi^{*}$. The results are in good agreement with the predictions from nuclear parton distribution functions, providing strong constraining power at small Bjorken-$x$.


Introduction
Inclusive Z-boson 1 production at hadron colliders [1] is an important benchmark process to test quantum chromodynamics (QCD), and can be factorised [2] as a product of the hard scattering and the initial state of the collision. State-of-the-art next-to-next-toleading order (NNLO) calculations in perturbative QCD (pQCD), together with next-tonext-to-leading logarithm (NNLL) resummation [3][4][5][6][7][8][9][10] and next-to-leading order (NLO) electroweak (EW) corrections [11][12][13], precisely describe the hard process, while the non-perturbative initial-state can be be modelled through parton distribution functions (PDF) [14,15] or nuclear PDFs (nPDF) [16][17][18][19][20][21]. As a result, Z-boson production at the Large Hadron Collider (LHC) carries valuable information in constraining the PDFs and nPDFs. Weakly coupled Z bosons and their leptonic decay final states have a negligible interaction with the nuclear medium and can be used as clean probes of nuclear-matter effects on the initial-state. The production of Z bosons is sensitive to only the initial-state, while hadronic probes are sensitive to both initial-and final-state nuclear matter effects. Therefore, together with hadronic probes Z-boson production can differentiate between effects of the initial-and final-state. Studies [16] also show that Z-boson production in proton-lead (pPb) collisions at the LHC is sensitive to heavier quark flavours. Improved information on the nuclear corrections is helpful to reduce proton PDF uncertainties and essential for distinguishing the distributions for the different parton flavours. Moreover, Z-boson production is an ideal process to probe transverse-momentum-dependent PDFs (TMDs) [22][23][24][25] when its transverse momentum is below ∼ 10 GeV, where non-perturbative QCD effects start to dominate [26].
The LHCb experiment published for the first time the inclusive Z production result in the forward region in pPb collisions at a nucleon-nucleon centre-of-mass energy √ s NN = 5.02 TeV [27]. The unique forward geometry coverage allows the LHCb detector to probe nPDFs at very small Bjorken-x (10 −4 < x < 10 −3 ). This measurement was followed by the inclusive and differential results in the central region from the CMS and ATLAS experiments [28,29] and in the forward region from the ALICE experiment [30] at √ s NN = 5.02 TeV. Inclusive and differential results are also measured at √ s NN = 8. 16 TeV by the CMS collaboration in the central region [31] and by the ALICE collaboration in the forward region [32]. These results are in good agreement with NLO pQCD predictions calculated with commonly used nPDFs such as EPPS09 [17] and nCTEQ15 [33]. This article presents the first differential measurement of the production cross-section, the forward-backward ratio of the production cross-sections and the nuclear modification factors of Z → µ + µ − production with the LHCb experiment using pPb collisions. The dataset was collected at √ s NN = 8. 16 TeV in 2016. The inclusive cross-sections in the proton-lead (forward) and lead-proton (backward) collisions are measured, together with the differential cross-sections as a function of the rapidity of the Z boson in the centre-ofmass frame (y * Z ), the transverse momentum (p Z T ) and an angular variable ϕ * [34]. The ratio between the cross-sections in forward and backward collisions in a common rapidity range (R FB ) and the nuclear modification factors (R pPb ) with respect to pp collisions are measured, both inclusively and differentially as functions of the above three variables. the free proton PDFs used inside the EPPS16 and nCTEQ15 nPDF sets, respectively. The muon rapidity acceptances in pp, in forward and in backward pPb collisions are different, which need to be corrected for the R FB and R pPb measurements. The PowhegBox generator with the CTEQ6.1 proton PDF is used to derive the corresponding correction factors.

Event selection
The online event selection is performed by a trigger, which consists of a hardware stage followed by a two-level software stage. The hardware trigger used in this analysis selects events containing at least one muon with p T greater than 500 MeV. In the first stage of the software trigger, at least one muon with p T greater than 1.3 GeV and p greater than 6.0 GeV is required. Alignment and calibration of the detector are performed in near real-time [55]. The same alignment and calibration information is propagated to the offline reconstruction, ensuring consistent and high-quality particle identification information between the trigger and the offline software. The present analysis is performed directly using candidates reconstructed at the trigger stage [56,57], owing to the identical performance of the online and offline reconstruction.
The second software-trigger stage and the offline reconstruction level selection require each event to have a PV reconstructed using at least four tracks measured in the vertex detector. For events with multiple PVs, the PV that has the smallest χ 2 IP with respect to the selected dimuon candidate is chosen, where χ 2 IP is defined as the difference between the vertex-fit χ 2 /ndf calculated with and without the two tracks from the dimuon candidate included in the vertex fit. Each identified muon candidate is required to have p T > 20 GeV, 2 < η < 4.5 and a good track fit. The two muon tracks of the Z boson candidate must form a good-quality vertex (vertex-fit χ 2 /ndf < 25), representing a tighter selection compared to the software trigger requirement. The dimuon invariant mass is required to be in the range between 60 and 120 GeV. The selected numbers of candidates are 268 and 166 for forward and backward collisions, respectively. The dimuon invariant mass of the selected candidates is shown in Fig. 1 for both forward and backward collisions. The red histogram shows the distributions of simulated Z → µ + µ − events generated using Pythia 8 with the CTEQ6L1 [40] PDF set, normalized to the number of observed candidates.

Data analysis
This analysis measures the cross-section of the Z → µ + µ − production in a fiducial region with 60 < m µ + µ − < 120 GeV where both muons satisfy p T > 20 GeV and 2 < η < 4.5. The differential cross-section in the fiducial region is defined as where N cand is the number of observed candidates after the selection in the fiducial region; ρ is the purity (the fraction of signal events); f FSR is the correction for final state radiation (FSR); L is the integrated luminosity; ϵ reco&sel , ϵ muon-id and ϵ trig are the efficiencies of reconstruction and selection, muon identification and trigger selection, respectively; and x can be y * Z , p Z T or ϕ * (applied throughout this section unless stated otherwise). The variable ϕ * is defined as where ∆η is the difference between the pseudorapidity of the two muons, ϕ acop = π − |∆ϕ| is the acoplanarity, with ∆ϕ being the difference between the azimuthal angles of the two muons. The angle ϕ * was introduced and studied in Refs. [34,58] as an observable complementary to the variable p Z T to study the shape of the Z boson transverse-momentum distribution with reduced sensitivity to experimental resolution effects.
The ratio of the Z → µ + µ − production cross-sections between forward and backward directions, R FB , is particularly sensitive to nuclear effects. In the case of pp collisions, the R FB is expected to be unity, because the Z boson production cross-section should be symmetric between forward and backward rapidity. However, in the case of pPb collisions, the nuclear modifications are different for forward (small Bjorken-x, i.e., x < 10 −3 ) and backward (large Bjorken-x, i.e., x > 10 −1 ) rapidity. As such, there is suppression at small Bjorken-x due to nuclear shadowing effects [59], and enhancement at large Bjorken-x because of the EMC effect [60]. Therefore, the R FB is an effective observable to probe nuclear-matter effects.
The forward-backward ratio is determined in the common rapidity range 2.5 < |y * Z | < 4.0 as where σ fw is the cross-section in the forward collisions, and σ bw in the backward collisions. The factor k FB (x), which corrects for the difference of the muon rapidity acceptance between forward and backward collisions, is determined as using PowhegBox with the proton PDF CTEQ6.1. Choosing the proton PDF instead of the nPDF avoids the acceptance correction factor being biased by nuclear modification information encoded in the nPDF. The prime symbol ("σ ′ ") in the equation indicates that this cross-section is calculated theoretically instead of from data measurements. The nuclear modification factor (R pPb ) for forward and backward collisions is measured as and respectively, where 208 is the number of binary nucleon-nucleon collisions in pPb collisions. The k pPb factor is to correct for the different y * µ acceptance between pPb and pp collisions and can be similarly calculated using PowhegBox with the proton PDF set CTEQ6.1, and for forward and backward collisions, respectively. By construction, the nuclear modification factor is expected to be unity if no nuclear matter effects are present, and can be used to extract nuclear modifications when a proton or neutron is confined in the Pb nucleus.
The dominant background contribution in this analysis is from QCD processes. This background has contributions from heavy flavour hadrons decaying to muons and charged hadrons mis-identified as muons. Contributions from other sources, such as tt, W + W − and Z → τ + τ − processes, are found to be negligible [61], and are not considered in the present work. Two methods are employed to estimate the QCD background, using pp collision samples weighted to forward or backward pPb collisions according to the multiplicity profiles.
The first method is a "same-sign" technique, based on the property that the yields of positively and negatively charged hadrons are equal in size in the final states. Therefore, the amounts of same-sign and opposite-sign muon-pair candidates that originated from charged hadrons (either mis-identified as muons, or decayed into muons) are the same. By selecting same-sign muon-pair candidates, the amount of QCD background contaminating the selected signal can be determined.
The second method is the ABCD-likelihood technique, see, e.g., Refs. [62][63][64], which is widely used for QCD background estimation in beyond-Standard-Model searches. In this approach two discriminating variables, the vertex-fit χ 2 /ndf and muon isolation variable I charged T , are considered. The vertex-fit χ 2 /ndf is used to verify if the pair of muon candidates originates from the same Z decay, otherwise it is classified as QCD background. The muon isolation variable is defined as the fraction of p T carried by the muon candidate over the sum of p T carried by all the tracks in a cone in the η − ϕ plane with a size of ∆R = ∆η 2 + ∆ϕ 2 < 0.5 around the muon candidate. This variable is used to assess if a muon candidate is surrounded by charged tracks, which characterises a mis-identified hadron or a muon from a charged-hadron decay. A requirement of I charged T > 0.7 can effectively separate the signal from the QCD background [27]. These two variables are used to separate the dimuon sample into four regions: While the ABCD method is more sensitive to QCD backgrounds from heavy-flavour decays, and the same-sign method to hadron mis-identification, there is a sizable overlap between the two methods. Consequently, the final QCD background is taken as the sum of the results from the two methods with a subtraction of this duplication. An estimation of the duplication of the two methods is performed by using the same-sign sample as input to the ABCD method. The resulting average purity is (99.70 ± 0.07)% and (99.75 ± 0.08)% for the forward and backward collisions, respectively. For the differential measurements as functions of y * Z , p Z T , and ϕ * , the background is estimated in the same way and found to be independent of these variables. The average purities given above for forward and backward collisions are used correspondingly for the differential measurements.
The efficiency cannot be directly derived from the pPb dataset because of its small statistical size. The high statistical pp collision dataset is instead chosen as the starting point for efficiency studies. The efficiencies are evaluated in the fiducial region in two steps. In the first step, a tag-and-probe method [65] is used to derive the efficiencies as a function of multiplicity, p µ T , η µ , p Z T and Z boson rapidity in the lab frame (y Z ) from both pp collision data and Z → µ + µ − simulation. The ratios between data and simulation are taken as correction functions. In the second step, the pPb simulation is weighted to the pPb collider data according to the corresponding multiplicity profiles. The correction functions derived in the first step are then used to further weight the pPb simulation to correct for the differences between data and simulation. After a series of reweightings, the pPb simulation can correctly reproduce the efficiencies of the pPb data. The final efficiency to be used for cross-section measurements can be then derived as the fraction of the weighted yields passing the corresponding selection criteria using the weighted pPb simulation.
The resulting average efficiencies (ϵ reco&sel , ϵ muon-id , and ϵ trig ) are shown in Table 1 for the forward and backward collisions. The much smaller ϵ reco&sel efficiency in the backward collisions is mainly due to the far larger contamination from charged tracks when the debris of the Pb nucleus travels towards the LHCb detector. The ϵ reco&sel efficiency for the differential measurements varies from 66.6% to 91.2%, depending on the y * Z , p Z T or ϕ * variables. The efficiencies ϵ muon-id and ϵ trig show negligible dependencies on these variables and are considered as constants.
The detector resolution can lead to migration of events from one interval to another in the distribution of an observable. The migration effect is studied by examining the ratio of the distributions of the y * Z , p Z T and ϕ * observables between generated and reconstructed events using simulation. It is found to be negligible for the observables y * Z and ϕ * , thanks to the excellent angular resolution of the LHCb detector. However, for the variable p Z T the migration effect can be as large as ∼ 10%. An unfolding correction is applied for the p Z T distributions using the RooUnfold [66] software package with an iterative Bayesian approach [67].
Final-state radiation effects can shift certain events out of the kinematic acceptance boundaries. For a fair comparison of the results with theoretical predictions, the FSR correction (f FSR ) is applied to the measured result, which is defined as the ratio of the cross-sections with and without the FSR effect turned on calculated using PowhegBox together with Pythia 8 [38]. The average value of the FSR correction is about 2.5%.
For the nuclear modification factor measurement, the pp reference cross-section at 8.16 TeV is interpolated using a third-order polynomial function from the pp measurements at 7, 8, and 13 TeV [61,[68][69][70] by the LHCb experiment. The resulting interpolated total cross-section at 8. 16 TeV is (98.1 ± 0.6) pb in the rapidity range 2.0 < y * Z < 4.5, and varies little with respect to the measured central value of 95 pb at 8 TeV [68]. The differential cross-section in intervals of y Z , p Z T and ϕ * variables is also interpolated in the same way for each interval. The method has been cross-checked via a linear extrapolation using the 7 and 8 TeV points leading to a negligible difference (< 0.1%) in the result.

Systematic uncertainties
Sources of systematic uncertainties considered are the background estimation, efficiency modelling, FSR correction, detector resolution correction, integrated luminosity, and theoretical corrections.
The uncertainty from the background estimation is less than 0.1%, and includes the statistical uncertainty of the pp collision dataset used in the estimation and the multiplicity weighting. The multiplicity weighting uncertainty is estimated by varying the weighting function within an envelope defined by the statistical fluctuations of the multiplicity profiles of the pp and pPb datasets. The efficiency uncertainty is about 2.5% on average and consists of statistical uncertainties from the pp collision dataset and the simulation samples of forward and backward collisions used in the tag-and-probe method, together Table 1: Observed number of candidates, input values and uncertainties for the overall crosssection and the theory corrections with uncertainties for the R FB and R pPb measurements.

Quantity
Forward 99.69 ± 0.07 99.75 ± 0.08 0.706 ± 0.002 1.518 ± 0.003 with the multiplicity weighting uncertainty. Here the multiplicity weighting uncertainty is estimated in the same way as for the background. The FSR correction uncertainty is estimated to be less than 0.1% as the difference of the default correction factors derived using PowhegBox with Pythia 8 with respect to an alternative approach using Photos [71]. The uncertainty of the pp reference cross-section at 8. 16 TeV is interpolated to be 0.6% on average as an error band enclosed by the upper and lower edges of the error bars of the 7, 8 and 13 TeV measurements. The luminosity uncertainty is about 2.5% [37]. For the differential measurements of the cross-section, R FB , and R pPb as a function of the p Z T , the detector-resolution correction is applied. The uncertainty is given by the difference between the Bayesian approach and other unfolding techniques [72,73], which are considered only in the differential measurements.
For the R FB and R pPb measurements, uncertainties from the y * µ acceptance corrections k FB and k pPb are considered, which include PDF uncertainties, and factorisation and renormalisation scale uncertainties.
The uncertainties of the theoretical predictions shown in the comparisons include the PDF and nPDF uncertainties, and the factorisation and renormalisation scale uncertainties. The (n)PDF uncertainties are estimated by varying the eigenvectors of the (n)PDF sets up and down simultaneously on both the p side and the Pb side. The factorisation and renormalisation scale uncertainties are obtained by varying the scale factors from 0.5 to 2.0. All the theoretical uncertainties are converted to a 68% confidence level when comparing with the measurements. Table 1 summarises the inputs and their uncertainties for the total cross-section, R FB and R pPb measurements. The uncertainties are evaluated separately for R FB and R pPb measurements because of the different y * Z coverage between pp, forward and backward pPb collisions, and are found to be identical to that for the cross-section measurement.
For the differential measurements of the cross-section, R FB , and R pPb as a function of the y * Z , p Z T and ϕ * variables, the uncertainty from the integrated luminosity is considered to be 100% correlated, while the uncertainties from other sources are considered to be partially correlated. Correlation matrices are evaluated for the statistical uncertainty and for the reconstruction and selection efficiency uncertainty, since they form the two largest sources of uncertainty. The statistical uncertainty correlation due to interval migration is evaluated using simulation by comparing the distributions at generator level and reconstruction level. As a result, a 10 -26% correlation can be observed for p Z T , a 2 -3% correlation for y * Z , and a negligibly small correlation for ϕ * . The correlation matrices for the reconstruction and selection efficiencies are calculated by varying the efficiencies in each interval independently within their uncertainties. Large correlations between different intervals are observed. The resulting correlation matrices for statistical uncertainty and the reconstruction and selection efficiency uncertainties are given in Appendix A.

Results and discussion
The total fiducial cross-section of Z → µ + µ − production in forward and backward collisions is measured to be σ fid Z→µ + µ − , pPb = 26.9 ± 1.6 ± 0.9 ± 0.7 nb , σ fid where the first uncertainty is statistical, the second is systematic and the third is due to the uncertainty in the luminosity determination. As a result of the different muon rapidity acceptances between forward and backward collisions, the y * Z coverage is also different. An additional acceptance requirement of 1.5 < y * Z < 4.0 for forward collisions and −4.0 < y * Z < −2.5 for backward collisions is applied on top of the fiducial volume. Figure 2 shows the measured results compared to the PowhegBox calculations using the CTEQ6.1 proton PDF set, EPPS16 nPDF and nCTEQ15 nPDF sets. For forward collisions, the measured value agrees with theoretical predictions with a much smaller uncertainty, which indicates a strong constraining power to the theoretical modelling of the nPDFs. For backward collisions, the measured value is slightly higher but statistically compatible with the theory predictions (difference below 2σ). This trend was also observed in the previous result at √ s NN = 5.02 TeV [27].
The measured differential fiducial cross-sections as a function of y * Z are shown in Figure 3 for forward and backward collisions, together with the PowhegBox calculations with CTEQ6.1, EPPS16 and nCTEQ15 (n)PDF sets. For forward collisions, the measured values show good agreement with the PowhegBox calculations, with a smaller uncertainty for the two intervals of 2.0 < y * Z < 3.0 compared to the theoretical calculations, which can be used to further constrain the nPDFs. For backward collisions, the uncertainty of the measurement is larger than that of the PowhegBox calculation, and the measured central value is higher than the prediction especially for the −3.5 < y * Z < −3.0 interval by about 2σ. However, the measurement and calculation are compatible within uncertainties.
The differential fiducial cross-section as a function of p Z T is shown in Figure 4 (a) and (b), together with the PowhegBox calculations. For forward collisions, the measured values give a smaller uncertainty compared to the PowhegBox calculation for low p Z T intervals, showing a strong constraining power. For backward collisions, the uncertainty of the measurement is larger than that of the PowhegBox calculations but statistically compatible. For p Z T above 30 GeV the central value of the measurement is slightly larger than the PowhegBox calculations, but statistically compatible with them for both forward and backward collisions. The same cross-section measurement as a function of p Z T at low transverse momentum, where non-perturbation effects start to dominate, is also provided with finer intervals, in order to aid theoretical TMD studies, as shown in Figure 4 (c) and (d). The differential fiducial cross-section as a function of the ϕ * variable is shown in Figure 5 together with the PowhegBox calculations. Since the angular observable ϕ * is equivalent to p Z T but not impacted by the momentum resolution, a similar conclusion can be drawn as for p Z T . The detailed numerical values and uncertainties of the differential cross-section as a function of y * Z , p Z T and ϕ * , together with the corresponding f FSR factors are given in Appendix B.
The cross-sections used to determine the forward-backward ratio are different with respect to the fiducial cross-sections. They are measured in the common Z-boson rapidity range 2.5 < |y * Z | < 4.0 with efficiency, purity and FSR corrections also determined in the same rapidity range as shown in Table 1, These corrections are found to be identical to those of the fiducial cross-section measurement. The resulting cross-sections for R FB measurements are 16.1 ± 1.5 nb and 13.4 ± 1.2 nb for forward and backward collisions, respectively. To eliminate the impact from different y * µ acceptances, the correction factor k FB for the overall R FB measurement is calculated to be 0.65 ± 0.02.
The value of R FB is then measured to be R FB = 0.78 ± 0.10, as shown in Figure 6. The measured R FB value is below unity, which is a reflection of the suppression due to,     e.g., nuclear shadowing at small Bjorken-x, together with an average enhancement at large Bjorken-x. The data is in agreement with the EPPS16 and nCTEQ15 predictions. The uncertainty of the measurement is smaller than the theoretical uncertainties using EPPS16 and nCTEQ15 nPDFs, which shows a constraining power on the nPDFs. In addition, the value of R FB is also measured differentially as a function of y * Z , p Z T and ϕ * . The corresponding k FB correction factors for the differential measurements are derived in the same way in separate intervals as shown in Table 2. Figure 6 (b), (c), and (d) show the above measured values of R FB as a function of y * Z , p Z T and ϕ * , respectively, together with the theory calculations. The R FB measurement as a function of y * Z shows a general suppression below unity for all three intervals, in good agreement with the nPDF predictions. For the two intervals in 2.5 < y * Z < 3.5 the measurements give a smaller uncertainty compared to the nCTEQ15 prediction which can be used to constrain the nPDFs. For the R FB measurements as a function of p Z T and ϕ * , a general suppression is expected, however, the significance of the suppression is weak given the large statistical fluctuation in the differential measurements. From theory a slightly stronger suppression at lower p Z T or ϕ * is expected, but given the precision of the current measurement this tendency is not visible. The R FB as a function of ϕ * shows a more stable trend than as a function of p Z T , which exhibits a larger fluctuation The detailed numerical values and uncertainties of the differential measurements of the R FB together with the corresponding f FSR factors are given in Appendix B. For the total nuclear modification factor measurements, because the y * Z coverage between pp and pPb collisions is different, the common y * Z acceptance range 2.0 < y * Z < 4.0 for forward rapidity and −4.0 < y * Z < −2.5 for backward rapidity is used. The y * µ acceptance correction factors are determined to be k fw. pPb = 0.706 ± 0.002 and k bw. pPb = 1.518 ± 0.003 for forward and backward rapidities, respectively. The cross-sections to be used in the R pPb calculations are measured to be 24.2 ± 1.9 nb and 13.4 ± 1.2 nb for forward and backward collisions, respectively, in the corresponding y * Z common acceptance ranges. The variables required to calculate these cross-sections are shown in Table 1. The corresponding pp reference cross-sections at 8. 16 TeV are interpolated in the corresponding y * Z common acceptance ranges as (95.18 ± 0.55) × 10 −3 nb and (79.10 ± 0.51) × 10 −3 nb for forward and backward rapidities, respectively.
The measured overall nuclear modification factors are R fw pPb = 0.94 ± 0.07 and R bw pPb = 1.20 ± 0.11 for forward and backward rapidities, respectively, where the total uncertainties are given. These results are shown in Figure 7 compared with the PowhegBox predictions using the EPPS16 and nCTEQ15 nPDF sets. The overall R pPb results show good compatibility between measurements and theoretical predictions. The backward rapidity result shows larger uncertainty compared to that of the nPDF sets. The measured central value is consistent with the prediction at a 2σ level. The forward rapidity result gives a higher precision than the EPPS16 and nCTEQ15 nPDF sets, and the central value is larger than the prediction, which shows a constraining power on the current nPDF sets.
For the differential measurements of the R pPb , the forward and backward cross-sections and the corresponding pp reference cross-sections are determined in the same way as for the total R pPb measurement but in intervals of y * Z , p Z T and ϕ * . The y * µ acceptance correction factors k pPb for the differential measurements are also derived using PowhegBox with the CTEQ6.1 PDF set considering theory uncertainties from PDF variations and the factorisation and renormalisation scales, as shown in Table 3. The resulting R pPb values as a function of y * Z , p Z T and ϕ * are shown in Figure 8. In general these results are compatible with nPDF predictions. For backward rapidity, larger uncertainties compared to the current nPDF predictions appear for all three observables. However, for forward rapidity the larger dataset gives a higher precision for certain intervals compared to the nPDF predictions. The detailed numerical values and uncertainties of the R pPb as a function of

Conclusion
This article presents a comprehensive analysis of Z → µ + µ − production using pPb collisions at 8. 16 TeV in the forward region recorded with the LHCb detector. The Z-boson production fiducial cross-section, R FB and R pPb are measured inclusively and differentially as a function of y * Z , p Z T , and ϕ * . The results are compatible with theoretical predictions from the EPPS16 and nCTEQ15 nPDFs. The precision of the measurements is significantly improved compared to the previous LHCb results using pPb collisions at a centre-of-mass energy of 5.02 TeV [27] collected in 2013. The measurements provide strong constraints on the nPDFs, especially at high rapidity corresponding to the small Bjorken-x region from 10 −4 to 10 −3 .                 B Differential cross-section, forward-backward ratio, and nuclear modification factors Table 20: Differential cross-section of Z → µ + µ − in intervals of y * Z for forward and backward collisions, together with the FSR correction. In the differential cross-section results, the first uncertainty is statistical, the second is systematic and the third is from integrated luminosity.  Table 21: Differential cross-section of Z → µ + µ − in intervals of p Z T for forward and backward collisions, together with the FSR correction. In the differential cross-section results, the first uncertainty is statistical, the second is systematic and the third is from integrated luminosity.

Appendices A Statistical and systematic correlation matrices
dσ/dp Z T [ nb/ GeV] f FSR Forward Table 22: Differential cross-section of Z → µ + µ − in intervals of p Z T for forward and backward collisions in the low p Z T region, together with the FSR correction. In the differential cross-section results, the first uncertainty is statistical, the second is systematic and the third is from integrated luminosity.
dσ/dp Z T [ nb/ GeV] f FSR Forward Table 23: Differential cross-section of Z → µ + µ − in intervals of ϕ * for forward and backward collisions, together with the FSR correction. In the differential cross-section results, the first uncertainty is statistical, the second is systematic and the third is from integrated luminosity.