EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

 CERN-EP-2022-089 LHCb-PAPER-2022-009 20 May 2022

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

LHCb collaboration

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

Submitted to JHEP

© 2024 CERN for the benefit of the LHCb collaboration. CC BY 4.0 licence.

##  1 Introduction

Inclusive $Z$-boson111In this paper, the label $Z$ boson is defined to include
contributions from virtual photons. 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-to-leading order (NNLO)
calculations in perturbative QCD (pQCD), together with next-to-next-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 with a well constrained final-
state, while hadronic probes are sensitive to both the initial- and final-
state 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 (${p}\mathrm{Pb}$)
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 much smaller than its mass [26].

The LHCb experiment published for the first time the inclusive $Z$ production
result in the forward region in ${p}\mathrm{Pb}$ collisions at a nucleon-
nucleon centre-of-mass energy
$\sqrt{s_{\scriptscriptstyle\text{NN}}}=5.02\text{\,Te\kern-1.00006ptV}$ [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,
ATLAS and ALICE experiments [28, 29, 30]. The ALICE and CMS collaborations
also presented inclusive and differential results in the central region [31,
32] at
$\sqrt{s_{\scriptscriptstyle\text{NN}}}=8.16\text{\,Te\kern-1.00006ptV}$.
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}\\!\rightarrow{\mu^{+}\mu^{-}}$
production with the LHCb experiment using ${p}\mathrm{Pb}$ collisions. The
dataset was collected at
$\sqrt{s_{\scriptscriptstyle\text{NN}}}=8.16\text{\,Te\kern-1.00006ptV}$ 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-
of-mass frame ($y^{*}_{{Z}}$), the transverse momentum
($p_{\mathrm{T}}^{{Z}}$) and an angular variable $\phi^{*}$ [34]. The ratio
between the cross-sections in forward and backward collisions in a common
rapidity range ($R_{\mathrm{FB}}$) and the nuclear modification factors
($R_{{p}\mathrm{Pb}}$) with respect to $pp$ collisions are measured, both
inclusively and differentially as functions of the above three variables.

##  2 Detector, datasets, and theoretical predictions

The LHCb detector [35, 36] is a single-arm forward spectrometer covering the
pseudorapidity ($\eta$) range from 2 to 5. The detector includes a high-
precision tracking system consisting of a silicon-strip vertex detector
surrounding the $pp$ interaction region, a large-area silicon-strip detector
located upstream of a dipole magnet with a bending power of about
$4{\mathrm{\,Tm}}$, and three stations of silicon-strip detectors and straw
drift tubes placed downstream of the magnet. The tracking system provides a
measurement of the momentum, $p$, of charged particles with a relative
uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV. Natural
units with $\hbar=c=1$ are used throughout. The impact parameter (IP), defined
as the minimum distance of a track to a primary vertex (PV), is measured with
a resolution of $(15+29/p_{\mathrm{T}})\,\upmu\text{m}$, where
$p_{\mathrm{T}}$ is the component of the momentum transverse to the beam, in
GeV. Muons are identified by a system composed of alternating layers of iron
and multiwire proportional chambers.

This analysis uses ${p}\mathrm{Pb}$ collision data at a centre-of-mass energy
per nucleon pair of $\sqrt{s_{\scriptscriptstyle\text{NN}}}$ =8.16 TeV
recorded in 2016 by the LHCb detector, with an energy of
$6.5\text{\,Te\kern-1.00006ptV}$ for the $p$ beam and
$2.56\text{\,Te\kern-1.00006ptV}$ per nucleon for the Pb beam. Since the
collision system is asymmetric, a detector acceptance of $2.0<\eta<4.5$ in the
laboratory frame is translated into a muon rapidity ($y^{*}_{\mu}$) acceptance
of $1.53<y^{*}_{\mu}<4.03$ in the centre-of-mass frame 222When the momentum of
a particle is much greater than its mass, its rapidity is approximately equal
to its pseudorapidity ($y\sim\eta$)., when the $p$ beam travels towards the
positive $\eta$ direction of the LHCb detector (forward collision). Similarly,
when the $p$ beam travels in the negative $\eta$ direction (backward
collision), the muon rapidity acceptance in the centre-of-mass frame is
$-4.97<y^{*}_{\mu}<-2.47$. The integrated luminosity is $12.2\pm
0.3\text{\,nb}^{-1}$ for forward collisions and $18.6\pm 0.5\text{\,nb}^{-1}$
for backward collisions. For the background and efficiency estimation, $pp$
collisions recorded in 2016 at $\sqrt{s}$ =13 TeV are used, corresponding to
an integrated luminosity of 1.7$\text{\,fb}^{-1}$.

Full detector simulation is required to model the detector acceptance effects
and selection efficiencies. The ${p}\mathrm{Pb}$ collisions are simulated by
embedding multiple minimum-bias events generated separately into a
${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ signal event, so that the multiplicity
profile in the simulated events agrees with that of the collision data. The
signal ${Z}\\!\rightarrow{\mu^{+}\mu^{-}}$ events are generated with Pythia 8
[37, 38] using a specific LHCb configuration and the CTEQ6L1 [39] PDF set with
a specific LHCb configuration [40] assuming $pp$ interactions with beam
momenta equal to the momenta per nucleon of the $p$ and Pb beams. Decays of
unstable particles are described by EvtGen [41], in which final-state
radiation is generated using Photos [42]. Minimum-bias events are generated
using the Epos event generator with the LHC model [43]. The interaction and
the response of the generated particles with the detector are implemented
using the Geant4 toolkit [44, 45], as described in Ref. [46].

The results of the analysis are compared to theory predictions determined at
fixed order in pQCD and given by the PowhegBox v2 [47, 48, 49, 50] generator,
which calculates the hard interaction using pQCD at NLO. The PowhegBox package
is interfaced with the Lhapdf package [51] to supply different (n)PDF sets for
describing the non-perturbative initial states. PDF sets CT14 [52] and CTEQ6.1
[53] are used for the $p$ side of the ${p}\mathrm{Pb}$ collisions, and nPDF
sets EPPS16 [17, 52] and nCTEQ15 [33] are used for the Pb side. Choosing CT14
and CTEQ6.1 PDF sets for the $p$ side is motivated by the fact that they are
the free proton PDFs used inside the EPPS16 and nCTEQ15 nPDF sets,
respectively. The muon rapidity acceptances in $pp$, in forward and in
backward ${p}\mathrm{Pb}$ collisions are different, which need to be corrected
for the $R_{\mathrm{FB}}$ and $R_{{p}\mathrm{Pb}}$ measurements. The PowhegBox
generator with the CTEQ6.1 proton PDF is used to derive the corresponding
correction factors.

##  3 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_{\mathrm{T}}$ greater than 500 MeV. In the first stage of the software
trigger, at least one muon with $p_{\mathrm{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 [54]. 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 [55, 56], 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 $\chi^{2}_{\text{IP}}$ with respect to the selected
dimuon candidate is chosen, where $\chi^{2}_{\text{IP}}$ is defined as the
difference between the vertex-fit $\chi^{2}$ 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_{\mathrm{T}}>20\text{\,Ge\kern-1.00006ptV}$, $2<\eta<4.5$ and a good track
fit. The two muon tracks of the $Z$ boson candidate must form a good-quality
vertex (vertex-fit $\chi^{2}<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 the
forward and backward collisions. The red line shows the distributions from
simulation generated using Pythia 8 with the CTEQ6L1 [39] PDF set, normalized
to the number of observed candidates.

(a)(b)

Figure 1: Dimuon invariant mass of the selected events for (a) forward and (b)
backward collisions, respectively. The red histograms show the distributions
from simulation normalised to the number of observed candidates.

##  4 Data analysis

This analysis measures the cross-section of the
${Z}\rightarrow{\mu^{+}\mu^{-}}$ production in a fiducial region with
$60<m_{{\mu^{+}\mu^{-}}}<120\text{\,Ge\kern-1.00006ptV}$ where both muons
satisfy $p_{\mathrm{T}}>20~{}\text{\,Ge\kern-1.00006ptV}$ and $2<\eta<4.5$.
The differential cross-section in the fiducial region is defined as

$\frac{\mathrm{d}\sigma^{\mathrm{fid}}_{{Z}\rightarrow{\mu^{+}\mu^{-}}}}{\mathrm{d}x}=\frac{\rho(x)\cdot
f_{\mathrm{FSR}}}{\mathcal{L}\cdot\epsilon^{\mathrm{reco\&sel}}(x)\cdot\epsilon^{\mathrm{muon{\text{-}}id}}(x)\cdot\epsilon^{\mathrm{trig}}(x)}\cdot\frac{\mathrm{d}N_{\mathrm{cand}}}{\mathrm{d}x}\,,$
(1)

where $N_{\mathrm{cand}}$ is the number of observed candidates after the
selection in the fiducial region; $\rho$ is the purity (the fraction of signal
events); $f_{\mathrm{FSR}}$ is the correction for final state radiation (FSR);
$\mathcal{L}$ is the integrated luminosity; $\epsilon^{\mathrm{reco\&sel}}$,
$\epsilon^{\mathrm{muon{\text{-}}id}}$ and $\epsilon^{\mathrm{trig}}$ are the
efficiencies of reconstruction and selection, muon identification and trigger
selection, respectively; and $x$ can be $y^{*}_{{Z}}$, $p_{\mathrm{T}}^{{Z}}$
or $\phi^{*}$ (applied throughout this section unless stated otherwise). The
variable $\phi^{*}$ is defined as

$\phi^{*}=\frac{\tan(\phi_{\mathrm{acop}}/2)}{\cosh(\Delta\eta/2)},$ (2)

where $\Delta\eta$ is the difference between the pseudorapidity of the two
muons, $\phi_{\mathrm{acop}}=\pi-|\Delta\phi|$ is the acoplanarity, with
$\Delta\phi$ being the difference between the azimuthal angles of the two
muons. The angle $\phi^{*}$ was first introduced in Ref. [34] as an observable
equivalent to $p_{\mathrm{T}}^{{Z}}$ to study the shape of the $Z$ boson
transverse-momentum distribution, but without being affected by the momentum
resolution of the detector.

The ratio of the ${Z}\rightarrow{\mu^{+}\mu^{-}}$ production cross-sections
between forward and backward directions, $R_{\mathrm{FB}}$, is particularly
sensitive to nuclear effects. In the case of $pp$ collisions, the
$R_{\mathrm{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 ${p}\mathrm{Pb}$ collisions, the nuclear modifications
are different for forward (small Bjorken-$x$) and backward (large Bjorken-$x$)
rapidity. As such, there is suppression at small Bjorken-$x$ due to nuclear
shadowing effects [57], and enhancement at large Bjorken-$x$ because of the
EMC effect [58]. Therefore, the $R_{\mathrm{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

$R_{\mathrm{FB}}(x)=\left.\frac{\mathrm{d}\sigma_{\mathrm{fw}}/\mathrm{d}x}{\mathrm{d}\sigma_{\mathrm{bw}}/\mathrm{d}x}\cdot
k_{\mathrm{FB}}(x)\right|_{2.5<|y^{*}_{{Z}}|<4.0}~{},$ (3)

where $\sigma_{\mathrm{fw}}$ is the cross-section in the forward collisions,
and $\sigma_{\mathrm{bw}}$ in the backward collisions. The factor
$k_{\mathrm{FB}}(x)$, which corrects for the difference of the muon rapidity
acceptance between forward and backward collisions, is determined as

$k_{\mathrm{FB}}(x)=\dfrac{\mathrm{d}\sigma^{\prime}_{(pp,-4.97<y^{*}_{\mu}<-2.47)}/\mathrm{d}x}{\mathrm{d}\sigma^{\prime}_{(pp,1.53<y^{*}_{\mu}<4.03)}/\mathrm{d}x}\Bigg{|}_{2.5<|y^{*}_{{Z}}|<4.0}~{},$
(4)

using PowhegBox with the proton PDF CTEQ6.1, without including any nuclear
modification information to preserve the nuclear matter effects from the data.
The prime symbol (“$\sigma^{\prime}$”) in the equation indicates that this
cross-section is calculated theoretically instead of from data measurements.

The nuclear modification factor ($R_{{p}\mathrm{Pb}}$) for forward and
backward collisions is measured as

$R_{{p}\mathrm{Pb}}^{\rm fw}(x)=k_{{p}\mathrm{Pb}}^{\rm
fw}(x)\cdot\dfrac{\mathrm{d}\sigma_{({p}\mathrm{Pb},1.53<y^{*}_{\mu}<4.03)}/\mathrm{d}x}{208\cdot\mathrm{d}\sigma_{(pp,2.0<y^{*}_{\mu}<4.5)}/\mathrm{d}x}$
(5)

and

$R_{{p}\mathrm{Pb}}^{\rm bw}(x)=k_{{p}\mathrm{Pb}}^{\rm
bw}(x)\cdot\dfrac{\mathrm{d}\sigma_{({p}\mathrm{Pb},-4.97<y^{*}_{\mu}<-2.47)}/\mathrm{d}x}{208\cdot\mathrm{d}\sigma_{(pp,-4.5<y^{*}_{\mu}<-2.0)}/\mathrm{d}x}~{},$
(6)

respectively, where 208 is the number of binary nucleon-nucleon collisions in
${p}\mathrm{Pb}$ collisions. The $k_{{p}\mathrm{Pb}}$ factor is to correct for
the different $y^{*}_{\mu}$ acceptance between ${p}\mathrm{Pb}$ and $pp$
collisions and can be similarly calculated using PowhegBox with the proton PDF
set CTEQ6.1,

$k_{{p}\mathrm{Pb}}^{\rm
fw}(x)=\dfrac{\mathrm{d}\sigma^{\prime}_{(pp,2.0<y^{*}_{\mu}<4.5)}/\mathrm{d}x}{\mathrm{d}\sigma^{\prime}_{(pp,1.53<y^{*}_{\mu}<4.03)}/\mathrm{d}x}$
(7)

and

$k_{{p}\mathrm{Pb}}^{\rm
bw}(x)=\dfrac{\mathrm{d}\sigma^{\prime}_{(pp,-4.5<y^{*}_{\mu}<-2.0)}/\mathrm{d}x}{\mathrm{d}\sigma^{\prime}_{(pp,-4.97<y^{*}_{\mu}<-2.47)}/\mathrm{d}x}$
(8)

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 ${t}{\overline{t}}$, $W^{+}$ $W^{-}$ and
${Z}\\!\rightarrow{\tau^{+}\tau^{-}}$ processes, are found to be negligible
[59], 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 ${p}\mathrm{Pb}$ 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., Ref. [60],
which is widely used for QCD background estimation in beyond-Standard-Model
searches. In this approach two discriminating variables, the vertex-fit
$\chi^{2}$ and muon isolation variable $I_{\mathrm{T}}^{\mathrm{charged}}$,
are considered. The vertex-fit $\chi^{2}$ 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_{\mathrm{T}}$ carried by the muon candidate over the sum of
$p_{\mathrm{T}}$ carried by all the tracks in a cone in the $\eta-\phi$ plane
with a size of $\Delta R=\sqrt{{\Delta\eta}^{2}+{\Delta\phi}^{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. These two variables are used to separate the
dimuon sample into four regions:

  * (A)

vertex-fit $\chi^{2}<25$ and at least one muon with
$I_{\mathrm{T}}^{\mathrm{charged}}>0.7$,

  * (B)

vertex-fit $\chi^{2}>70$ and at least one muon with
$I_{\mathrm{T}}^{\mathrm{charged}}>0.7$,

  * (C)

vertex-fit $\chi^{2}<25$ and both muons with
$I_{\mathrm{T}}^{\mathrm{charged}}<0.7$, and

  * (D)

vertex-fit $\chi^{2}>70$ and both muons with
$I_{\mathrm{T}}^{\mathrm{charged}}<0.7$,

where ($\mathrm{A+C}$) is the signal selection fiducial volume. The background
contamination in the region ($\mathrm{A+C}$) is therefore estimated as
$N_{\mathrm{bkg}}(\mathrm{A+C})=N(\mathrm{C})/N(\mathrm{D})\cdot
N(\mathrm{B+D})$.

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\pm 0.07)\%$ and $(99.75\pm
0.08)\%$ for the forward and backward collisions, respectively. For the
differential measurements as functions of $y^{*}_{{Z}}$,
$p_{\mathrm{T}}^{{Z}}$, and $\phi^{*}$, 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 efficiencies are evaluated in the fiducial region in two steps. In the
first step, a tag-and-probe method [61] is used to derive the efficiencies as
a function of multiplicity, $p_{\mathrm{T}}^{\mu}$, ${\eta}^{\mu}$,
$p_{\mathrm{T}}^{{Z}}$ and $Z$ boson rapidity in the lab frame ($y_{{Z}}$)
from both $pp$ collision data and ${Z}\rightarrow{\mu^{+}\mu^{-}}$ simulation.
The ratios between data and simulation are taken as correction functions. In
the second step, the ${p}\mathrm{Pb}$ simulation is weighted to the
${p}\mathrm{Pb}$ collider data according to the corresponding multiplicity
profiles. The correction functions derived in the first step are then used to
further weight the ${p}\mathrm{Pb}$ simulation to correct for the differences
between data and simulation. After a series of reweightings, the
${p}\mathrm{Pb}$ simulation can correctly reproduce the efficiencies of the
${p}\mathrm{Pb}$ 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
${p}\mathrm{Pb}$ simulation.

Table 1: Observed number of candidates, input values and uncertainties for the overall cross-section and the theory corrections with uncertainties for the $R_{\mathrm{FB}}$ and $R_{{p}\mathrm{Pb}}$ measurements. Quantity | Forward | Backward
---|---|---
$N_{\text{cand}}$ (for $\sigma^{\mathrm{fid}}$) | 268 | 166
$N_{\text{cand}}$ (for $R_{\mathrm{FB}}$) | $160$ | $166$
$N_{\text{cand}}$ (for $R_{{p}\mathrm{Pb}}$) | $241$ | $166$
$\rho\,[\%]$ | $99.69\pm 0.07$ | $99.75\pm 0.08$
$\epsilon^{\mathrm{reco\&sel}}$ [%] | $87.2\phantom{0}\pm 2.9$ | $72.0\phantom{0}\pm 2.5$
$\epsilon^{\mathrm{muon{\text{-}}id}}$ [%] | $97.3\phantom{0}\pm 0.3$ | $97.3\phantom{0}\pm 0.3$
$\epsilon^{\mathrm{trig}}$ [%] | $98.3\phantom{0}\pm 0.6$ | $97.1\phantom{0}\pm 0.6$
$\mathcal{L}$ $[\text{\,nb}^{-1}]$ | $12.2\phantom{0}\pm 0.3$ | $18.6\phantom{0}\pm 0.5$
$f_{\mathrm{FSR}}$ | $1.02\phantom{0}\pm 0.01$ | $1.02\phantom{0}\pm 0.01$
$k_{\mathrm{FB}}$ (for $R_{\mathrm{FB}}$) | $0.65\phantom{0}\pm 0.02$ |  00–
$k_{{p}\mathrm{Pb}}$ (for $R_{{p}\mathrm{Pb}}$) | $0.706\pm 0.002$ | $1.518\pm 0.003$

The resulting average efficiencies ($\epsilon^{\mathrm{reco\&sel}}$,
$\epsilon^{\mathrm{muon{\text{-}}id}}$, and $\epsilon^{\mathrm{trig}}$) are
shown in Table 1 for the forward and backward collisions. The much smaller
$\epsilon^{\mathrm{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
$\epsilon^{\mathrm{reco\&sel}}$ efficiency for the differential measurements
varies from $(66.6\pm 4.1)\%$ to $(91.2\pm 2.9)\%$, depending on the
$y^{*}_{{Z}}$, $p_{\mathrm{T}}^{{Z}}$ or $\phi^{*}$ variables. The
efficiencies $\epsilon^{\mathrm{muon{\text{-}}id}}$ and
$\epsilon^{\mathrm{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_{\mathrm{T}}^{{Z}}$ and $\phi^{*}$ observables between generated and
reconstructed events using simulation. It is found to be negligible for the
observables $y^{*}_{{Z}}$ and $\phi^{*}$, thanks to the excellent angular
resolution of the LHCb detector. However, for the variable
$p_{\mathrm{T}}^{{Z}}$ the migration effect can be as large as $\sim 10\%$. An
unfolding correction is applied for the $p_{\mathrm{T}}^{{Z}}$ distributions
using the RooUnfold [62] software package with an iterative Bayesian approach
[63].

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_{\mathrm{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 [37]. 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\text{\,Te\kern-1.00006ptV}$ is interpolated using a third-
order polynomial function from the $pp$ measurements at 7, 8, and 13 TeV [59,
64, 65] by the LHCb experiment. The resulting interpolated total cross-section
at 8.16 TeV is $(98.1\pm 0.6)\times 10^{-3}\text{\,nb}$ in the rapidity range
$2.0<y^{*}_{{Z}}<4.5$. The differential cross-section in intervals of
$y_{{Z}}$, $p_{\mathrm{T}}^{{Z}}$ and $\phi^{*}$ 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.

##  5 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 weighitng
uncertainty is estimated by varying the weighitng function within an envelope
defined by the statistical fluctuations of the multiplicity profiles of the
$pp$ and ${p}\mathrm{Pb}$ 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 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
correction factors between using the EPPS16 and nCTEQ15 nPDF sets. The
uncertainty of the $pp$ reference cross-section at 8.16 TeV is interpolated 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\%$. For
the differential measurements of the cross-section, $R_{\mathrm{FB}}$, and
$R_{{p}\mathrm{Pb}}$ as a function of the $p_{\mathrm{T}}^{{Z}}$, the
detector-resolution correction is applied. The uncertainty is given by the
difference between the Bayesian approach and other unfolding techniques [66,
67], which are considered only in the differential measurements.

The uncertainties of the theoretical predictions shown in the comparisons with
experimental data include the same uncertainties as those considered through
the $y^{*}_{\mu}$ acceptance corrections in the data measurements of
$R_{\mathrm{FB}}$ and $R_{{p}\mathrm{Pb}}$. In addition, the PDF and nPDF
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_{\mathrm{FB}}$ and $R_{{p}\mathrm{Pb}}$ measurements. The
uncertainties are evaluated separately for $R_{\mathrm{FB}}$ and
$R_{{p}\mathrm{Pb}}$ measurements because of the different $y^{*}_{{Z}}$
coverage between $pp$, forward and backward ${p}\mathrm{Pb}$ collisions, and
are found to be identical to that for the cross-section measurement.

For the differential measurements of the cross-section, $R_{\mathrm{FB}}$, and
$R_{{p}\mathrm{Pb}}$ as a function of the $y^{*}_{{Z}}$,
$p_{\mathrm{T}}^{{Z}}$ and $\phi^{*}$ 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_{\mathrm{T}}^{{Z}}$, a $2$ – $3\%$
correlation for $y^{*}_{{Z}}$, and a negligibly small correlation for
$\phi^{*}$. 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.

##  6 Results and discussion

The total fiducial cross-section of ${Z}\rightarrow{\mu^{+}\mu^{-}}$
production in forward and backward collisions is measured to be

$\displaystyle\sigma^{\mathrm{fid}}_{{Z}\rightarrow{\mu^{+}\mu^{-}},~{}{p}\mathrm{Pb}}$
$\displaystyle=~{}26.9\pm 1.6\pm 0.9\pm 0.7~{}\text{\,nb}~{},$
$\displaystyle\sigma^{\mathrm{fid}}_{{Z}\rightarrow{\mu^{+}\mu^{-}},~{}\mathrm{Pb}{p}}$
$\displaystyle=~{}13.4\pm 1.0\pm 0.5\pm 0.3~{}\text{\,nb}~{},$

where the first uncertainty is statistical, the second is systematic and the
third is due to the uncertainty in the luminosity determination. Because of
the different $y^{*}_{{Z}}$ acceptances between forward and backward
collisions, 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. This trend was also observed in the previous result at
$\sqrt{s_{\scriptscriptstyle\text{NN}}}=5.02\text{\,Te\kern-1.00006ptV}$ [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. However, the measurement
and calculation are compatible within uncertainties.

The differential fiducial cross-section as a function of
$p_{\mathrm{T}}^{{Z}}$ 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_{\mathrm{T}}^{{Z}}$ 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_{\mathrm{T}}^{{Z}}$ above $30\text{\,Ge\kern-1.00006ptV}$ 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_{\mathrm{T}}^{{Z}}$ at
low transverse momentum with finer intervals is provided, 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 $\phi^{*}$ variable is shown in
Figure 5 together with the PowhegBox calculations. Since the angular
observable $\phi^{*}$ is equivalent to $p_{\mathrm{T}}^{{Z}}$ but not impacted
by the momentum resolution, a similar conclusion can be drawn as for
$p_{\mathrm{T}}^{{Z}}$. The detailed numerical values and uncertainties of the
differential cross-section as a function of $y^{*}_{{Z}}$,
$p_{\mathrm{T}}^{{Z}}$ and $\phi^{*}$, together with the corresponding
$f_{\mathrm{FSR}}$ factors are given in Appendix B.

Figure 2:  The measured overall ${Z}\rightarrow{\mu^{+}\mu^{-}}$ production
fiducial cross-section compared to the PowhegBox prediction using CTEQ6.1,
EPPS16 and nCTEQ15 (n)PDF sets, for forward and backward collisions,
respectively.

(a)(b)

Figure 3:  The measured differential fiducial cross-section as a function of
$y^{*}_{{Z}}$ for (a) forward and (b) backward collisions. The theoretical
predictions are calculated using PowhegBox with CTEQ6.1, EPPS16 and nCTEQ15
(n)PDF sets.

(a)(b)
(c)(d)

Figure 4:  The measured differential fiducial cross-section as a function of
$p_{\mathrm{T}}^{{Z}}$ for (a) forward and (b) backward collisions, and the
corresponding version for (c) forward and (d) backward collisions with fine
intervals at low $p_{\mathrm{T}}^{{Z}}$. The theoretical predictions are
calculated using PowhegBox with CTEQ6.1, EPPS16 and nCTEQ15 (n)PDF sets.

(a)(b)

Figure 5:  The measured differential fiducial cross-section as a function of
$\phi^{*}$ for (a) forward and (b) backward collisions. The theoretical
predictions are calculated using PowhegBox with CTEQ6.1, EPPS16 and nCTEQ15
(n)PDF sets.

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_{\mathrm{FB}}$
measurements are $15.7\pm 1.4\text{\,nb}$ and $13.1\pm 1.2\text{\,nb}$ for
forward and backward collisions, respectively. To eliminate the impact from
different $y^{*}_{\mu}$ acceptances, the correction factor $k_{\mathrm{FB}}$
for the overall $R_{\mathrm{FB}}$ measurement is calculated to be $0.65\pm
0.02$.

The value of $R_{\mathrm{FB}}$ is then measured to be $R_{\mathrm{FB}}=0.78\pm
0.10$, as shown in Figure 6. The measured $R_{\mathrm{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,
with a better agreement as compared to the nCTEQ15 nPDF. 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_{\mathrm{FB}}$ is also measured differentially as
a function of $y^{*}_{{Z}}$, $p_{\mathrm{T}}^{{Z}}$ and $\phi^{*}$. The
corresponding $k_{\mathrm{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_{\mathrm{FB}}$ as a function of $y^{*}_{{Z}}$, $p_{\mathrm{T}}^{{Z}}$ and
$\phi^{*}$, respectively, together with the theory calculations. The
$R_{\mathrm{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_{\mathrm{FB}}$
measurements as a function of $p_{\mathrm{T}}^{{Z}}$ and $\phi^{*}$, a general
suppression is expected and observed. From theory a slightly stronger
suppression at lower $p_{\mathrm{T}}^{{Z}}$ or $\phi^{*}$ is expected, but
given the precision of the current measurement this tendency is not visible.
The $R_{\mathrm{FB}}$ as a function of $\phi^{*}$ shows a more stable trend
than as a function of $p_{\mathrm{T}}^{{Z}}$, which exhibits a larger
fluctuation The detailed numerical values and uncertainties of the
differential measurements of the $R_{\mathrm{FB}}$ together with the
corresponding $f_{\mathrm{FSR}}$ factors are given in Appendix B.

(a)(b)
(c)(d)

Figure 6:  The measured forward-backward ratio ($R_{\mathrm{FB}}$), (a) for the overall measurement, (b) as a function of $y^{*}_{{Z}}$, (c) as a function of $p_{\mathrm{T}}^{{Z}}$, and (d) as a function of $\phi^{*}$, together with the PowhegBox predictions using CTEQ6.1, EPPS16 and nCTEQ15 (n)PDF sets.  Table 2: The $y^{*}_{\mu}$ acceptance correction factors ($k_{\mathrm{FB}}$) for $R_{\mathrm{FB}}$ measured in intervals of $y^{*}_{{Z}}$, $p_{\mathrm{T}}^{{Z}}$, and $\phi^{*}$. $y^{*}_{{Z}}$ | $[2.5,3.0]$ | $[3.0,3.5]$ | $[3.5,4.0]$ |  |  |
---|---|---|---|---|---|---
$k_{\mathrm{FB}}(y^{*}_{{Z}})$ | $0.29\pm 0.01$ | $0.94\pm 0.04$ | $2.73\pm 0.13$ |  |  |
$p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-0.80005ptV]}$ | $[0,4.6]$ | $[4.6,10.5]$ | $[10.5,19]$ | $[19,34]$ | $[34,63]$ | $[63,270]$
$k_{\mathrm{FB}}(p_{\mathrm{T}}^{{Z}})$ | $0.72\pm 0.03$ | $0.65\pm 0.03$ | $0.63\pm 0.03$ | $0.59\pm 0.03$ | $0.58\pm 0.02$ | $0.46\pm 0.02$
$\phi^{*}$ | $[0,0.02]$ | $[0.02,0.05]$ | $[0.05,0.1]$ | $[0.1,0.2]$ | $[0.2,0.4]$ | $[0.4,4]$
$k_{\mathrm{FB}}(\phi^{*})$ | $0.70\pm 0.03$ | $0.67\pm 0.03$ | $0.65\pm 0.03$ | $0.62\pm 0.03$ | $0.59\pm 0.03$ | $0.49\pm 0.02$

For the total nuclear modification factor measurements, because the
$y^{*}_{{Z}}$ coverage between $pp$ and ${p}\mathrm{Pb}$ 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^{*}_{\mu}$ acceptance correction factors are determined to be
$k_{{p}\mathrm{Pb}}^{\rm fw.}=0.706\pm 0.002$ and $k_{{p}\mathrm{Pb}}^{\rm
bw.}=1.518\pm 0.003$ for forward and backward rapidities, respectively. The
cross-sections to be used in the $R_{{p}\mathrm{Pb}}$ calculations are
measured to be $23.7\pm 2.1\text{\,nb}$ and $13.1\pm 1.2\text{\,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\text{\,Te\kern-1.00006ptV}$ are interpolated in the
corresponding $y^{*}_{{Z}}$ common acceptance ranges as $(95.18\pm 0.55)\times
10^{-3}\text{\,nb}$ and $(79.10\pm 0.51)\times 10^{-3}\text{\,nb}$ for forward
and backward rapidities, respectively.

The measured overall nuclear modification factors are $R_{{p}\mathrm{Pb}}^{\rm
fw}=0.94\pm 0.07$ and $R_{{p}\mathrm{Pb}}^{\rm bw}=1.21\pm 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_{{p}\mathrm{Pb}}$ 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$\sigma$ 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_{{p}\mathrm{Pb}}$, the forward and
backward cross-sections and the corresponding $pp$ reference cross-sections
are determined in the same way as for the total $R_{{p}\mathrm{Pb}}$
measurement but in intervals of $y^{*}_{{Z}}$, $p_{\mathrm{T}}^{{Z}}$ and
$\phi^{*}$. The $y^{*}_{\mu}$ acceptance correction factors
$k_{{p}\mathrm{Pb}}$ 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_{{p}\mathrm{Pb}}$ values as a function of $y^{*}_{{Z}}$,
$p_{\mathrm{T}}^{{Z}}$ and $\phi^{*}$ 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_{{p}\mathrm{Pb}}$ as a function
of $y^{*}_{{Z}}$, $p_{\mathrm{T}}^{{Z}}$, and $\phi^{*}$, together with the
corresponding $f_{\mathrm{FSR}}$ factors are given in Appendix B.

Figure 7:  Measurements of the nuclear modification factor ($R_{{p}\mathrm{Pb}}$) compared to the PowhegBox predictions using the EPPS16 and nCTEQ15 (n)PDF sets. Table 3: The $y^{*}_{\mu}$ acceptance correction factors ($k_{{p}\mathrm{Pb}}$) for $R_{{p}\mathrm{Pb}}$ measured in intervals of $y^{*}_{{Z}}$, $p_{\mathrm{T}}^{{Z}}$, and $\phi^{*}$ for forward and backward collisions. | $y^{*}_{{Z}}$ | $[2.0,2.5]$ | $[2.5,3.0]$ | $[3.0,3.5]$ | $[3.5,4.0]$ |  |
---|---|---|---|---|---|---|---
$k_{{p}\mathrm{Pb}}$ ($y^{*}_{{Z}}$) | forward | $0.3713$ | $0.7468$ | $1.2482$ | $2.1875$ |  |
|  | $\pm 0.0005$ | $\pm 0.0004$ | $\pm 0.0010$ | $\pm 0.0004$ |  |
| $y^{*}_{{Z}}$ | $[-4.0,-3.5]$ | $[-3.5,-3.0]$ | $[-3.0,-2.5]$ |  |  |
$k_{{p}\mathrm{Pb}}$ ($y^{*}_{{Z}}$) | backward | $0.7935$ | $1.3177$ | $2.5525$ |  |  |
|  | $\pm 0.0001$ | $\pm 0.0001$ | $\pm 0.0009$ |  |  |
|  $p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-0.80005ptV]}$ | $[0,4.6]$ | $[4.6,10.5]$ | $[10.5,19]$ | $[19,34]$ | $[34,63]$ | $[63,270]$
$k_{{p}\mathrm{Pb}}$ ($p_{\mathrm{T}}^{{Z}}$) | forward | $0.7821$ | $0.7212$ | $0.6903$ | $0.6613$ | $0.6136$ | $0.5455$
|  | $\pm 0.0056$ | $\pm 0.0049$ | $\pm 0.0045$ | $\pm 0.0040$ | $\pm 0.0035$ | $\pm 0.0027$
| backward | $1.6073$ | $1.7735$ | $1.8341$ | $1.9387$ | $2.0742$ | $2.8031$
|  | $\pm 0.0034$ | $\pm 0.0039$ | $\pm 0.0040$ | $\pm 0.0045$ | $\pm 0.0047$ | $\pm 0.0069$
| $\phi^{*}$ | $[0,0.02]$ | $[0.02,0.05]$ | $[0.05,0.1]$ | $[0.1,0.2]$ | $[0.2,0.4]$ | $[0.4,4]$
$k_{{p}\mathrm{Pb}}$ ($\phi^{*}$) | forward | $0.7702$ | $0.7322$ | $0.7080$ | $0.6804$ | $0.6461$ | $0.5756$
|  | $\pm 0.0054$ | $\pm 0.0049$ | $\pm 0.0047$ | $\pm 0.0043$ | $\pm 0.0040$ | $\pm 0.0032$
| backward | $1.6400$ | $1.7113$ | $1.8105$ | $1.8645$ | $1.9741$ | $2.4238$
|  | $\pm 0.0036$ | $\pm 0.0036$ | $\pm 0.0041$ | $\pm 0.043$ | $\pm 0.0047$ | $\pm 0.0065$

(a)(b)
(c)(d)
(e)(f)

Figure 8:  Nuclear modification factors ($R_{{p}\mathrm{Pb}}$) as a function
of (top row) $y^{*}_{{Z}}$, (middle row) $p_{\mathrm{T}}^{{Z}}$ and (bottom
row) $\phi^{*}$, together with the PowhegBox prediction using EPPS16 and
nCTEQ15 nPDF sets, where the left column is for forward collisions and the
right column is for backward collisions.

##  7 Conclusion

This article presents a comprehensive analysis of
${Z}\rightarrow{\mu^{+}\mu^{-}}$ production using ${p}\mathrm{Pb}$ collisions
at 8.16 TeV in the forward region recorded with the LHCb detector. The
$Z$-boson production fiducial cross-section, $R_{\mathrm{FB}}$ and
$R_{{p}\mathrm{Pb}}$ are measured inclusively and differentially as a function
of $y^{*}_{{Z}}$, $p_{\mathrm{T}}^{{Z}}$, and $\phi^{*}$. The results are
compatible with theoretical predictions from the EPPS16 and nCTEQ15 nPDFs. The
measurements are also compatible with previous results from various
experiments at 5.02 TeV. The precision of the measurements is significantly
improved compared to the previous LHCb results using ${p}\mathrm{Pb}$
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}$.

## Appendices

##  Appendix A Statistical and systematic correlation matrices

Table 4: Correlation matrix of statistical uncertainty in $y^{*}_{{Z}}$ intervals for forward collisions. $y^{*}_{{Z}}$ | Index | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---|---
$[1.5,2.0]$ | 1 | 1.00 |  |  |  |
$[2.0,2.5]$ | 2 | 0.01 | 1.00 |  |  |
$[2.5,3.0]$ | 3 | 0.00 | 0.01 | 1.00 |  |
$[3.0,3.5]$ | 4 | 0.00 | 0.00 | 0.01 | 1.00 |
$[3.5,4.0]$ | 5 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00
Table 5: Correlation matrix of statistical uncertainty in $y^{*}_{{Z}}$ intervals for backward collisions. $y^{*}_{{Z}}$ | Index | 1 | 2 | 3
---|---|---|---|---
$[-3.0,-2.5]$ | 1 | 1.00 |  |
$[-3.5,-3.0]$ | 2 | 0.01 | 1.00 |
$[-4.0,-3.5]$ | 3 | 0.00 | 0.01 | 1.00
Table 6: Correlation matrix of statistical uncertainty in $p_{\mathrm{T}}^{{Z}}$ intervals for forward collisions. $p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-1.00006ptV]}$ | Index | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---|---
$[0.0,4.6]$ | 1 | 1.00 |  |  |  |  |
$[4.6,10.5]$ | 2 | 0.93 | 1.00 |  |  |  |
$[10.5,19.0]$ | 3 | 0.00 | 0.05 | 1.00 |  |  |
$[19.0,34.0]$ | 4 | 0.00 | 0.00 | 0.04 | 1.00 |  |
$[34.0,63.0]$ | 5 | 0.00 | 0.00 | 0.00 | 0.03 | 1.00 |
$[63.0,270.0]$ | 6 | 0.00 | 0.00 | 0.00 | 0.00 | 0.02 | 1.00
Table 7: Correlation matrix of statistical uncertainty in $p_{\mathrm{T}}^{{Z}}$ intervals for backward collisions. $p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-1.00006ptV]}$ | Index | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---|---
$[0.0,4.6]$ | 1 | 1.00 |  |  |  |  |
$[4.6,10.5]$ | 2 | 0.08 | 1.00 |  |  |  |
$[10.5,19.0]$ | 3 | 0.00 | 0.04 | 1.00 |  |  |
$[19.0,34.0]$ | 4 | 0.00 | 0.00 | 0.03 | 1.00 |  |
$[34.0,63.0]$ | 5 | 0.00 | 0.00 | 0.00 | 0.02 | 1.00 |
$[63.0,270.0]$ | 6 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 1.00
Table 8: Correlation matrix of statistical uncertainty in $\phi^{*}$ intervals for forward collisions. $\phi^{*}$ | Index | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---|---
$[0.00,0.02]$ | 1 | 1.00 |  |  |  |  |
$[0.02,0.05]$ | 2 | 0.01 | 1.00 |  |  |  |
$[0.05,0.10]$ | 3 | 0.00 | 0.00 | 1.00 |  |  |
$[0.10,0.20]$ | 4 | 0.00 | 0.00 | 0.00 | 1.00 |  |
$[0.20,0.40]$ | 5 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 |
$[0.40,4.00]$ | 6 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00
Table 9: Correlation matrix of statistical uncertainty in $\phi^{*}$ intervals for backward collisions. $\phi^{*}$ | Index | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---|---
$[0.00,0.02]$ | 1 | 1.00 |  |  |  |  |
$[0.02,0.05]$ | 2 | 0.00 | 1.00 |  |  |  |
$[0.05,0.10]$ | 3 | 0.00 | 0.00 | 1.00 |  |  |
$[0.10,0.20]$ | 4 | 0.00 | 0.00 | 0.00 | 1.00 |  |
$[0.20,0.40]$ | 5 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 |
$[0.40,4.00]$ | 6 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00
Table 10: Correlation matrix of statistical uncertainty in $p_{\mathrm{T}}^{{Z}}$ intervals for forward collisions in the low $p_{\mathrm{T}}^{{Z}}$ range. $p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-1.00006ptV]}$ | Index | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---|---
$[0.0,2.2]$ | 1 | 1.00 |  |  |  |  |  |  |
$[2.2,4.6]$ | 2 | 0.24 | 1.00 |  |  |  |  |  |
$[4.6,7.2]$ | 3 | 0.04 | 0.26 | 1.00 |  |  |  |  |
$[7.2,8.7]$ | 4 | 0.01 | 0.04 | 0.25 | 1.00 |  |  |  |
$[8.7,10.5]$ | 5 | 0.01 | 0.01 | 0.04 | 0.24 | 1.00 |  |  |
$[10.5,12.8]$ | 6 | 0.00 | 0.01 | 0.01 | 0.03 | 0.21 | 1.00 |  |
$[12.8,15.4]$ | 7 | 0.00 | 0.00 | 0.01 | 0.01 | 0.03 | 0.20 | 1.00 |
$[15.4,19.0]$ | 8 | 0.00 | 0.00 | 0.00 | 0.01 | 0.01 | 0.02 | 0.18 | 1.00
Table 11: Correlation matrix of statistical uncertainty in $p_{\mathrm{T}}^{{Z}}$ intervals for backward collisions in the low $p_{\mathrm{T}}^{{Z}}$ range. $p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-1.00006ptV]}$ | Index | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---|---
$[0.0,2.2]$ | 1 | 1.00 |  |  |  |  |  |  |
$[2.2,4.6]$ | 2 | 0.20 | 1.00 |  |  |  |  |  |
$[4.6,7.2]$ | 3 | 0.03 | 0.22 | 1.00 |  |  |  |  |
$[7.2,8.7]$ | 4 | 0.01 | 0.03 | 0.21 | 1.00 |  |  |  |
$[8.7,10.5]$ | 5 | 0.01 | 0.01 | 0.02 | 0.20 | 1.00 |  |  |
$[10.5,12.8]$ | 6 | 0.00 | 0.00 | 0.01 | 0.02 | 0.18 | 1.00 |  |
$[12.8,15.4]$ | 7 | 0.00 | 0.00 | 0.00 | 0.01 | 0.02 | 0.17 | 1.00 |
$[15.4,19.0]$ | 8 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 0.02 | 0.14 | 1.00
Table 12: Correlation matrix of efficiency uncertainty in $y^{*}_{{Z}}$ intervals for forward collisions. $y^{*}_{{Z}}$ | Index | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---|---
$[2.0,2.5]$ | 1 | 1.00 |  |  |  |
$[2.5,3.0]$ | 2 | 0.71 | 1.00 |  |  |
$[3.0,3.5]$ | 3 | 0.37 | 0.72 | 1.00 |  |
$[3.5,4.0]$ | 4 | 0.03 | 0.24 | 0.71 | 1.00 |
$[4.0,4.5]$ | 5 | 0.00 | 0.05 | 0.45 | 0.79 | 1.00
Table 13: Correlation matrix of efficiency uncertainty in $y^{*}_{{Z}}$ intervals for backward collisions. $y^{*}_{{Z}}$ | Index | 1 | 2 | 3
---|---|---|---|---
$[-3.0,-2.5]$ | 1 | 1.00 |  |
$[-3.5,-3.0]$ | 2 | 0.64 | 1.00 |
$[-4.0,-3.5]$ | 3 | 0.48 | 0.75 | 1.00
Table 14: Correlation matrix of efficiency uncertainty in $p_{\mathrm{T}}^{{Z}}$ intervals for forward collisions. $p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-1.00006ptV]}$ | Index | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---|---
$[0.0,4.6]$ | 1 | 1.00 |  |  |  |  |
$[4.6,10.5]$ | 2 | 0.88 | 1.00 |  |  |  |
$[10.5,19.0]$ | 3 | 0.87 | 0.89 | 1.00 |  |  |
$[19.0,34.0]$ | 4 | 0.92 | 0.89 | 0.91 | 1.00 |  |
$[34.0,63.0]$ | 5 | 0.85 | 0.80 | 0.79 | 0.87 | 1.00 |
$[63.0,270.0]$ | 6 | 0.75 | 0.82 | 0.82 | 0.82 | 0.74 | 1.00
Table 15: Correlation matrix of efficiency uncertainty in $p_{\mathrm{T}}^{{Z}}$ intervals for backward collisions. $p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-1.00006ptV]}$ | Index | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---|---
$[0.0,4.6]$ | 1 | 1.00 |  |  |  |  |
$[4.6,10.5]$ | 2 | 0.85 | 1.00 |  |  |  |
$[10.5,19.0]$ | 3 | 0.91 | 0.91 | 1.00 |  |  |
$[19.0,34.0]$ | 4 | 0.83 | 0.94 | 0.86 | 1.00 |  |
$[34.0,63.0]$ | 5 | 0.68 | 0.84 | 0.65 | 0.80 | 1.00 |
$[63.0,270.0]$ | 6 | 0.50 | 0.71 | 0.65 | 0.70 | 0.77 | 1.00
Table 16: Correlation matrix of efficiency uncertainty in $\phi^{*}$ intervals for forward collisions. $\phi^{*}$ | Index | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---|---
$[0.00,0.02]$ | 1 | 1.00 |  |  |  |  |
$[0.02,0.05]$ | 2 | 0.83 | 1.00 |  |  |  |
$[0.05,0.10]$ | 3 | 0.92 | 0.90 | 1.00 |  |  |
$[0.10,0.20]$ | 4 | 0.83 | 0.85 | 0.80 | 1.00 |  |
$[0.20,0.40]$ | 5 | 0.85 | 0.71 | 0.77 | 0.73 | 1.00 |
$[0.40,4.00]$ | 6 | 0.85 | 0.89 | 0.85 | 0.81 | 0.82 | 1.00
Table 17: Correlation matrix of efficiency uncertainty in $\phi^{*}$ intervals for backward collisions. $\phi^{*}$ | Index | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---|---
$[0.00,0.02]$ | 1 | 1.00 |  |  |  |  |
$[0.02,0.05]$ | 2 | 0.90 | 1.00 |  |  |  |
$[0.05,0.10]$ | 3 | 0.88 | 0.80 | 1.00 |  |  |
$[0.10,0.20]$ | 4 | 0.80 | 0.79 | 0.77 | 1.00 |  |
$[0.20,0.40]$ | 5 | 0.76 | 0.70 | 0.79 | 0.58 | 1.00 |
$[0.40,4.00]$ | 6 | 0.53 | 0.60 | 0.55 | 0.75 | 0.54 | 1.00
Table 18: Correlation matrix of efficiency uncertainty in
$p_{\mathrm{T}}^{{Z}}$ intervals for forward collisions in the low
$p_{\mathrm{T}}^{{Z}}$ range.

$p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-0.80005ptV]}$ | Index | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---|---
$[0.0,2.2]$ | 1 | 1.00 |  |  |  |  |  |  |
$[2.2,4.6]$ | 2 | 0.57 | 1.00 |  |  |  |  |  |
$[4.6,7.2]$ | 3 | 0.71 | 0.83 | 1.00 |  |  |  |  |
$[7.2,8.7]$ | 4 | 0.57 | 0.67 | 0.68 | 1.00 |  |  |  |
$[8.7,10.5]$ | 5 | 0.61 | 0.71 | 0.80 | 0.68 | 1.00 |  |  |
$[10.5,12.8]$ | 6 | 0.72 | 0.73 | 0.83 | 0.64 | 0.85 | 1.00 |  |
$[12.8,15.4]$ | 7 | 0.66 | 0.77 | 0.84 | 0.62 | 0.71 | 0.85 | 1.00 |
$[15.4,19.0]$ | 8 | 0.58 | 0.78 | 0.78 | 0.46 | 0.76 | 0.77 | 0.70 | 1.00

Table 19: Correlation matrix of efficiency uncertainty in
$p_{\mathrm{T}}^{{Z}}$ intervals for backward collisions in the low
$p_{\mathrm{T}}^{{Z}}$ range.

$p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-0.80005ptV]}$ | Index | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---|---
$[0.0,2.2]$ | 1 | 1.00 |  |  |  |  |  |  |
$[2.2,4.6]$ | 2 | 0.82 | 1.00 |  |  |  |  |  |
$[4.6,7.2]$ | 3 | 0.73 | 0.76 | 1.00 |  |  |  |  |
$[7.2,8.7]$ | 4 | 0.78 | 0.66 | 0.70 | 1.00 |  |  |  |
$[8.7,10.5]$ | 5 | 0.48 | 0.74 | 0.58 | 0.48 | 1.00 |  |  |
$[10.5,12.8]$ | 6 | 0.86 | 0.84 | 0.69 | 0.61 | 0.58 | 1.00 |  |
$[12.8,15.4]$ | 7 | 0.72 | 0.72 | 0.67 | 0.88 | 0.55 | 0.52 | 1.00 |
$[15.4,19.0]$ | 8 | 0.79 | 0.78 | 0.81 | 0.81 | 0.55 | 0.74 | 0.71 | 1.00

##  Appendix B Differential cross-section, forward-backward ratio, and nuclear
modification factors

Table 20:  Differential cross-section of ${Z}\rightarrow{\mu^{+}\mu^{-}}$ 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. | $y^{*}_{{Z}}$ | $d\sigma/dy^{*}_{{Z}}\,[\text{\,nb}]$ | $f_{\mathrm{FSR}}$
---|---|---|---
Forward | $[1.5,2.0]$ | $5.60\pm 1.08\pm 0.20\pm 0.15$ | $1.017\pm 0.002$
| $[2.0,2.5]$ | $16.19\pm 1.80\pm 0.53\pm 0.43$ | $1.020\pm 0.001$
| $[2.5,3.0]$ | $20.36\pm 2.02\pm 0.72\pm 0.53$ | $1.023\pm 0.011$
| $[3.0,3.5]$ | $10.05\pm 1.42\pm 0.36\pm 0.26$ | $1.030\pm 0.006$
| $[3.5,4.0]$ | $1.59\pm 0.56\pm 0.06\pm 0.04$ | $1.014\pm 0.002$
Backward | $[-4.0,-3.5]$ | $4.87\pm 0.89\pm 0.17\pm 0.18$ | $1.028\pm 0.002$
| $[-3.5,-3.0]$ | $14.02\pm 1.49\pm 0.49\pm 0.53$ | $1.018\pm 0.007$
| $[-3.0,-2.5]$ | $7.83\pm 1.13\pm 0.30\pm 0.30$ | $1.020\pm 0.002$
Table 21:  Differential cross-section of ${Z}\rightarrow{\mu^{+}\mu^{-}}$ in intervals of $p_{\mathrm{T}}^{{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. |  $p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-1.00006ptV]}$ | $d\sigma/dp_{\mathrm{T}}^{{Z}}\,[\text{\,nb}/\text{\,Ge\kern-1.00006ptV}]$ | $f_{\mathrm{FSR}}$
---|---|---|---
Forward | $[0.0,4.6]$ | $1.441\pm 0.209\pm 0.046\pm 0.038$ | $1.068\pm 0.009$
| $[4.6,10.5]$ | $1.086\pm 0.156\pm 0.035\pm 0.029$ | $1.033\pm 0.007$
| $[10.5,19.0]$ | $0.612\pm 0.092\pm 0.020\pm 0.016$ | $0.994\pm 0.008$
| $[19.0,34.0]$ | $0.272\pm 0.046\pm 0.010\pm 0.007$ | $0.988\pm 0.007$
| $[34.0,63.0]$ | $0.124\pm 0.023\pm 0.005\pm 0.003$ | $1.034\pm 0.010$
| $[63.0,270.0]$ | $(5.9\pm 1.9\pm 0.4\pm 0.2)\times 10^{-3}$ | $1.043\pm 0.035$
Backward | $[0.0,4.6]$ | $1.097\pm 0.162\pm 0.042\pm 0.027$ | $1.082\pm 0.021$
| $[4.6,10.5]$ | $0.434\pm 0.096\pm 0.016\pm 0.011$ | $1.025\pm 0.013$
| $[10.5,19.0]$ | $0.218\pm 0.051\pm 0.009\pm 0.005$ | $0.976\pm 0.019$
| $[19.0,34.0]$ | $0.201\pm 0.040\pm 0.009\pm 0.005$ | $0.981\pm 0.023$
| $[34.0,63.0]$ | $0.031\pm 0.011\pm 0.002\pm 0.001$ | $1.028\pm 0.030$
| $[63.0,270.0]$ | $(1.3\pm 0.8\pm 0.1\pm 0.1)\times 10^{-3}$ | $1.064\pm 0.011$
Table 22:  Differential cross-section of ${Z}\rightarrow{\mu^{+}\mu^{-}}$ in intervals of $p_{\mathrm{T}}^{{Z}}$ for forward and backward collisions in the low $p_{\mathrm{T}}^{{Z}}$ 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. |  $p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-1.00006ptV]}$ | $d\sigma/dp_{\mathrm{T}}^{{Z}}\,[\text{\,nb}/\text{\,Ge\kern-1.00006ptV}]$ | $f_{\mathrm{FSR}}$
---|---|---|---
Forward | $[0.0,2.2]$ | $1.463\pm 0.352\pm 0.059\pm 0.038$ | $1.079\pm 0.028$
| $[2.2,4.6]$ | $1.450\pm 0.347\pm 0.045\pm 0.038$ | $1.064\pm 0.003$
| $[4.6,7.2]$ | $1.414\pm 0.310\pm 0.044\pm 0.037$ | $1.048\pm 0.004$
| $[7.2,8.7]$ | $0.742\pm 0.344\pm 0.027\pm 0.019$ | $1.035\pm 0.020$
| $[8.7,10.5]$ | $0.906\pm 0.261\pm 0.044\pm 0.024$ | $1.005\pm 0.038$
| $[10.5,12.8]$ | $0.901\pm 0.227\pm 0.037\pm 0.024$ | $1.003\pm 0.026$
| $[12.8,15.4]$ | $0.534\pm 0.167\pm 0.020\pm 0.014$ | $1.000\pm 0.020$
| $[15.4,19.0]$ | $0.475\pm 0.130\pm 0.022\pm 0.012$ | $0.981\pm 0.031$
Backward | $[0.0,2.2]$ | $0.813\pm 0.223\pm 0.029\pm 0.020$ | $1.107\pm 0.014$
| $[2.2,4.6]$ | $1.352\pm 0.259\pm 0.062\pm 0.033$ | $1.072\pm 0.033$
| $[4.6,7.2]$ | $0.596\pm 0.166\pm 0.024\pm 0.015$ | $1.037\pm 0.021$
| $[7.2,8.7]$ | $0.438\pm 0.187\pm 0.024\pm 0.011$ | $1.025\pm 0.046$
| $[8.7,10.5]$ | $0.220\pm 0.119\pm 0.011\pm 0.005$ | $0.999\pm 0.036$
| $[10.5,12.8]$ | $0.232\pm 0.104\pm 0.008\pm 0.006$ | $0.991\pm 0.008$
| $[12.8,15.4]$ | $0.168\pm 0.083\pm 0.008\pm 0.004$ | $0.965\pm 0.035$
| $[15.4,19.0]$ | $0.250\pm 0.083\pm 0.012\pm 0.006$ | $0.970\pm 0.034$
Table 23:  Differential cross-section of ${Z}\rightarrow{\mu^{+}\mu^{-}}$ in intervals of $\phi^{*}$ 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. | $\phi^{*}$ | $d\sigma/d\phi^{*}\,[\text{\,nb}]$ | $f_{\mathrm{FSR}}$
---|---|---|---
Forward | $[0.0,0.0]$ | $261.288\pm 36.234\pm 10.959\pm 6.865$ | $1.031\pm 0.028$
| $[0.0,0.1]$ | $179.805\pm 24.698\pm\phantom{1}6.125\pm 4.724$ | $1.044\pm 0.012$
| $[0.1,0.1]$ | $104.654\pm 14.375\pm\phantom{1}3.886\pm 2.750$ | $1.010\pm 0.018$
| $[0.1,0.2]$ | $46.938\pm\phantom{1}6.847\pm\phantom{1}1.765\pm 1.233$ | $1.018\pm 0.017$
| $[0.2,0.4]$ | $17.902\pm\phantom{1}2.984\pm\phantom{1}0.656\pm 0.470$ | $1.004\pm 0.002$
| $[0.4,4.0]$ | $0.786\pm\phantom{1}0.151\pm\phantom{1}0.042\pm 0.021$ | $1.024\pm 0.024$
Backward | $[0.0,0.0]$ | $173.372\pm 26.752\pm\phantom{1}6.029\pm 4.292$ | $1.053\pm 0.006$
| $[0.0,0.1]$ | $99.962\pm 16.216\pm\phantom{1}3.547\pm 2.475$ | $1.007\pm 0.009$
| $[0.1,0.1]$ | $46.233\pm\phantom{1}8.585\pm\phantom{1}1.672\pm 1.145$ | $1.015\pm 0.010$
| $[0.1,0.2]$ | $24.080\pm\phantom{1}4.396\pm\phantom{1}0.863\pm 0.596$ | $1.014\pm 0.001$
| $[0.2,0.4]$ | $8.552\pm\phantom{1}1.866\pm\phantom{1}0.355\pm 0.212$ | $1.015\pm 0.015$
| $[0.4,4.0]$ | $0.165\pm\phantom{1}0.062\pm\phantom{1}0.009\pm 0.004$ | $1.023\pm 0.018$
Table 24: The forward-backward ratio ($R_{\mathrm{FB}}$) in intervals of $|y^{*}_{{Z}}|$, together with the FSR correction. For the $R_{\mathrm{FB}}$ results, the first uncertainty is statistical, the second is systematic and the third is from integrated luminosity. $|y^{*}_{{Z}}|$ | $R_{\mathrm{FB}}$ |  $f_{\mathrm{FSR}}$ (forward) |  $f_{\mathrm{FSR}}$ (backward)
---|---|---|---
$[2.5,3.0]$ | $0.75\pm 0.13\pm 0.27\pm 0.26$ | $1.023\pm 0.011$ | $1.020\pm 0.002$
$[3.0,3.5]$ | $0.67\pm 0.12\pm 0.24\pm 0.24$ | $1.030\pm 0.006$ | $1.018\pm 0.007$
$[3.5,4.0]$ | $0.89\pm 0.36\pm 0.71\pm 0.71$ | $1.014\pm 0.002$ | $1.028\pm 0.002$
Table 25: The forward-backward ratio ($R_{\mathrm{FB}}$) in intervals of $p_{\mathrm{T}}^{{Z}}$, together with the FSR correction. For the $R_{\mathrm{FB}}$ results, the first uncertainty is statistical, the second is systematic and the third is from integrated luminosity. $p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-1.00006ptV]}$ | $R_{\mathrm{FB}}$ |  $f_{\mathrm{FSR}}$ (forward) |  $f_{\mathrm{FSR}}$ (backward)
---|---|---|---
$[0.0,4.6]$ | $0.67\pm 0.15\pm 0.04\pm 0.02$ | $1.068\pm 0.009$ | $1.082\pm 0.021$
$[4.6,10.5]$ | $0.92\pm 0.27\pm 0.06\pm 0.03$ | $1.033\pm 0.007$ | $1.025\pm 0.013$
$[10.5,19.0]$ | $1.25\pm 0.37\pm 0.08\pm 0.05$ | $0.994\pm 0.008$ | $0.976\pm 0.019$
$[19.0,34.0]$ | $0.44\pm 0.14\pm 0.03\pm 0.02$ | $0.988\pm 0.007$ | $0.981\pm 0.023$
$[34.0,63.0]$ | $1.03\pm 0.46\pm 0.08\pm 0.04$ | $1.034\pm 0.010$ | $1.028\pm 0.030$
$[63.0,270.0]$ | $0.74\pm 0.59\pm 0.09\pm 0.03$ | $1.043\pm 0.035$ | $1.064\pm 0.011$
Table 26: The forward-backward ratio ($R_{\mathrm{FB}}$) in intervals of $\phi^{*}$, together with the FSR correction. For the $R_{\mathrm{FB}}$ results, the first uncertainty is statistical, the second is systematic and the third is from integrated luminosity. $\phi^{*}$ | $R_{\mathrm{FB}}$ |  $f_{\mathrm{FSR}}$ (forward) |  $f_{\mathrm{FSR}}$ (backward)
---|---|---|---
$[0.00,0.02]$ | $0.71\pm 0.16\pm 0.05\pm 0.03$ | $1.031\pm 0.028$ | $1.053\pm 0.006$
$[0.02,0.05]$ | $0.82\pm 0.19\pm 0.05\pm 0.03$ | $1.044\pm 0.012$ | $1.007\pm 0.009$
$[0.05,0.10]$ | $0.86\pm 0.22\pm 0.06\pm 0.03$ | $1.010\pm 0.018$ | $1.015\pm 0.010$
$[0.10,0.20]$ | $0.75\pm 0.19\pm 0.05\pm 0.03$ | $1.018\pm 0.017$ | $1.014\pm 0.001$
$[0.20,0.40]$ | $0.65\pm 0.21\pm 0.05\pm 0.02$ | $1.004\pm 0.002$ | $1.015\pm 0.015$
$[0.40,4.00]$ | $0.94\pm 0.46\pm 0.08\pm 0.03$ | $1.024\pm 0.024$ | $1.023\pm 0.018$
Table 27: The nuclear modification factors ($R_{{p}\mathrm{Pb}}$) in intervals of $y^{*}_{{Z}}$ for forward and backward collisions, together with the FSR correction. For the $R_{{p}\mathrm{Pb}}$ results, the first uncertainty is statistical, the second is systematic and the third is from integrated luminosity. | $y^{*}_{{Z}}$ | $R_{{p}\mathrm{Pb}}$ | $f_{\mathrm{FSR}}$
---|---|---|---
Forward | $[2.0,2.5]$ | $0.87\pm 0.10\pm 0.03\pm 0.02$ | $1.020\pm 0.001$
| $[2.5,3.0]$ | $0.98\pm 0.10\pm 0.04\pm 0.03$ | $1.023\pm 0.011$
| $[3.0,3.5]$ | $0.89\pm 0.13\pm 0.03\pm 0.02$ | $1.030\pm 0.006$
| $[3.5,4.0]$ | $0.83\pm 0.29\pm 0.03\pm 0.02$ | $1.014\pm 0.002$
Backward | $[-4.0,-3.5]$ | $0.92\pm 0.17\pm 0.03\pm 0.03$ | $1.028\pm 0.002$
| $[-3.5,-3.0]$ | $1.31\pm 0.14\pm 0.05\pm 0.05$ | $1.018\pm 0.007$
| $[-3.0,-2.5]$ | $1.29\pm 0.19\pm 0.05\pm 0.05$ | $1.020\pm 0.002$
Table 28: The nuclear modification factors ($R_{{p}\mathrm{Pb}}$) in intervals of $p_{\mathrm{T}}^{{Z}}$ for forward and backward collisions, together with the FSR correction. For the $R_{{p}\mathrm{Pb}}$ results, the first uncertainty is statistical, the second is systematic and the third is from integrated luminosity. |  $p_{\mathrm{T}}^{{Z}}$ $\mathrm{[Ge\kern-1.00006ptV]}$ | $R_{{p}\mathrm{Pb}}$ | $f_{\mathrm{FSR}}$
---|---|---|---
Forward | $[0.0,4.6]$ | $1.04\pm 0.15\pm 0.04\pm 0.03$ | $1.068\pm 0.009$
| $[4.6,10.5]$ | $0.77\pm 0.11\pm 0.03\pm 0.02$ | $1.033\pm 0.007$
| $[10.5,19.0]$ | $0.83\pm 0.13\pm 0.03\pm 0.02$ | $0.994\pm 0.008$
| $[19.0,34.0]$ | $0.89\pm 0.15\pm 0.03\pm 0.02$ | $0.988\pm 0.007$
| $[34.0,63.0]$ | $1.34\pm 0.25\pm 0.06\pm 0.04$ | $1.034\pm 0.010$
| $[63.0,270.0]$ | $1.41\pm 0.44\pm 0.10\pm 0.04$ | $1.043\pm 0.035$
Backward | $[0.0,4.6]$ | $1.63\pm 0.24\pm 0.07\pm 0.04$ | $1.082\pm 0.021$
| $[4.6,10.5]$ | $0.76\pm 0.17\pm 0.03\pm 0.02$ | $1.025\pm 0.013$
| $[10.5,19.0]$ | $0.79\pm 0.18\pm 0.03\pm 0.02$ | $0.976\pm 0.019$
| $[19.0,34.0]$ | $1.94\pm 0.39\pm 0.09\pm 0.05$ | $0.981\pm 0.023$
| $[34.0,63.0]$ | $1.13\pm 0.40\pm 0.06\pm 0.03$ | $1.028\pm 0.030$
| $[63.0,270.0]$ | $1.63\pm 0.98\pm 0.11\pm 0.04$ | $1.064\pm 0.011$
Table 29: The nuclear modification factors ($R_{{p}\mathrm{Pb}}$) in intervals of $\phi^{*}$ for forward and backward collisions, together with the FSR correction. For the $R_{{p}\mathrm{Pb}}$ results, the first uncertainty is statistical, the second is systematic and the third is from integrated luminosity. | $\phi^{*}$ | $R_{{p}\mathrm{Pb}}$ | $f_{\mathrm{FSR}}$
---|---|---|---
Forward | $[0.00,0.02]$ | $0.94\pm 0.13\pm 0.04\pm 0.02$ | $1.031\pm 0.028$
| $[0.02,0.05]$ | $0.82\pm 0.11\pm 0.03\pm 0.02$ | $1.044\pm 0.012$
| $[0.05,0.10]$ | $0.83\pm 0.11\pm 0.03\pm 0.02$ | $1.010\pm 0.018$
| $[0.10,0.20]$ | $0.87\pm 0.13\pm 0.04\pm 0.02$ | $1.018\pm 0.017$
| $[0.20,0.40]$ | $1.07\pm 0.18\pm 0.05\pm 0.03$ | $1.004\pm 0.002$
| $[0.40,4.00]$ | $1.52\pm 0.29\pm 0.09\pm 0.04$ | $1.024\pm 0.024$
Backward | $[0.00,0.02]$ | $1.32\pm 0.20\pm 0.05\pm 0.03$ | $1.053\pm 0.006$
| $[0.02,0.05]$ | $1.07\pm 0.17\pm 0.04\pm 0.03$ | $1.007\pm 0.009$
| $[0.05,0.10]$ | $0.94\pm 0.17\pm 0.04\pm 0.02$ | $1.015\pm 0.010$
| $[0.10,0.20]$ | $1.22\pm 0.22\pm 0.05\pm 0.03$ | $1.014\pm 0.001$
| $[0.20,0.40]$ | $1.57\pm 0.34\pm 0.07\pm 0.04$ | $1.015\pm 0.015$
| $[0.40,4.00]$ | $1.34\pm 0.51\pm 0.08\pm 0.03$ | $1.023\pm 0.018$

## References

  * [1] S. D. Drell and T.-M. Yan, _Massive lepton pair production in hadron-hadron collisions at high-energies_ , Phys. Rev. Lett. 25 (1970) 316, Erratum ibid. 25 (1970) 902
  * [2] J. C. Collins, D. E. Soper, and G. F. Sterman, _Factorization of hard processes in QCD_ , Adv. Ser. Direct. High Energy Phys. 5 (1989) 1, arXiv:hep-ph/0409313
  * [3] G. A. Ladinsky and C. P. Yuan, _The Nonperturbative regime in QCD resummation for gauge boson production at hadron colliders_ , Phys. Rev. D50 (1994) R4239, arXiv:hep-ph/9311341
  * [4] C. Balazs, J.-w. Qiu, and C. P. Yuan, _Effects of QCD resummation on distributions of leptons from the decay of electroweak vector bosons_ , Phys. Lett. B355 (1995) 548, arXiv:hep-ph/9505203
  * [5] C. Balazs and C. P. Yuan, _Testing multiple gluon dynamics at the Tevatron_ , Phys. Rev. Lett. 79 (1997) 2398, arXiv:hep-ph/9703405
  * [6] C. Balazs and C. P. Yuan, _Soft gluon effects on lepton pairs at hadron colliders_ , Phys. Rev. D56 (1997) 5558, arXiv:hep-ph/9704258
  * [7] F. Landry, R. Brock, P. M. Nadolsky, and C. P. Yuan, _Tevatron Run-1 $Z$ boson data and Collins-Soper-Sterman resummation formalism_, Phys. Rev. D67 (2003) 073016, arXiv:hep-ph/0212159
  * [8] A. V. Konychev and P. M. Nadolsky, _Universality of the Collins-Soper-Sterman nonperturbative function in gauge boson production_ , Phys. Lett. B633 (2006) 710, arXiv:hep-ph/0506225
  * [9] S. Camarda et al., _DYTurbo: Fast predictions for Drell-Yan processes_ , Eur. Phys. J. C80 (2020) 251, arXiv:1910.07049, [Erratum: Eur. Phys. J. C80, 440 (2020)]
  * [10] S. Camarda, L. Cieri, and G. Ferrera, _Drell–Yan lepton-pair production: qT resummation at N3LL accuracy and fiducial cross sections at N3LO_ , Phys. Rev. D104 (2021) L111503, arXiv:2103.04974
  * [11] R. Hamberg, W. L. van Neerven, and T. Matsuura, _A complete calculation of the order $\alpha_{s}^{2}$ correction to the Drell-Yan $K$-factor_, Nucl. Phys. B359 (1991) 343, Erratum ibid. 644 (2002) 403
  * [12] S. Catani et al., _Vector boson production at hadron colliders: a fully exclusive QCD calculation at NNLO_ , Phys. Rev. Lett. 103 (2009) 082001, arXiv:0903.2120
  * [13] K. Melnikov and F. Petriello, _Electroweak gauge boson production at hadron colliders through $O(\alpha_{s}^{2})$_, Phys. Rev. D74 (2006) 114017, arXiv:hep-ph/0609070
  * [14] J. Rojo et al., _The PDF4LHC report on PDFs and LHC data: Results from Run I and preparation for Run II_ , J. Phys. G42 (2015) 103103, arXiv:1507.00556
  * [15] J. Butterworth et al., _PDF4LHC recommendations for LHC Run II_ , J. Phys. G43 (2016) 023001, arXiv:1510.03865
  * [16] A. Kusina et al., _Vector boson production in pPb and PbPb collisions at the LHC and its impact on nCTEQ15 PDFs_ , Eur. Phys. J. C77 (2017) 488, arXiv:1610.02925
  * [17] K. J. Eskola, P. Paakkinen, H. Paukkunen, and C. A. Salgado, _EPPS16: Nuclear parton distributions with LHC data_ , Eur. Phys. J. C77 (2017) 163, arXiv:1612.05741
  * [18] D. de Florian and R. Sassot, _Nuclear parton distributions at next-to-leading order_ , Phys. Rev. D69 (2004) 074028, arXiv:hep-ph/0311227
  * [19] M. Hirai, S. Kumano, and T.-H. Nagai, _Determination of nuclear parton distribution functions and their uncertainties in next-to-leading order_ , Phys. Rev. C76 (2007) 065207, arXiv:0709.3038
  * [20] S. Atashbar Tehrani, _Nuclear parton densities and their uncertainties at the next-to-leading order_ , Phys. Rev. C86 (2012) 064301
  * [21] H. Khanpour and S. Atashbar Tehrani, _Global analysis of nuclear parton distribution functions and their uncertainties at next-to-next-to-leading order_ , Phys. Rev. D93 (2016) 014026, arXiv:1601.00939
  * [22] S. J. Brodsky, D. S. Hwang, and I. Schmidt, _Final state interactions and single spin asymmetries in semi-inclusive deep inelastic scattering_ , Phys. Lett. B530 (2002) 99, arXiv:hep-ph/0201296
  * [23] J. C. Collins, _Leading twist single transverse-spin asymmetries: Drell-Yan and deep inelastic scattering_ , Phys. Lett. B536 (2002) 43, arXiv:hep-ph/0204004
  * [24] D. Boer, P. J. Mulders, and F. Pijlman, _Universality of T odd effects in single spin and azimuthal asymmetries_ , Nucl. Phys. B667 (2003) 201, arXiv:hep-ph/0303034
  * [25] A. V. Belitsky, X. Ji, and F. Yuan, _Final state interactions and gauge invariant parton distributions_ , Nucl. Phys. B656 (2003) 165, arXiv:hep-ph/0208038
  * [26] J.-C. Peng and J.-W. Qiu, _Novel phenomenology of parton distributions from the Drell–Yan process_ , Prog. Part. Nucl. Phys. 76 (2014) 43, arXiv:1401.0934
  * [27] LHCb collaboration, R. Aaij et al., _Observation of $Z$ production in proton-lead collisions at LHCb_, JHEP 09 (2014) 030, arXiv:1406.2885
  * [28] CMS collaboration, V. Khachatryan et al., _Study of Z boson production in pPb collisions at $\sqrt{s_{\scriptscriptstyle\text{NN}}}=5.02$ TeV_, Phys. Lett. B759 (2016) 36, arXiv:1512.06461
  * [29] ATLAS collaboration, G. Aad et al., _$Z$ boson production in $p+$Pb collisions at $\sqrt{s_{\scriptscriptstyle\text{NN}}}=5.02$ TeV measured with the ATLAS detector_, Phys. Rev. C92 (2015) 044915, arXiv:1507.06232
  * [30] ALICE collaboration, J. Adam et al., _W and Z boson production in p-Pb collisions at $\sqrt{s_{\scriptscriptstyle\text{NN}}}$ = 5.02 TeV_, JHEP 02 (2017) 077, arXiv:1611.03002
  * [31] ALICE collaboration, S. Acharya et al., _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_, JHEP 09 (2020) 076, arXiv:2005.11126
  * [32] CMS collaboration, A. M. Sirunyan et al., _Study of Drell-Yan dimuon production in proton-lead collisions at $\sqrt{s_{\mathrm{NN}}}=$ 8.16 TeV_, JHEP 05 (2021) 182, arXiv:2102.13648
  * [33] K. Kovarik et al., _nCTEQ15 - Global analysis of nuclear parton distributions with uncertainties in the CTEQ framework_ , Phys. Rev. D93 (2016) 085037, arXiv:1509.00792
  * [34] M. Vesterinen and T. R. Wyatt, _A Novel Technique for Studying the Z Boson Transverse Momentum Distribution at Hadron Colliders_ , Nucl. Instrum. Meth. A602 (2009) 432, arXiv:0807.4956
  * [35] LHCb collaboration, A. A. Alves Jr. et al., _The LHCb detector at the LHC_, JINST 3 (2008) S08005
  * [36] LHCb collaboration, R. Aaij et al., _LHCb detector performance_ , Int. J. Mod. Phys. A30 (2015) 1530022, arXiv:1412.6352
  * [37] T. Sjöstrand, S. Mrenna, and P. Skands, _A brief introduction to PYTHIA 8.1_ , Comput. Phys. Commun. 178 (2008) 852, arXiv:0710.3820
  * [38] T. Sjöstrand, S. Mrenna, and P. Skands, _PYTHIA 6.4 physics and manual_ , JHEP 05 (2006) 026, arXiv:hep-ph/0603175
  * [39] J. Pumplin et al., _New generation of parton distributions with uncertainties from global QCD analysis_ , JHEP 07 (2002) 012, arXiv:hep-ph/0201195
  * [40] I. Belyaev et al., _Handling of the generation of primary events in Gauss, the LHCb simulation framework_ , J. Phys. Conf. Ser. 331 (2011) 032047
  * [41] D. J. Lange, _The EvtGen particle decay simulation package_ , Nucl. Instrum. Meth. A462 (2001) 152
  * [42] P. Golonka and Z. Was, _PHOTOS Monte Carlo: A precision tool for QED corrections in $Z$ and $W$ decays_, Eur. Phys. J. C45 (2006) 97, arXiv:hep-ph/0506026
  * [43] T. Pierog et al., _EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider_ , Phys. Rev. C92 (2015) 034906, arXiv:1306.0121
  * [44] Geant4 collaboration, J. Allison et al., _Geant4 developments and applications_ , IEEE Trans. Nucl. Sci. 53 (2006) 270
  * [45] Geant4 collaboration, S. Agostinelli et al., _Geant4: A simulation toolkit_ , Nucl. Instrum. Meth. A506 (2003) 250
  * [46] M. Clemencic et al., _The LHCb simulation application, Gauss: Design, evolution and experience_, J. Phys. Conf. Ser. 331 (2011) 032023
  * [47] P. Nason, _A New method for combining NLO QCD with shower Monte Carlo algorithms_ , JHEP 11 (2004) 040, arXiv:hep-ph/0409146
  * [48] S. Frixione, P. Nason, and C. Oleari, _Matching NLO QCD computations with Parton Shower simulations: the POWHEG method_ , JHEP 11 (2007) 070, arXiv:0709.2092
  * [49] S. Alioli, P. Nason, C. Oleari, and E. Re, _A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX_ , JHEP 06 (2010) 043, arXiv:1002.2581
  * [50] S. Alioli, P. Nason, C. Oleari, and E. Re, _NLO vector-boson production matched with shower in POWHEG_ , JHEP 07 (2008) 060, arXiv:0805.4802
  * [51] A. Buckley et al., _LHAPDF6: parton density access in the LHC precision era_ , Eur. Phys. J. C75 (2015) 132, arXiv:1412.7420
  * [52] S. Dulat et al., _New parton distribution functions from a global analysis of quantum chromodynamics_ , Phys. Rev. D93 (2016) 033006, arXiv:1506.07443
  * [53] D. Stump et al., _Inclusive jet production, parton distributions, and the search for new physics_ , JHEP 10 (2003) 046, arXiv:hep-ph/0303013
  * [54] G. Dujany and B. Storaci, _Real-time alignment and calibration of the LHCb Detector in Run II_ , J. Phys. Conf. Ser. 664 (2015) 082010
  * [55] R. Aaij et al., _The LHCb trigger and its performance in 2011_, JINST 8 (2013) P04022, arXiv:1211.3055
  * [56] R. Aaij et al., _Tesla: an application for real-time data analysis in High Energy Physics_ , Comput. Phys. Commun. 208 (2016) 35, arXiv:1604.05596
  * [57] S. J. Brodsky and H. J. Lu, _Shadowing and Antishadowing of Nuclear Structure Functions_ , Phys. Rev. Lett. 64 (1990) 1342
  * [58] European Muon collaboration, J. J. Aubert et al., _The ratio of the nucleon structure functions $F_{2}^{N}$ for iron and deuterium_, Phys. Lett. B123 (1983) 275
  * [59] LHCb collaboration, R. Aaij et al., _Measurement of the forward $Z$ boson cross-section in $p$ $p$ collisions at $\sqrt{s}=$7 TeV_, JHEP 08 (2015) 039, arXiv:1505.07024
  * [60] CDF collaboration, F. Abe et al., _A Measurement of $\sigma B(W\rightarrow e\nu)$ and $\sigma B(Z^{0}\rightarrow e^{+}e-)$ in $\bar{p}p$ collisions at $\sqrt{s}=1800$ GeV_, Phys. Rev. D44 (1991) 29
  * [61] LHCb collaboration, R. Aaij et al., _Measurement of the track reconstruction efficiency at LHCb_ , JINST 10 (2015) P02007, arXiv:1408.1251
  * [62] T. Adye, _Unfolding algorithms and tests using RooUnfold_ , in PHYSTAT 2011,  (Geneva), 313–318, CERN, 2011, arXiv:1105.1160
  * [63] G. D’Agostini, _A Multidimensional unfolding method based on Bayes’ theorem_ , Nucl. Instrum. Meth. A362 (1995) 487
  * [64] LHCb collaboration, R. Aaij et al., _Measurement of forward $W$ and $Z$ boson production in $p$ $p$ collisions at $\sqrt{s}=$8 TeV_, JHEP 01 (2016) 155, arXiv:1511.08039
  * [65] LHCb collaboration, R. Aaij et al., _Measurement of the forward $Z$ boson production cross-section in $p$ $p$ collisions at $\sqrt{s}=$13 TeV_, JHEP 09 (2016) 136, arXiv:1607.06495
  * [66] A. Höcker and V. Kartvelishvili, _SVD approach to data unfolding_ , Nucl. Instrum. Meth. A372 (1996) 469, arXiv:hep-ph/9509307
  * [67] G. Cowan, Statistical data analysis, Oxford university press, 1998


# Dual-cycle Constrained Bijective VAE-GAN For Tagged-to-Cine Magnetic
Resonance Image Synthesis

###### Abstract

Tagged magnetic resonance imaging (MRI) is a widely used imaging technique for
measuring tissue deformation in moving organs. Due to tagged MRI’s intrinsic
low anatomical resolution, another matching set of cine MRI with higher
resolution is sometimes acquired in the same scanning session to facilitate
tissue segmentation, thus adding extra time and cost. To mitigate this, in
this work, we propose a novel dual-cycle constrained bijective VAE-GAN
approach to carry out tagged-to-cine MR image synthesis. Our method is based
on a variational autoencoder backbone with cycle reconstruction constrained
adversarial training to yield accurate and realistic cine MR images given
tagged MR images. Our framework has been trained, validated, and tested using
1,768, 416, and 1,560 subject-independent paired slices of tagged and cine MRI
from twenty healthy subjects, respectively, demonstrating superior performance
over the comparison methods. Our method can potentially be used to reduce the
extra acquisition time and cost, while maintaining the same workflow for
further motion analyses.

Index Terms—  Tagged MRI, image synthesis, deep learning, Generative
adversarial networks.

##  1 Introduction

Assessment of an internal organ’s deformation has been an important topic in
both medical imaging research and clinical practice. One of the most popular
modalities to measure internal motion is tagged magnetic resonance (MR)
imaging (MRI) with spatially encoded tag patterns [1]. Tagged MR images with
horizontal and vertical tag patterns are acquired as the internal organ, such
as the tongue, moves while the tag patterns deform together with the organ;
recording all deformation information in the deformed tag patterns. These tag
deformations can then be analyzed with motion extraction methods including
Harmonic Phase (HARP) based methods, recovering motion fields in either two-
dimensional (2D) or three-dimensional (3D) spaces [2, 3]. Tagged MRI is an
established method for internal organ motion imaging [4] that continues to
enjoy widespread use.

In practice, due to the nature of tagged MRI’s encoding method, image
resolution on the anatomy of acquired images is relatively low. Although
motion fields are computed in each whole MR slice, it is usually necessary to
separate the motion field to leave only the part within the organ of interest.
As a result, another separate set of dynamic cine MR images without tags are
sometimes acquired alongside the tagged MR images as a matching pair [5]. Cine
MR images have higher resolution than the tagged MR images and are used for
segmenting the region of interest, such as the whole tongue [6], in which to
compute the motion field. In addition, cine MR images are used to assess
surface motion of the organs. However, this extra acquisition can double the
time and cost.

To address this, we systematically investigate a variational autoencoders
(VAEs)-based generative model with adversarial training objectives to
synthesize cine MRI from tagged MRI as shown in Fig. 1. To the best of our
knowledge, this is the first attempt at tagged-to-cine MR image synthesis.
Although the L1 loss in VAEs uses pixel-wise pair similarity to provide a good
baseline, the results are often blurry [7, 8, 9]. Therefore, an additional
generative adversarial network (GAN) [10, 11] is further incorporated into our
model to enrich the texture, while mitigating potential differences in texture
or contrast within real cine MR images. Considering that image translation is
largely unconstrained and even ill-posed [12], the anatomical structure of a
specific input-output pair is not necessarily consistent if we only align the
distributions with GANs. To address this, we introduce an unpaired cycle
reconstruction constraint [12, 13] into our VAE-based model so that our
proposed dual-cycle constrained bijective VAE-GAN can maintain the anatomical
structure.

Both quantitative and qualitative evaluation results using a total of 3,774
paired slices of tagged and cine MRI from twenty healthy subjects show the
validity of our proposed dual-cycle constrained bijective VAE-GAN framework
and its superiority to conventional VAEs and GANs based image style
translation methods.

##  2 Methodology

Given the paired tagged and cine MR images $\\{x^{t},x^{c}\\}$, we propose to
learn a parameterized mapping $f:x^{t}\rightarrow\tilde{x}^{t\rightarrow c}$
from the tagged MR images $x^{t}$ to the generated cine MR images
$\tilde{x}^{t\rightarrow c}$ to closely resemble the real cine MR images
$x^{c}$.

###  2.1 VAE-based synthesis

A straightforward baseline structure of $f$ can be the VAE, which is
constructed with an encoder $Enc$ and a decoder $Dec$. The VAE first maps
$x^{t}$ to a code $z$ in a latent space via the $Enc$, and then decodes $z$ to
reconstruct the target image $\tilde{x}^{t\rightarrow c}$ via the $Dec$. For
the reconstruction error, we adopt the $L_{1}$ loss, which can be formulated
as $L_{1}({\tilde{x}^{t\rightarrow c}},x^{c})=|{\tilde{x}^{t\rightarrow
c}}-{\mathbf{x}}^{c}|$. The KL-divergence $L_{KL}$ between the latent code $z$
and the prior standard Gaussian is modeled with the standard
reparameterization trick [14]. The training objective is to minimize

$L_{VAE}=L_{1}({\widetilde{x}^{t\rightarrow c}},x^{c})+\alpha
L_{KL}(Enc(x^{t})||\mathcal{N}(0,I)),$ (1)

where $\alpha$ is a balancing weight for the KL divergence term to penalize
the deviation of the distribution of the latent code from the prior
distribution.

A limitation of VAEs, however, is that the generated images tend to be blurry
[7, 8] due to the injected noise, limited expressiveness of the inference
models, or imperfect element-wise criteria such as the $L_{1}$ or $L_{2}$ loss
[15]. Although recent studies [14] have greatly improved the predicted log-
likelihood, the VAE image generation quality still lags behind GANs [7, 16,
17].

In order to enforce perceptual realism with respect to the anatomic structure
and improve the quality of the generated textures, we adopt the additional
adversarial training procedure [18].

###  2.2 Dual-cycle constrained bijective VAE-GAN

The GAN training aims to learn a mapping $f$ with the output $f(x^{t})$ that
is indistinguishable from real cine MR images by an adversary trained to
classify $f(x^{t})$ apart from $x^{c}$. A possible solution to combine the
objectives of GAN and VAE is following the VAE-GAN framework [15] to configure
a discriminator to discriminate real samples from both the reconstructions and
the generated examples with sampling Gaussian noise. Similarly, the Pix2Pix
[19] proposes to use the real-fake pair as input to the discriminator to
enhance the stability of GANs. In theory, the GAN’s objective can induce an
output distribution over $f(x^{t})$ that matches the empirical distribution of
$x^{c}$. The optimal $f$ thereby translates the domain of $x^{t}$ to a domain
$f(x^{t})$ distributed identically to $x^{c}$.

However, such a translation does not guarantee that the individual input
$x^{t}$ and output $\tilde{x}^{t\rightarrow c}=f(x^{t})$ are paired up in a
meaningful way. There are infinitely many mappings $f$ that will induce the
same distribution over $f(x^{t})$ [12]. While the generated cine MR images may
be visually pleasing and realistic, the anatomical shape or structure may not
be necessarily consistent with the input tagged MR images. We note that the
objective of our tagged-to-cine MR image synthesis is to facilitate the
segmentation of the internal organ or to observe surface motion of the organ;
therefore, we introduce a stricter requirement with respect to the structural
consistency over conventional image translation approaches [12, 13].

Fig. 1: Illustration of the proposed dual-cycle constrained bijective VAE-GAN.
Encoder, decoder, and discriminator have two parallel modules for $x^{t}$ and
$x^{c}$, respectively. Note that only the gray masked subnets are used for
testing.

Recently, cycleGAN and its variants [12, 13] have shown outstanding
performance in many unpaired image translation tasks. They enforce the
generated output with $f:x^{t}\rightarrow\tilde{x}^{t\rightarrow c}$ can be
mapped back to the original image with an additional inverse translator
$f^{-1}:x^{c}\rightarrow\tilde{x}^{c\rightarrow t}$. With two bijective
mappings, the distribution of $\tilde{x}^{t\rightarrow c}$ can be largely
constrained and therefore the anatomic structure can be better maintained
compared with the conventional approaches. Accordingly, we propose to
introduce the unpaired cycle constraint [12, 13] for our VAE-based pair-wise
synthesis. Our framework, as illustrated in Fig. 1, is based on a dual-branch
VAEs and cycle-constrained GANs.

Fig. 2: Comparison of different tagged-to-cine MR generation methods,
including the vanilla VAE*, VAE+vanilla GAN [15]*, Pix2Pix [19]*, and our
proposed dual-cycle constrained bijective VAE-GAN. * indicates the first
attempt at tagged-to-cine MR image synthesis.

Detailed structure: The two encoders (i.e., $Enc_{t}$ and $Enc_{c}$) have the
same structure and take tagged MR images $x^{t}$ or cine MR images $x^{c}$ as
input, respectively. The two decoders (i.e., $Dec_{t}$ and $Dec_{c}$) mirror
the structure of the encoder responsible for tagged or cine MRI generation.
The latent space of both tagged and cine MRI is shared [13].

Similarly, there are two discriminators (i.e., $Dis_{t}$ and $Dis_{c}$) for
the adversarial training of tagged or cine MRI. For tagged MR images sampled
from the tagged MRI domain, $Dis_{t}$ should output true, while for images
generated by $Dec_{t}$, it should output false.

$Dec_{t}$ can generate two types of images: 1) images from the reconstruction
stream $\tilde{x}^{t\rightarrow t}=Dec_{t}(z\sim Enc_{t}(x^{t}))$ and 2)
images from the translation stream $\tilde{x}^{t\rightarrow t}=Dec_{t}(z\sim
Enc_{c}(x^{c}))$. Since the reconstruction stream is trained in a supervised
manner, it suffices that we only apply the adversarial training to images from
the translation stream, $\tilde{x}^{c\rightarrow t}$, given by

$\displaystyle L_{1}^{t}(\tilde{x}^{t\rightarrow
t},x^{t})=|Dec_{t}(Enc_{t}(x^{t}))-x^{t}|$ (2) $\displaystyle
L_{Dis_{t}}=\text{log}(Dis_{t}(x^{t}))+\text{log}(1-Dis_{t}(\tilde{x}^{c\rightarrow
t})).$ (3)

We apply similar processing to $Dec_{c}$, where $Dis_{c}$ is trained to output
true for real images sampled from cine MRI and false for images generated from
$Dec_{c}$, given by

$\displaystyle L_{1}^{c}(\tilde{x}^{c\rightarrow
c},x^{c})=|Dec_{c}(Enc_{c}(x^{c}))-x^{c}|$ (4) $\displaystyle
L_{Dis_{c}}=\text{log}(Dis_{c}(x^{c}))+\text{log}(1-Dis_{c}(\tilde{x}^{t\rightarrow
c})).$ (5)

Training Strategy: We jointly optimize the learning objectives of our backbone
VAE and the adversarial cycle-reconstruction streams:

${}^{\text{\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ min}}_{Enc_{t},Dec_{t}}\leavevmode\nobreak\
L_{VAE}^{t}+\beta L_{1}^{t}(\tilde{x}^{t\rightarrow t},x^{t})+\lambda
L_{Dis_{t}},$ (6) ${}^{\text{\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
min}}_{Enc_{c},Dec_{c}}\leavevmode\nobreak\ L_{VAE}^{c}+\beta
L_{1}^{c}(\tilde{x}^{c\rightarrow c},x^{t})+\lambda L_{Dis_{c}},$ (7)
$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
{}^{\text{max}}_{Dis_{t}}\leavevmode\nobreak\ L_{Dis_{t}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
^{\text{max}}_{Dis_{c}}\leavevmode\nobreak\ L_{Dis_{c}},$ (8)

where $L_{VAE}^{t}=L_{1}({\widetilde{x}^{t\rightarrow c}},x^{c})+\alpha
L_{KL}(Enc_{t}(x^{t})||\mathcal{N}(0,I))$, and
$L_{VAE}^{c}=L_{1}({\widetilde{x}^{c\rightarrow t}},x^{t})+\alpha
L_{KL}(Enc_{c}(x^{c})||\mathcal{N}(0,I))$ are the backbone VAE objectives for
$x^{t}$ and $x^{c}$, respectively. The hyper-parameters $\beta$ and $\lambda$
control the weights of the cycle-based reconstruction objective term and the
discriminate term, respectively. Inheriting from GAN, the $Enc$-$Dec$ and
$Dis$ play the round-based min-max adversarial game to improve each other to
find a saddle point. In practice, we usually sample a similar number of
$x^{t}$ and $x^{c}$ in training [13].

Testing Tagged-to-Cine Translation: Targeting our tagged-to-cine MR synthesis
objective, after training, we can obtain the translation functions by
assembling a subset of the subnetworks, i.e., $Enc_{t}$ and $Dec_{c}$.
Therefore, given tagged MR images $x^{t}$, we can generate its corresponding
$\tilde{x^{t\rightarrow c}}=Dec_{c}(Enc_{t}(x^{t}))$. The involved subnets are
indicated with the gray mask in Fig. 1.

##  3 EXPERIMENTS AND RESULTS

###  3.1 Data Acquisition and Subjects

For the experiments performed in this study, a total of 3,774 paired tagged
and cine MR images from twenty healthy subjects were acquired while speaking
an utterance, “a souk” in line with a periodic metronome-like sound. MRI
scanning was performed on a Siemens 3.0T TIM Trio system with a 12-channel
head coil and a 4-channel neck coil using a segmented gradient echo sequence
[3]. The field of view was 240$\times$240 mm with an in-plane resolution of
1.87$\times$1.87 mm and a slice thickness of 6 mm. The image sequence was
obtained at the rate of 26 fps. Both cine and tagged MRI are in the same
spatiotemporal coordinate space.

###  3.2 Experimental Setup

We used 10 subjects (1,768 slice pairs) for training, 2 subjects (416 slice
pairs) for hyper-parameter validation, and 8 subjects (1,560 slice pairs) for
evaluation. Tagged MR images with horizontal tag patterns were used in our
evaluation. For fair comparison, we resized the tagged and cine MR images to
256$\times$256 and adopted the $Enc$, $Dec$, and $Dis$ backbones from [19] for
all of the methods.

Our framework was implemented using the PyTorch deep learning toolbox. The
training was performed on an NVIDIA V100 GPU, which took about 4 hours. In
practice, translating one tagged MR image to a cine MR image in testing with
only $Enc_{t}$ and $Dec_{c}$ took about 0.1 seconds.

In order to align the absolute value of each loss, we set different weights
for each part with grid searching on the validation data. Specifically, we set
weight $\alpha=1$, $\beta=1$, and $\lambda=0.5$. We used Adam optimizer for
training. The learning rate was set at $lr_{Enc,Dec}=1\mathrm{e}{-3}$ and
$lr_{Dis}=1\mathrm{e}{-4}$ and the momentum was set at 0.5.

###  3.3 Qualitative Evaluations

The synthesis results using VAE [14], VAE-GAN [15], Pix2Pix [19], and our
proposed method are shown in Fig. 2. The proposed framework successfully
synthesized the cine MR images, which are consistent with the target ground
truth cine MR images. The VAE only model simply relies on the pixel-wise
matching, thus resulting in relatively blurry outputs. There is a relatively
large difference in contrast and texture within the tongue region between the
real cine MR images, which could affect subsequent analyses. Thus, the added-
on adversarial objective is expected to enrich the fine-grained tissue texture
[18]. VAE-GAN [15] adopts the vanilla GAN, which is difficult to optimize, and
often leads to the well-known problem of mode collapse. By using the pair-wise
inputs to the discriminator, the model collapse in Pix2Pix [19] was relatively
alleviated. Although the generative outputs were visually realistic, the
unconstrained image translation yielded distortion in the anatomical
structure. In contrast, our proposed dual-cycle constrained bijective VAE-GAN
with the cycle reconstruction constraint generated visually pleasing results
with better structural consistency, which is important for subsequent
analyses.

Table 1: Numerical comparisons of four methods in testing across 1,560 slice
pairs. The best and the second best results are bold and underlined,
respectively.

Methods | L1 $\downarrow$SSIM $\uparrow$ | PSNR $\uparrow$ | IS $\uparrow$ |
---|---|---|---|---
VAEonly |  148.8$\pm$0.1 | 0.9486$\pm$0.0014 | 32.68$\pm$0.06 | 11.13$\pm$0.11
VAE-GAN [15] | 151.4$\pm$0.2 | 0.9507$\pm$0.0013 | 34.14$\pm$0.05 | 11.68$\pm$0.17
Pix2Pix [19] | 150.2$\pm$0.2 |  0.9612$\pm$0.0011 |  36.81$\pm$0.07 |  13.77$\pm$0.15
Proposed |  149.6$\pm$0.2 |  0.9746$\pm$0.0015 |  38.72$\pm$0.07 |  15.58$\pm$0.12

###  3.4 Quantitative Evaluations

The synthesized images were expected to have realistic-looking textures, and
to be structurally coherent with its corresponding ground-truth $x^{c}$. For
quantitative evaluation, we adopted widely used evaluation metrics: mean L1
error, structural similarity index measure (SSIM), peak signal-to-noise ratio
(PSNR), and inception score (IS) [20].

Table 1 lists numerical comparisons between the proposed framework, VAE [14],
VAE-GAN [15], and Pix2Pix [19] for the 8 testing subjects. The proposed dual-
cycle constrained bijective VAE-GAN outperformed the other three comparison
methods with respect to SSIM, PSNR, and IS. We note that all of the compared
methods have the L1 minimization objective.

##  4 Conclusion and Future Work

We presented a novel framework to synthesize cine MRI from its paired tagged
MRI. We systematically investigated VAE with the GAN objectives, and proposed
a novel dual-cycle constrained bijective VAE-GAN. The adversarial training
introduced the realistic texture, while mitigating potential differences in
contrast and texture within the real cine MR images. In addition, importantly,
the qualitative evaluation shows the anatomical structure of the tongue was
better maintained with the additional cycle reconstruction constraint. Our
experimental results demonstrated that our method surpassed the comparison
methods as quantitatively and qualitatively assessed. The synthesized cine MRI
can potentially be used for further segmenting the tongue and observing
surface motion, and, in our future work, we will develop a segmentation
network in conjunction with the proposed network using the synthesized cine MR
images.

##  5 COMPLIANCE WITH ETHICAL STANDARDS

The subjects signed an informed consent form, and data were collected in
accordance with the protocol approved by the Institutional Review Board (IRB)
of the University of Maryland, Baltimore. The retrospective use of the data
was granted by IRB of Massachusetts General Hospital.

##  6 ACKNOWLEDGMENTS

This work is supported by NIH R01DC014717, R01DC018511, R01CA133015, and
P41EB022544.

## References

  * [1] Nael F Osman, Elliot R McVeigh, and Jerry L Prince,  “Imaging heart motion using harmonic phase mri,”  IEEE TMI, vol. 19, no. 3, pp. 186–202, 2000.
  * [2] Fangxu Xing, Jonghye Woo, Arnold D Gomez, Dzung L Pham, Philip V Bayly, Maureen Stone, and Jerry L Prince,  “Phase vector incompressible registration algorithm for motion estimation from tagged magnetic resonance images,”  IEEE TMI, vol. 36, no. 10, pp. 2116–2128, 2017.
  * [3] Fangxu Xing, Jonghye Woo, Junghoon Lee, Emi Z Murano, Maureen Stone, and Jerry L Prince,  “Analysis of 3-D tongue motion from tagged and cine magnetic resonance images,”  Journal of Speech, Language, and Hearing Research, vol. 59, no. 3, pp. 468–479, 2016.
  * [4] Caroline Petitjean, Nicolas Rougon, and Philippe Cluzel,  “Assessment of myocardial function: a review of quantification methods and results using tagged mri,”  Journal of Cardiovascular Magnetic Resonance, vol. 7, no. 2, pp. 501–516, 2005.
  * [5] Vijay Parthasarathy, Jerry L Prince, Maureen Stone, Emi Z Murano, and Moriel NessAiver,  “Measuring tongue motion from tagged cine-mri using harmonic phase (harp) processing,”  The Journal of the Acoustical Society of America, pp. 491–504, 2007\.
  * [6] Junghoon Lee, Jonghye Woo, Fangxu Xing, Emi Z Murano, Maureen Stone, and Jerry L Prince,  “Semi-automatic segmentation of the tongue for 3D motion analysis with dynamic mri,”  in ISBI, 2013, p. 1465.
  * [7] Ian Goodfellow,  “Nips 2016 tutorial: Generative adversarial networks,”  arXiv:1701.00160, 2016.
  * [8] Xiaofeng Liu, BVK Vijaya Kumar, Ping Jia, and Jane You,  “Hard negative generation for identity-disentangled facial expression recognition,”  Pattern Recognition, vol. 88, pp. 1–12, 2019.
  * [9] Xiaofeng Liu, Xing Fangxu, Yang Chao, C.C.-Jay Kuo, Georges El Fakhri, and Jonghye Woo,  “Symmetric-constrained irregular structure inpainting for brain mri registration with tumor pathology,”  in MICCAI BrainLes, 2020.
  * [10] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio,  “Generative adversarial nets,”  in NeurIPS, 2014, pp. 2672–2680.
  * [11] Karim Armanious, Chenming Jiang, Marc Fischer, Thomas Küstner, Tobias Hepp, Konstantin Nikolaou, Sergios Gatidis, and Bin Yang,  “Medgan: Medical image translation using gans,”  Computerized Medical Imaging and Graphics, vol. 79, pp. 101684, 2020\.
  * [12] Jun-Yan Zhu, Taesung Park, Phillip Isola, and Alexei A Efros,  “Unpaired image-to-image translation using cycle-consistent adversarial networks,”  in ICCV, 2017, pp. 2223–2232.
  * [13] Ming-Yu Liu, Thomas Breuel, and Jan Kautz,  “Unsupervised image-to-image translation networks,”  in NeurIPS, 2017, pp. 700–708.
  * [14] Durk P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and Max Welling,  “Improved variational inference with inverse autoregressive flow,”  in NeurIPS, 2016, pp. 4743–4751.
  * [15] Anders Boesen Lindbo Larsen, Søren Kaae Sønderby, Hugo Larochelle, and Ole Winther,  “Autoencoding beyond pixels using a learned similarity metric,”  in International conference on machine learning. PMLR, 2016, pp. 1558–1566.
  * [16] Xiaofeng Liu, Tong Che, Yiqun Lu, Chao Yang, Site Li, and Jane You,  “Auto3d: Novel view synthesis through unsupervisely learned variational viewpoint and global 3d representation,”  in ECCV, 2020, pp. 1358–1376.
  * [17] Xiaofeng Liu, Site Li, Lingsheng Kong, Wanqing Xie, Ping Jia, Jane You, and BVK Kumar,  “Feature-level frankenstein: Eliminating variations for discriminative recognition,”  in CVPR, 2019, pp. 637–646.
  * [18] Ian Goodfellow, Yoshua Bengio, Aaron Courville, and Yoshua Bengio,  Deep Learning, vol. 1,  MIT press Cambridge, 2016.
  * [19] Phillip Isola, Jun-Yan Zhu, Tinghui Zhou, and Alexei A Efros,  “Image-to-image translation with conditional adversarial networks,”  in CVPR, 2017, pp. 1125–1134.
  * [20] Yuhang Song, Chao Yang, Zhe Lin, Xiaofeng Liu, Qin Huang, Hao Li, and C-C Jay Kuo,  “Contextual-based image inpainting: Infer, match, and translate,”  in ECCV, 2018, pp. 3–19.