\section{Introduction}

In the framework of collinear factorization \cite{Collins:1989gx}, the differential cross section
for a given process in hadronic collisions can be calculated in terms of
convolutions of process-independent parton distribution functions (PDFs) and process-specific perturbatively calculable coefficient functions. By combining data for inclusive deep-inelastic scattering (DIS) from HERA, fixed-target experiments, or neutrino scattering, with various LHC data, including e.g. (di-)jet, top-pair and electroweak-boson production, proton PDFs have by now reached accuracies at the level of a few percent or better over wide kinematic regions. This is vital for precision studies of processes in the Standard-Model and beyond (see \cite{Gao:2017yyd} and references therein). 

One may also apply the same factorization to processes involving high-energy nuclei, such as e+A or p+A. This provides information on nuclear parton distribution functions (nPDFs). Being fundamental properties of atomic nuclei, nPDFs are of much importance for our general understanding of the strong interactions. At  the same time, extracting ``cold-nuclear matter'' effects from collisions with a small projectile and a target nucleus also provides a baseline for quark-gluon plasma studies. 

The backbone for analyses of nPDFs have been fixed-target DIS e+A data, often accompanied by Drell-Yan (DY) dilepton production data. In addition to DIS with charged lepton beams, also a good amount of neutrino DIS data are nowadays available that provide some additional sensitivity to the flavor dependence of nuclear effects. In recent years, data from p+Pb collisions from the LHC have fulfilled the promise of extending the kinematic reach and further constraining the flavor dependence of the nPDFs. Prime examples are the $\mathrm{D}$-meson production at forward rapidities measured by LHCb \cite{LHCb:2017yua} which is particularly sensitive to small-$x$ gluon nPDFs \cite{Eskola:2019bgf}, and $W^{\pm}$ and $Z$ boson production data from ATLAS \cite{Aad:2015gta} and CMS \cite{Khachatryan:2015hha, Khachatryan:2015pzs}, which are primarily sensitive to different quark-flavor combinations at moderate values of $x$ and provide potential constraints for flavor dependence \cite{Paukkunen:2010qg} and also for strange quarks and gluons \cite{Kusina:2020lyz}. Furthermore, dijet production \cite{CMS:2018jpl} can offer further constraints on gluon nPDFs in the intermediate $x$-region \cite{Eskola:2019dui}. Other possible observables include inclusive direct-photon and pion production in p+Pb collisions at the LHC. The former has been recently measured by ATLAS \cite{ATLAS:2019ery} and for the latter there are already several data sets available from ALICE \cite{ALICE:2018vhm, ALICE:2016dei, ALICE:2021est}.

Numerous nPDF analyses are available by now which can be classified in terms of the data used in the analysis and of their perturbative precision. Early ``global'' analyses include the EPS09~\cite{Eskola:2009uj} and DSSZ~\cite{deFlorian:2011fp} sets, which were based on charged-lepton and neutrino DIS and DY data from fixed-target experiments, and inclusive pion production data from d+Au collisions at RHIC. The EPPS16 analysis \cite{Eskola:2016oht} was the first to include also data from p+Pb collisions at the LHC by analyzing Run~I data for dijets and $W^{\pm}$ and $Z$ boson production. The original nCTEQ15 analysis \cite{Kovarik:2015cma} did not include any LHC data, but has recently been extended to include e.g. Run~II $W^{\pm}$ and $Z$ data from p+Pb collisions (nCTEQ15wz analysis \cite{Kusina:2020lyz}). Similar data sets have been considered also by the NNPDF collaboration in their most recent analysis nNNPDF2.0 \cite{AbdulKhalek:2020yuc}. The nCTEQ15 analysis has been recently updated to contain also single-inclusive hadron production \cite{Duwentaster:2021ioo} and also other active groups have recently prepared new analyses. The EPPS21 \cite{Eskola:2021nhw} and nNNPDF3.0 \cite{Khalek:2022zqe} analyses both include a significant amount of new LHC data including e.g. dijet, $W^{\pm}$ and $Z$ boson and inclusive D-meson production from the LHCb. All analyses mentioned so far have been performed at the next-to-leading-order (NLO) of perturbative QCD. There are a few recent analyses that have been performed at next-to-next-to-leading order (NNLO): nNNPDF1.0 \cite{AbdulKhalek:2019mzd}, TUJU19 \cite{Walt:2019slu} and KSASG20 \cite{Khanpour:2020zyu}. So far, these have only considered charged-lepton and neutrino DIS, and fixed-target DY data.

In this work we present new nuclear PDF analyses at NLO and NNLO, using the same framework as in our previous analysis, TUJU19, but now including in addition to neutral-current DIS with a lepton beam and charged-current neutrino DIS data also new electroweak (EW)-boson production data from the LHC. This is the first time where LHC data are employed in a full NNLO analysis for nuclear PDFs. We include the LHC data also for a proton ``baseline'' fit. In this way, we obtain a fully consistent way of computing cross sections for p+$A$ collisions at NLO or NNLO. The resulting nPDFs may also be readily used for further cross section calculations at NNLO in nuclear collisions. As an application we study EW-boson production in Pb+Pb collisions and compare our results to recent ATLAS and CMS data \cite{ATLAS:2019ibd, ATLAS:2019maq,CMS:2021kvd} which, contrary to expectations, have previously been found to be difficult to describe within a factorization-based NNLO calculation \cite{Eskola:2020lee}. 

Our paper is organized as follows. We describe the theoretical framework for our analysis in Sec. \ref{sec:theory}, then discuss the analysis procedure in Sec. \ref{sec:analysis} and the selection of experimental data in Sec. \ref{sec:expdata}. The results of the analysis are presented in Sec. \ref{sec:results}. We discuss the application to Pb+Pb collisions at the LHC in Sec. \ref{sec:applications}. Finally, we summarize our work in Sec. \ref{sec:summary}, where also an outlook toward future developments is given.

\section{Theoretical framework}
\label{sec:theory}

The theoretical framework is very similar to that adopted in our earlier analysis, TUJU19 \cite{Walt:2019slu}, which includes also a detailed description of neutral-current (NC) and charged-current (CC) DIS processes. Here we give an overview of the calculational framework that we have set up to use the new LHC data.

\subsection{Drell-Yan and $W^{\pm}$, $Z$ boson production processes}
\label{sec:crosssections}

As new constraints we include LHC data for inclusive hadroproduction of electro-weak (EW) bosons, $W^{\pm}$ and $Z/\gamma^{*}$. Often such processes are referred to as Drell-Yan (DY) processes, originally relating to reactions where two hadrons collide and form a highly virtual photon that decays to a lepton pair \cite{Drell:1970wh}. The signature of such an event is a lepton pair ($l\bar{l}$ with $l=e,\,\mu,\,\tau$) of the same flavor with a large invariant mass. In leading order (LO) such a process can take place only by annihilation of a quark-antiquark pair of the same flavor,
\begin{equation}
q\bar{q} \rightarrow \gamma^* \rightarrow l\bar{l}\,.
\end{equation}
 
At high enough energies, similar scattering may form also other EW bosons, namely $Z$ and $W^{\pm}$. As the former has the same quantum numbers as a photon, the production and decay channels are the same as in traditional DY. In the latter case, however, the LO process requires e.g. a $u\bar{d}$ initial state and the leptonic decay channel to include production of a charged lepton and a corresponding neutrino,
\begin{align}
q\bar{q} &\rightarrow Z \rightarrow l\bar{l}\,,\quad \text{and} \notag \\
q_i\bar{q_j} &\rightarrow W^{\pm} \rightarrow l\nu_{l}\,.
\end{align}
As the (anti)neutrinos from the $W^{-(+)}$ decays cannot be directly measured in the detectors and only the momentum of the charged lepton is known, some further kinematical cuts are required to increase the sensitivity to the signal process.

A factorization theorem for DY has been explicitly proved to all orders in perturbation theory \cite{Bodwin:1984hc, Collins:1985ue, Collins:1988ig}. In order to later be able to assess the impact of DY (or EW-boson production) data on our analysis, it is useful to write down the double differential DY cross section at LO \cite{ellis_qcd_1996, sterman_handbook_2001}:
\begin{equation}
\frac{\mathrm{d}^2\sigma}{\mathrm{d}M^2 \mathrm{d}y} = \frac{4\pi\alpha^2 \tau}{3 N_C M^4} \sum_q e_q^2 \left[ q(x_1) \bar{q}(x_2)+\bar{q}(x_1) q(x_2) \right]
\label{eq-dy}
\end{equation}
where $N_C$ is the number of colors, $\alpha$ the fine structure constant, $\tau=M^2/s$ (with $\sqrt{s}$ the center-of-mass energy), and we have introduced the variables
\begin{align}
M^2 &= (p_1+p_2)^2 \equiv \hat{s} = x_1 x_2 s & &(\text{invariant mass})\,, \notag \\
y &= \frac{1}{2}\log\frac{x_1}{x_2} & &(\text{rapidity}).
\label{eq-dy-variables}
\end{align}
Here $p_{1,2}$ are momenta of the incoming partons carrying momentum fraction $x_{1,2}$ of the total hadron momentum. In terms of $M^2$ and $y$ we have (again, to LO):
\begin{align}
&x_1=\frac{M}{\sqrt{s}}\exp(y) & &\mathrm{and} & &x_2=\frac{M}{\sqrt{s}}\exp(-y)\,.
\label{eq-dy-variables2}
\end{align}
These relations can be used to estimate the sensitivity to different regions of momentum fraction when studying EW boson production at specific kinematics. Often in the experimental analysis, however, pseudorapidity $\eta$ is used instead of $y$ as the former does not require information of the mass of the given system. 

In Eq.~(\ref{eq-dy}), $q(x_1)$ and $\bar{q}(x_2)$ are the quark and antiquark PDFs, whose scale dependence we have omitted. Comparing to the LO DIS cross section formula, we notice that in DY we always have a combination of quark and antiquark PDFs. Thus we expect DY processes to provide more information on the sea quark densities.

The well-known factorized all-order expression for the Drell-Yan cross section is based on convolutions of the PDFs with perturbative hard-scattering functions. Beyond LO, there are also contributions from processes other than $q\bar{q}^{(\prime)}$ annihilation. For example, starting from NLO, also contributions from gluon-(anti)quark initial states arise, and at NNLO and beyond gluon-gluon scattering participates. Denoting the partonic hard-scattering function for a given partonic channel $ab\to {\mathrm{EW\,boson}}+X$ by $\omega_{ab}$, we have the perturbative expansion
\begin{align}
\omega_{ab}\,=\,\omega^{(0)}_{ab}+\frac{\alpha_s}{\pi}\,\omega^{(1)}_{ab}
+\left(\frac{\alpha_s}{\pi}\right)^2 \omega^{(2)}_{ab}+{\cal O}(\alpha_s^3)\,,
\end{align}
where $\alpha_s$ is the strong coupling evaluated at some renormalization scale $\mu_R$.
The NLO coefficient functions $\omega^{(1)}_{ab}$ have been known for a long time~\cite{Altarelli:1979ub}. The NNLO corrections $\omega^{(2)}_{ab}$ are also fully available~\cite{Hamberg:1990np,Harlander:2002wh,Anastasiou:2003yy,Anastasiou:2003ds,Catani:2009sm,Gavin:2010az}, and their inclusion is a new feature of our TUJU21 analysis.

\subsection{Input parameterization}

Apart from the perturbative hard-scattering functions, another key ingredient in a global analysis of collinear PDFs is their (Dokshitzer–Gribov–Lipatov–Altarelli–Parisi) DGLAP evolution \cite{Gribov:1972ri, Lipatov:1974qm,  Altarelli:1977zs, Dokshitzer:1977sg}, which describes how the PDFs depend on the factorization scale. The perturbative order of the splitting kernels to be used in the evolution equations needs to be chosen in accordance with the order in the hard-scattering functions. As only the perturbative scale evolution of the PDFs is given by the evolution equations, a non-perturbative input at some initial scale is required to obtain a PDF set. In this analysis the baseline parton distributions of a proton are parametrized similarly as in \cite{Walt:2019slu}:
\begin{equation}
xf^{p}_i\left(x,Q_0^2 \right) = c_0\,x^{c_1} (1-x)^{c_2} \left(1+c_3\,x + c_4\,x^2 \right),
\label{pdf-parameterization}
\end{equation}
where the index $i$ runs over parton flavor $i=g, \,d_{\mathrm{v}}, \,u_{\mathrm{v}}, \,\bar{u}, \,\bar{d}, \,\bar{s}$, with the subscript ``$\mathrm{v}$'' refering to the up or down valence-quark contribution. We have kept the flavour dependence in the valence sector but for the sea quarks we had to assume $\bar{u} = \bar{d} = \bar{s} = s$ as it turned out that within the applied data and the adopted framework it is difficult to have a converged fit without such a condition. We define the input parameterization at the charm mass threshold, $Q_0^2=m_c^2=1.69\,\mathrm{GeV^2}$. We do not include any intrinsic charm content at this scale and apply the FONLL-A(C) general-mass variable-flavor-number-scheme \cite{Cacciari:1998it, Forte:2010ta} to calculate the heavy-quark cross sections at $Q>Q_0$ in our NLO (NNLO) analysis. 

To obtain the PDFs for protons bound in a nucleus, we add dependence on the nuclear mass number $A$ to the parameters $c_i$ by re-parameterizing them as
\begin{equation}
c_k\,\rightarrow c_k(A) = c_{k,0}+c_{k,1}\left( 1 - A^{-c_{k,2}} \right),
\label{coeff-A}
\end{equation}
where $k={0,\dots,4}$. A similar form has been used also in the nCTEQ15 analysis \cite{Kovarik:2015cma, Kusina:2020lyz}. It is worth pointing out that with $A=1$ the $A$-dependent right-hand part of Eq.~(\ref{coeff-A}) vanishes and the free proton PDFs are recovered when indentifying $c_{k,0} = c_{k}$. In order to obtain the full nuclear PDF, we need to add also the contribution from neutrons in the nucleus. The PDFs of neutrons are not fitted separately, but are determined from the proton PDFs based on isospin symmetry, giving $u^{\mathrm{n}/A} = d^{\mathrm{p}/A}$ and $d^{\mathrm{n}/A} = u^{\mathrm{p}/A}$, and likewise for the light antiquarks. Then, an average PDF for a nucleon bound in a nucleus with $Z$ protons and $A-Z$ neutrons is obtained as
\begin{equation}
f_i^{\,N/A} \left( x,Q^{\,2} \right) = \frac{Z\cdot f_i^{\,p/A}+ (A-Z)\cdot f_i^{\,n/A}}{A}\,.
\label{nucleon}
\end{equation}


\section{Analysis procedure}
\label{sec:analysis}

\subsection{Minimization procedure}

At the heart of PDF fitting is the minimization of the difference between the experimental data and the calculated cross sections, providing the optimal PDF set within the adopted framework. In this work, as well as in our previous analysis \cite{Walt:2019slu}, this is done by minimizing ${\large{\chi}}^2$ defined as
\begin{align}
{\large{\chi}}^2\left(\textbf{m},\,\textbf{b}\right) =& \sum_i \, \frac{\left[ m_i - \Sigma_{\alpha}\,\gamma^i_{\alpha}\,\mu_i\,b_{\alpha} - \mu_i \right]^2}{\left( \delta_{i,\mathrm{stat}} \sqrt{\mu_i m_i} \right)^2 + \left( \delta_{i,\mathrm{uncorr}} m_i \right)^2 } \notag \\
+& \sum_{\alpha} b^2_{\alpha}.
\label{eq-chi2-scaled}
\end{align}
Here, $\mu_i$ is the value of the measured data point for a given observable, whereas $m_i$ is the actual theoretical value calculated using DGLAP-evolved PDFs with given parameters $\{c_k\}$. The uncertainties are represented by the $\delta_{i,\mathrm{stat}}$ and $\delta_{i,\mathrm{uncorr}}$, which are the relative statistical and uncorrelated systematic uncertainties, respectively, as well as by the $\gamma^{i}_{\alpha} \mu_i$, which are the correlated errors. Furthermore, the $b_{\alpha}$ are the so-called nuisance parameters that are determined during the fitting. To account for the correlated systematic uncertainties for each data set, we allow for shifts of the calculated cross section within the quoted uncertainty, penalizing such shifts by the additional $b^2_{\alpha}$ contributions to ${\large{\chi}}^2$.

The minimization of ${\large{\chi}}^2$, Eq.~(\ref{eq-chi2-scaled}), provides a central set of PDFs with the parameter values allowing the best description of the used data. Since the experimental data always contain several uncertainties, as listed above, a separate error analysis needs to be performed to study how these uncertainties propagate into the fitted PDFs. In this QCD analysis the Hessian method \cite{Pumplin:2001ct,Pumplin:2000vx} is used for the analysis of the uncertainties. The Hessian error analysis is performed assuming a quadratic expansion of the function ${\large{\chi}}^2 = {\large{\chi}}_0^2 + \Delta {\large{\chi}}^2$ around its global minimum. Here, ${\large{\chi}}_0^2$ is the value of the function at the global minimum (with the best-fit parameters $\{k_0\}$) and $\Delta {\large{\chi}}^2$ is the displacement from the minimum \cite{Pumplin:2001ct,Pumplin:2000vx}. During the performed error analysis $\Delta {\large{\chi}}^2$ defines the tolerance criterion determining the allowed growth of ${\large{\chi}}^2$. As argued in our previous work \cite{Walt:2019slu}, for the proton baseline with $13$ free fit parameters we select $\Delta {\large{\chi}}^2 = 20$ and for the nuclear PDF error analysis we choose $\Delta {\large{\chi}}^2 = 50$ for our 16 free parameters.

\subsection{The fitting framework}

Our global analyses of the baseline proton and nuclear PDFs are performed with the \textsc{xFitter} \cite{Zenaiev:2016jnq,Bertone:2017tig} tool. The main goal of the \textsc{xFitter} project is to provide an open-source tool to fit proton PDFs with varied theoretical assumptions. In order to perform a nuclear PDF analysis several modifications of the code were performed in the context of the work presented in Ref.~\cite{Walt:2019slu}. 

When going from DIS to collisions between two hadrons where cross section calculations involve convolutions of more than one PDF, the required calculations become computationally expensive, and even more so when increasing the perturbative precision to NNLO. Therefore, to implement data for DY and EW-boson production beyond LO, it is necessary to prepare fast interpolation grids allowing for efficient comparisons with the data when the PDF parameters are iterated. Standard tools to handle such convolutions with interpolation grids, such as \textsc{APPLgrid} \cite{Carli:2010rw}, can be linked with \textsc{xFitter}, facilitating fast cross section calculations using the pre-computed grids. To prepare the actual grids used in the fitting, a code capable of calculating the given production cross section needs to be linked with an interpolation code. Depending on the available data files, further processing in \textsc{xFitter} requires grids in a ROOT format \cite{Brun:1997pa}. It is also possible to use $K$-factors when going from NLO to NNLO cross sections to reduce the computing effort. We have used such $K$-factors in our proton baseline fit in cases where they were available in xFitter for a given data set. The settings and PDFs used to calculate the $K$-factors can be found from the references provided in Table~\ref{tab-expdatady-proton}.

The schematic overview of the fitting routine and the required tool set (\textsc{xFitter} and all other modules) as applied for the \textit{nuclear} PDF fit is shown in Fig.~\ref{fig-dy-nPDF}. For this, a new convolution routine had to be implemented in \textsc{xFitter} where one first needs to prepare fast interpolation grids by calculating the observables in the relevant kinematic region, using an existing PDF set, e.g. in the LHAPDF6 format \cite{Buckley:2014ana}. In this work, MCFM 8.3 \cite{Falkowski:2012cu, Campbell:2015qma, Boughezal:2016wmq} has been used to calculate the interpolation grids for EW boson production after suitable modification to handle asymmetric collisions as needed for p+Pb. The obtained grids were then used as the input for the optimization routine. This setup provides the possibility to use fast interpolation grids in a plain text format generated with \textsc{MCFM} up to order NNLO, without relying on approximative $K$-factors. During the fitting procedure the actual theoretical predictions were obtained by convoluting the fitted nPDFs based on updated parameters with the pre-calculated differential cross sections provided in form of grids. For the convolution step one needs to specify the proton PDF baseline that was used for the generation of the fast interpolation grids in MCFM. In our case, we are adopting our own proton PDF baseline prepared in the form of an LHAPDF6 library. 
\begin{figure*}[h!]
\centering
\includegraphics[width=0.9\textwidth]{xFitter-Schema_DY_nPDFs.png}
\caption{Schematic view of the high-level \textsc{xFitter} functionalities as relevant for Drell-Yan and $W^{\pm}$, $Z$ boson production processes, with the newly implemented convolution routine in \textsc{xFitter} as used for the nPDF fit TUJU21. \textsc{xFitter} logo credited from~\cite{xfitter_link}.}
\label{fig-dy-nPDF}
\end{figure*}


\subsection{Experimental data}
\label{sec:expdata}

We build upon our previous analysis, TUJU19, including all the same charged-lepton and neutrino DIS data as before. On top of these we include data for DY and EW boson production for both our proton baseline and the nuclear PDFs. 

The DIS data used for the proton baseline fit are summarized in Table~\ref{tab-expdata-proton}, also showing the resulting ${\large{\chi}}^2$ values at NLO and NNLO obtained in the analysis. The newly included experimental data from DY and $W^{\pm}$ and $Z$-boson production for the proton baseline are listed in Table \ref{tab-expdatady-proton}. For the experimental proton data used, the fast interpolation grids were publicly available and the details of grid generation for each proton PDF data set can be found from the references provided in Table~\ref{tab-expdata-proton}. These data include high- and low-mass DY data from ATLAS at $\sqrt{s} = 7~\text{TeV}$ \cite{Aad:2013iua, Aad:2014qja}, EW boson production data from ATLAS at $\sqrt{s} = 7~\text{TeV}$ \cite{Aad:2012sb}, and increased-luminosity data for these from ATLAS at $\sqrt{s} = 7~\text{TeV}$ \cite{Aaboud:2016btc}. In addition, also $W^{\pm}$ production cross sections from CMS at $\sqrt{s} = 8~\text{TeV}$ have been included. In total the newly included data sets consist of 134 data points.
\begin{center}
\renewcommand{\arraystretch}{1.25}
\begin{table}[tb!]
\caption{Summary of experimental DIS data used to determine our proton PDF baseline. In the last two columns the ${\large{\chi}}^2$ values at NLO and NNLO obtained in our analysis are provided.}
\label{tab-expdata-proton}
\scriptsize
\begin{tabular}{l|l|c|r|c|r|r}
\toprule
  Exp. & Data set & Year & Ref. & $N_{\mathrm{dp}}$  & ${\large{\chi}}^2$ NLO  & ${\large{\chi}}^2$ NNLO  \\ \hline \hline
  BCDMS & F2p 100GeV & 1996 & \cite{Benvenuti:1989rh} & 83 & 96.30  & 93.69 \\
        & F2p 120GeV & & & 90   & 70.54 & 68.70 \\ 
        & F2p 200GeV & & & 79   & 91.81  & 86.32 \\ 
        & F2p 280GeV  & & & 75   & 67.52  & 69.71 \\ \hline
  HERA 1+2 & NCep 920 & 2015 & \cite{Abramowicz:2015mha} & 377 & 459.71 & 482.23 \\\ 
           & NCep 820 & & &	70    & 72.91 & 73.47\\
           & NCep 575 & & &	254   & 222.64 & 231.35 \\
           & NCep 460 & & &	204   & 218.84 & 225.68 \\
           & NCem & & &	159   & 227.80  & 232.79 \\
           & CCep & & &	39   & 46.59  & 43.17 \\
           & CCem & & &	42   & 60.49 & 63.60 \\ \hline
  NMC-97 & NCep & 1997 & \cite{Arneodo:1996qe}  & 100 & 117.72 & 111.31 \\ \hline \hline            
   \multicolumn{4}{r|}{\textbf{In total:}} & 1559 &  &  \\
\toprule  
\end{tabular} 
\end{table}
\end{center}

\begin{center}
\renewcommand{\arraystretch}{1.25}
\begin{table}[tb!]
\caption{Summary of experimental DY and $W^{\pm},\,Z$ boson production data used to determine our proton PDFs. In the last two columns the ${\large{\chi}}^2$ values at NLO and NNLO obtained in our analysis are provided.}
\label{tab-expdatady-proton}
\scriptsize
\begin{tabular}{l|l|c|r|c|r|r}
\toprule
  Exp. & Data set & Year & Ref. & $N_{\mathrm{dp}}$  & ${\large{\chi}}^2$ NLO  & ${\large{\chi}}^2$ NNLO  \\ \hline \hline
  ATLAS DY & high mass DY & 2013 & \cite{Aad:2013iua} & 13 & 10.91  & 11.54 \\
        & low mass DY & 2014 & \cite{Aad:2014qja} &	 8  & 22.86  & 8.68 \\ \hline  
  ATLAS $W^{\pm},\,Z$ &  $W^+$ lepton $\eta$ & 2012 & \cite{Aad:2012sb} & 11 & 15.77 & 14.00 \\ 
            & $W^-$ lepton $\eta$ & & & 11  & 7.98 & 8.54 \\ 
             & $Z$ $y$ & & &	8  & 4.10 & 2.71 \\
        & $W^+$ lepton $y$ & 2016 & \cite{Aaboud:2016btc} & 11  & 19.57 & 10.71 \\  
             & $W^-$ lepton $y$ & & &	 11 & 11.11 & 11.82 \\ 
             & high mass CC $Z$ $y$ & & &	6 & 7.61 & 6.05 \\ 
            & high mass CF $Z$ $y$ & & &	6 & 3.92 & 3.95  \\ 
            &  low mass $Z$ $y$ & & &	6 & 41.32 & 23.84 \\
            & peak CC $Z$ $y$ & & & 12 & 45.76 & 14.40 \\
           & peak CF $Z$ $y$ & & & 9 & 21.85 & 7.55 \\ \hline
  CMS $W^{\pm}$ & $W^+$ $\sigma$ 8 TeV & 2016 & \cite{Khachatryan:2016pev} & 11 & 10.64 & 4.81 \\ 
           &  $W^-$ $\sigma$ 8 TeV & & &	11   & 8.14 & 9.27 \\ \hline \hline
   \multicolumn{4}{r|}{\textbf{In total:}} & 134 &  &    \\
\toprule  
\end{tabular} 
\end{table}
\end{center}

Table \ref{tab-expdata} provides a list of nuclear-DIS data as also used in the TUJU19 analysis, but with the ${\large{\chi}}^2$ values obtained in this analysis. The input data files and the fast interpolation grids used for the fitting procedure laid out in Fig.~\ref{fig-dy-nPDF} for nPDFs have been collected and prepared as part of this analysis for the newly included data points summarized in Table \ref{tab-expdatady-nPDFs}. The added data include Run~I measurements of $Z$ boson production in p+Pb collisions at the LHC by ATLAS \cite{Aad:2015gta} and CMS \cite{Khachatryan:2015pzs} at $\sqrt{s_{\mathrm{NN}}} = 5.02~\text{TeV}$ and the more recent Run~II measurement of $W^{\pm}$ boson production in p+Pb collisions at $\sqrt{s_{\mathrm{NN}}} = 8.16~\text{TeV}$ by CMS \cite{CMS:2019leu}. In total these add 74 data points to the nPDF analysis. In addition to these, there are more EW-boson data available from the LHC experiments. In particular, there are data from ALICE \cite{ALICE:2016rzo} and LHCb \cite{LHCb:2014jgh} that extend the kinematic reach but have only two data points per set and suffer from large statistical uncertainties. There are also {Run-I} data for $W^{\pm}$ production from CMS \cite{CMS:2015ehw} at the same kinematics than the more recent data but with significantly larger uncertainties. There are also fixed-target DY data available e.g. from E772 \cite{Alde:1990im}, E866 \cite{NuSea:1999egr} and CLAS \cite{CLAS:2019vsb} experiments that have been used previously in similar analyses. As the grid computations at NNLO are computationally very heavy, we included only the datasets that we expect to provide the strongest constraints for the nPDFs in the selected set up. In section \ref{sec:DYpPb} we consider the recent Run-II CMS measurement for $Z/\gamma^*$ production for which it has proven difficult to obtain $\chi^2/N_{\mathrm{dp}}$ values close to unity in a NLO QCD analysis \cite{Eskola:2021nhw, Khalek:2022zqe}. As we have discussed in the context of the TUJU19 nPDF analysis, some authors have found tension between the neutrino- and charged-lepton DIS data. To check for potential tension with the new LHC data we have also performed fits without any neutrino-DIS data and found that the new data was equally well described as when the neutrino data were included. Thus we do not find any tension between these data sets.
\begin{center}
\renewcommand{\arraystretch}{1.25}
\begin{table}[tb!]
\caption{{Summary of experimental DIS data used to determine the nuclear PDFs. In the last two columns the ${\large{\chi}}^2$ values at NLO and NNLO obtained in our analysis are provided.}}
\label{tab-expdata}
\scriptsize
\begin{tabular}{l|l|c|r|c|r|r}
\toprule
  Nucleus & Exp. & Year & Ref. & $N_{\mathrm{dp}}$  & ${\large{\chi}}^2$ NLO  & ${\large{\chi}}^2$ NNLO  \\   
\hline
\hline
  D & NMC 97 & 1996  & \cite{Arneodo:1996qe} &  120 & 151.61  & 121.52  \\ 
       & EMC 90 & 1989  & \cite{Arneodo:1989sy} & 21 & 24.31 & 22.89  \\ \hline
  He/D  & HERMES  & 2002    & \cite{Airapetian:2002fx} & 7 &  6.79 & 8.92 \\ 
      & NMC 95, re.  & 1995   & \cite{Amaudruz:1995tq} &  13 & 10.67  & 10.42  \\   
      & SLAC E139 & 1994  & \cite{Gomez:1993ri} & 11 & 6.47  & 4.42 \\ \hline
   Li/D  & NMC 95 & 1995  & \cite{Arneodo:1995cs} &  12 & 9.10  & 9.00  \\ \hline
   Be/D  &  SLAC E139 & 1994  & \cite{Gomez:1993ri} & 10 & 11.58  & 11.51  \\ \hline    
   Be/C   & NMC 96 & 1996  & \cite{Arneodo:1996rv} & 14 & 13.56  & 16.06  \\ \hline
   C    & EMC 90 & 1989  & \cite{Arneodo:1989sy} & 17 &  13.41 & 13.44  \\
   C/D  & FNAL E665 & 1995   & \cite{Adams:1995is} & 3  & 2.00  & 1.83 \\ 
        & SLAC E139 & 1994   & \cite{Gomez:1993ri} & 6 &  20.69 & 13.86  \\
        & EMC 88 & 1988  & \cite{Ashman:1988bf} & 9 & 3.70  & 4.22  \\  
        & NMC 95, re.  & 1995    & \cite{Amaudruz:1995tq} & 13 &  34.96 & 19.49  \\ \hline
   C/Li  & NMC 95, re. & 1995    & \cite{Amaudruz:1995tq} & 10 & 7.77  & 10.18 \\ \hline
   N/D  & HERMES &  2002    & \cite{Airapetian:2002fx} & 1 &  0.95 & 2.08  \\ \hline
   Al/D     & SLAC E139 & 1994  & \cite{Gomez:1993ri} &  10 & 18.49  &  9.49 \\\hline    
   Al/C    & NMC 96 &  1996   & \cite{Arneodo:1996rv} & 14   & 7.29  & 6.23   \\\hline 
   Ca &    EMC 90  & 1989  & \cite{Arneodo:1989sy} & 19 & 13.41  & 13.44  \\
   Ca/D  & NMC 95, re. & 1995    & \cite{Amaudruz:1995tq} & 12 & 34.75  & 21.82 \\ 
     &    FNAL E665 &  1995  & \cite{Adams:1995is} & 3  &  1.84 & 2.41  \\
     &    SLAC E139  &  1994  & \cite{Gomez:1993ri} & 6  & 16.74  & 9.01 \\                   
   Ca/Li  & NMC 95, re. & 1995  & \cite{Amaudruz:1995tq} &  10 & 1.45  &  1.33  \\ \hline
   Ca/C  & NMC 95, re. & 1995  & \cite{Amaudruz:1995tq}  & 10  & 9.35  & 8.00 \\ 
     &  NMC 96 &  1996  & \cite{Arneodo:1996rv} & 14  & 8.45  &  6.42  \\ \hline 
   Fe  &    SLAC E140  & 1993 & \cite{Dasu:1993vk} & 2  & 0.14  & 0.04 \\
   Fe/D  &    SLAC E139  & 1994  & \cite{Gomez:1993ri} & 14 & 44.53  & 32.07  \\
   Fe/C  & NMC 96 & 1996  & \cite{Arneodo:1996rv} & 14  &  11.17 & 9.93 \\ \hline     
   $\nu$ Fe    & CDHSW  & 1991 & \cite{Berge:1989hr}  &  464  & 404.26  &  358.19  \\ 
   $\bar{\nu}$ Fe    & CDHSW  & 1991  & \cite{Berge:1989hr} & 462  & 439.53  &  395.99  \\ \hline
   Cu/D  & EMC 88  &  1988 & \cite{Ashman:1988bf}  & 9  & 8.38  & 5.71  \\
        & EMC 93 & 1993   & \cite{Ashman:1992kv} & 19  & 26.38  &  12.58 \\ \hline
   Kr/D  & HERMES  & 2002   & \cite{Airapetian:2002fx}  & 1   & 2.02  &  2.02 \\ \hline   
   Ag/D  & SLAC E139  & 1994  & \cite{Gomez:1993ri} & 6  & 21.37  & 18.80 \\ \hline   
   Sn/D  & EMC 88 &  1988 & \cite{Ashman:1988bf}  & 8 & 14.37 & 13.98 \\ \hline   
   Sn/C  & NMC 96 & 1996 & \cite{Arneodo:1996rv}  & 14 & 6.48  & 8.52 \\
        & NMC 96, $Q^2$ dep. & 1996  & \cite{Arneodo:1996ru} & 134 & 76.13  & 75.03  \\ \hline     
   Xe/D   & FNAL E665 & 1992  & \cite{Adams:1992nf} &  3 & 1.64  & 1.34 \\ \hline
   Au/D  & SLAC E139  & 1994  & \cite{Gomez:1993ri} &  11 & 16.89  &  18.66  \\ \hline 
   Pb/D  & FNAL E665  & 1995  & \cite{Adams:1995is} & 2  &  8.28 & 7.72   \\ \hline 
   Pb/C   & NMC 96 &  1996  & \cite{Arneodo:1996rv} & 14   & 8.32  &  5.42  \\ \hline
   $\nu$ Pb   & CHORUS  & 2005 & \cite{Onengut:2005kv} &  405   & 259.48  &  237.85 \\ 
   $\bar{\nu}$ Pb   & CHORUS  & 2005  & \cite{Onengut:2005kv} & 405  & 356.01  & 352.09  \\ \hline \hline 
   \multicolumn{4}{r|}{\textbf{In total:}} & 2336 &  &  \\
\toprule
\end{tabular}
\end{table}
\end{center}

\begin{center}
\renewcommand{\arraystretch}{1.25}
\begin{table}[tb!]
\caption{Summary of experimental $W^{\pm},\,Z$ boson production data from LHC p+Pb collisions in Run I and Run II used to determine our nuclear PDFs. In the last two columns the ${\large{\chi}}^2$ values at NLO and NNLO obtained in our analysis are provided.}
\label{tab-expdatady-nPDFs}
\scriptsize
\begin{tabular}{l|l|l|c|r|c|r|r}
\toprule
  Nucleus & Proc. & Exp. & Year & Ref. & $N_{\mathrm{dp}}$  & ${\large{\chi}}^2$ NLO  & ${\large{\chi}}^2$ NNLO  \\ 
\hline
\hline
pPb & $Z$  & LHC Run I ATLAS  &  2015 & \cite{Aad:2015gta}  & 14  &  16.40 & 12.82 \\  
& $Z$  & LHC Run I CMS  & 2015  & \cite{Khachatryan:2015pzs} & 12  &  8.76 & 7.30 \\
& $W^{-}$  & LHC Run II CMS  &  2019 & \cite{CMS:2019leu}  & 24  & 39.58 & 42.83 \\   
& $W^{+}$  & LHC Run II CMS  &   &   & 24  &  41.08 & 39.07 \\ \hline \hline  
   \multicolumn{5}{r|}{\textbf{In total:}} & 74 &  &    \\
\end{tabular}
\end{table}
\end{center}

\section{Results}
\label{sec:results}

\subsection{Proton baseline}
\label{sec:protonpdfs}

Analyses of nuclear PDFs have often been performed by using an existing proton PDF set as a baseline for the nuclear modifications. In this work, however, we have fitted the proton PDFs using the same setup as for the nuclear PDFs. This ensures that all assumptions like sum rules, parton flavor decomposition, etc., as well as all parameters like coupling constants and quark masses, and also further settings like e.g. the heavy flavor mass scheme, are applied in a consistent way. Furthermore, this paves the way for a future combined analysis of proton and nuclear PDFs.

In this section the updated free proton PDF sets in TUJU21 are compared to those of our earlier TUJU19 analysis, which were determined using DIS data only. The free proton PDFs used as a baseline for the nuclear part of the QCD analysis were updated by including experimental data for DY, $W^{\pm},\,Z$ boson production processes taken by the ATLAS and CMS collaborations at the LHC; see Sec. \ref{sec:expdata} for details.
The comparisons are presented in Fig. \ref{proton_TUJU21_NLO} for NLO and in Fig. \ref{proton_TUJU21_NNLO} for NNLO.

The impact of the newly added LHC data is rather mild at NLO. We observe that the uncertainties of the PDFs have become slightly smaller in some cases, especially for the valence quarks. At NNLO, the resulting distributions for the valence quarks, and especially for gluons are slightly decreased with respect to our previous analysis. The results obtained for the updated free proton PDF baseline confirm that DY, and $W^{\pm},\,Z$ boson production data can be accommodated together with the DIS data and provide further constraints in a global analysis. Using these data to pin down the proton PDFs in the same framework as the nuclear ones will ensure that the baseline is well constrained in the region where new data are included for the nPDFs.

The parameters for the input distributions for our best fit of the proton baseline are collected in Appendix~\ref{app-pdf-params}.

\begin{figure*}[tb!]
     \begin{center}
          \subfigure{
              \includegraphics[width=0.237\linewidth]{proton-19vs21_nlo_f_g_1_69.pdf}} 
              %
          \subfigure{    
              \includegraphics[width=0.237\linewidth]{proton-19vs21_nlo_f_ubar_1_69.pdf}} 
             
          \subfigure{                           
              \includegraphics[width=0.237\linewidth]{proton-19vs21_nlo_f_uv_1_69.pdf}} 
             
          \subfigure{                                
              \includegraphics[width=0.237\linewidth]{proton-19vs21_nlo_f_dv_1_69.pdf}} 
              %
          \subfigure{              
              \includegraphics[width=0.237\textwidth]{proton-19vs21_nlo_f_g_100.pdf}} 
              %
          \subfigure{              
              \includegraphics[width=0.237\textwidth]{proton-19vs21_nlo_f_ubar_100.pdf}} 
             
          \subfigure{                        
              \includegraphics[width=0.237\textwidth]{proton-19vs21_nlo_f_uv_100.pdf}} 
             
          \subfigure{                                   
              \includegraphics[width=0.237\textwidth]{proton-19vs21_nlo_f_dv_100.pdf}} 
              %
                            %
          \end{center} 
    \caption{Proton baseline PDFs in TUJU21 at NLO compared to the previous TUJU19 results, shown at the initial scale $Q_0^2=1.69\,\mathrm{GeV}^2$ (upper panels) and at $Q^2=100\,\mathrm{GeV}^2$ (lower panels) after DGLAP evolution.}
\label{proton_TUJU21_NLO}    
    \end{figure*}

\begin{figure*}[tb!]
     \begin{center}
          \subfigure{
              \includegraphics[width=0.237\textwidth]{proton-19vs21_nnlo_f_g_1_69.pdf}} 
              %
          \subfigure{    
              \includegraphics[width=0.237\textwidth]{proton-19vs21_nnlo_f_ubar_1_69.pdf}} 
             
          \subfigure{                           
              \includegraphics[width=0.237\textwidth]{proton-19vs21_nnlo_f_uv_1_69.pdf}} 
             
          \subfigure{                                
              \includegraphics[width=0.237\textwidth]{proton-19vs21_nnlo_f_dv_1_69.pdf}} 
              %
          \subfigure{              
              \includegraphics[width=0.237\textwidth]{proton-19vs21_nnlo_f_g_100.pdf}} 
              %
          \subfigure{              
              \includegraphics[width=0.237\textwidth]{proton-19vs21_nnlo_f_ubar_100.pdf}} 
             
          \subfigure{                        
              \includegraphics[width=0.237\textwidth]{proton-19vs21_nnlo_f_uv_100.pdf}} 
             
          \subfigure{                                   
              \includegraphics[width=0.237\textwidth]{proton-19vs21_nnlo_f_dv_100.pdf}} 
              %
                            %
          \end{center} 
    \caption{Same as for Fig. \ref{proton_TUJU21_NLO}, but at NNLO.}
\label{proton_TUJU21_NNLO}    
    \end{figure*}


\subsection{Nuclear PDFs}

\begin{figure*}[b!]
     \begin{center}
          \subfigure{
              \includegraphics[width=0.237\linewidth]{Pb-19vs21_nlo_f_g_1_69.pdf}} 
              %
          \subfigure{    
              \includegraphics[width=0.237\linewidth]{Pb-19vs21_nlo_f_ubar_1_69.pdf}} 
             
          \subfigure{                           
              \includegraphics[width=0.237\linewidth]{Pb-19vs21_nlo_f_uv_1_69.pdf}} 
             
          \subfigure{                                
              \includegraphics[width=0.237\linewidth]{Pb-19vs21_nlo_f_dv_1_69.pdf}} 
              %
          \subfigure{              
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nlo_f_g_100.pdf}} 
              %
          \subfigure{              
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nlo_f_ubar_100.pdf}} 
             
          \subfigure{                        
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nlo_f_uv_100.pdf}} 
             
          \subfigure{                                   
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nlo_f_dv_100.pdf}} 
              %
          \subfigure{              
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nlo_R_g_100.pdf}} 
              %
          \subfigure{              
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nlo_R_ubar_100.pdf}} 
             
          \subfigure{                        
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nlo_R_uv_100.pdf}} 
             
          \subfigure{                                   
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nlo_R_dv_100.pdf}} 
                            %
                            %
          \end{center} 
    \caption{NLO nuclear parton distribution functions in TUJU21 for a lead nucleus, compared to the previous TUJU19 results, shown at the initial scale $Q_0^2=1.69\,\mathrm{GeV}^2$ (upper panels) and at $Q^2=100\,\mathrm{GeV}^2$ after DGLAP evolution (center panels). The lower panels show
    the corresponding ratios of PDFs for a proton bound in lead over the free proton PDFs.}
\label{Pb_TUJU21_NLO}    
    \end{figure*}

\begin{figure*}[b!]
     \begin{center}
          \subfigure{
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nnlo_f_g_1_69.pdf}} 
              %
          \subfigure{    
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nnlo_f_ubar_1_69.pdf}} 
             
          \subfigure{                           
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nnlo_f_uv_1_69.pdf}} 
             
          \subfigure{                                
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nnlo_f_dv_1_69.pdf}} 
              %
          \subfigure{              
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nnlo_f_g_100.pdf}} 
              %
          \subfigure{              
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nnlo_f_ubar_100.pdf}} 
             
          \subfigure{                        
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nnlo_f_uv_100.pdf}} 
             
          \subfigure{                                   
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nnlo_f_dv_100.pdf}} 
              %
          \subfigure{              
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nnlo_R_g_100.pdf}} 
              %
          \subfigure{              
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nnlo_R_ubar_100.pdf}} 
             
          \subfigure{                        
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nnlo_R_uv_100.pdf}} 
             
          \subfigure{                                   
              \includegraphics[width=0.237\textwidth]{Pb-19vs21_nnlo_R_dv_100.pdf}} 
                            %
          \end{center} 
    \caption{Same as for Fig. \ref{Pb_TUJU21_NLO}, but at NNLO.}
\label{Pb_TUJU21_NNLO}    
    \end{figure*}

The resulting nuclear PDFs, referred to as TUJU21, are presented in Figs. \ref{Pb_TUJU21_NLO} at NLO and \ref{Pb_TUJU21_NNLO} at NNLO for the bound-proton-in-lead-nucleus PDFs, including also ratios to our baseline free proton PDFs. We also compare to  the nPDFs of our previous DIS-only analysis TUJU19. At NLO, the largest differences between the two analyses occur for gluons and sea quarks. For gluons the small-$x$ suppression is significantly milder than in TUJU19, along with a slightly reduced uncertainty. Also, the strong antishadowing enhancement at intermediate values of $x$ we found previously is now tamed to a more moderate $\sim 10\%$ effect. At the initial scale of the fit the sea quark nPDFs are now slightly lower in the small-$x$ region, with a somewhat smaller uncertainty band, but have remained very similar at larger $x$. Because of the larger gluon nPDF in the updated fit, the sea quark distributions become larger at higher scales through scale evolution. Overall, the gluon uncertainties are likely still underestimated due to the rather rigid form of the input parameterization. For the valence quarks the resulting nPDFs are very similar as in our previous analysis, suggesting that the added EW-boson data do not provide significant constraints for the valence sector.

\begin{figure*}[tb!]
\begin{center}
\includegraphics[width=\textwidth]{Chi2_NLO-vs-NNLO.png}
\caption{Comparison of ${\large{\chi}}^2$ values divided by the individual number of data points per data set, $N_{\mathrm{dp}}$, at NLO and NNLO. The ``ideal'' value ${\large{\chi}}^2/N_{\mathrm{dp}}\,=\,1.0$ is marked by the horizontal black dotted line. The bars in the diagram corresponding to ${\large{\chi}}^2/N_{\mathrm{dp}}\,> \,3.0$ have been truncated for the purpose of a clearer representation, which is symbolised by the dashed light-grey line. The newly included data for $Z$ and $W^{\pm}$ boson production from LHC Run I and Run II are shown on the far right-hand-side.}
\label{fig-Chi2}
\end{center}
\end{figure*}

At NNLO the changes with respect to our previous analysis are clearly milder. The uncertainties have now become slightly larger for gluons and smaller for sea quarks, but otherwise these are well consistent with our previous analysis. Also for the valence quarks the differences are small, with uncertainties slightly reduced. The previously observed opposite behavior of the nuclear modifications for $u_{\mathrm{v}}$ and $d_{\mathrm{v}}$ is now less pronounced. Even though there is no significant reduction in the resulting uncertainty bands, the mutual agreement between the nuclear effects found in our NLO and NNLO fits suggests that such effects are now better captured than in the DIS-only fit.

The parameters for the input distributions for our best fit of nPDFs are also collected in Appendix~\ref{app-pdf-params}. The error sets, covering the allowed modifications of each parameter within the quoted tolerance, are included as part of the resulting LHAPDF grids. In case of nuclear PDFs, the first 33 sets reflect the central result and the uncertainties in the nuclear PDF fit and the last 26 quantify the uncertainty in the underlying proton baseline analysis.


\subsection{Comparison to data}

An overview of the resulting ${\large{\chi}}^2$ values, divided by the number of data points, $N_\mathrm{dp}$, is shown in Fig. \ref{fig-Chi2} for NLO and NNLO. Since the optimal values for the PDF parameters are obtained by minimizing ${\large{\chi}}^2$, its value is an indicator for the quality of the fit, with ${\large{\chi}}^2/N_\mathrm{dp} \approx 1$ in the optimal case. We have recalled the definition of ${\large{\chi}}^2$ used in this work in Eq. (\ref{eq-chi2-scaled}); further details can be found in Ref.~\cite{Walt:2019slu}. Values above ${\large{\chi}}^2/N_{\mathrm{dp}}>3.0$ have been truncated in Fig. \ref{fig-Chi2} for better representation, but the actual numbers are given in Table \ref{tab-expdata}.

As for the (neutrino) DIS data the fit results are very similar to those shown in our earlier analysis, we limit in the present paper the comparisons to experimental data to the new EW boson data. Our results are shown in Figs.~\ref{fig-Zboson} and \ref{fig-Wboson}. In all cases the experimental error bars are the quadratic combinations of statistical and systematical (correlated and uncorrelated) uncertainties. Figure~\ref{fig-Zboson} shows the comparisons to $Z$ boson production data measured by ATLAS \cite{Aad:2015gta} and CMS \cite{Khachatryan:2015pzs} in the p+Pb run at $\sqrt{s_{\mathrm{NN}}} = 5.02~\text{TeV}$ as a function of $Z$ boson rapidity, $y_Z$, in the center-of-mass frame, at NLO and NNLO. In each case we find a good fit, although at NLO the comparison to the ATLAS data results in ${\large{\chi}}^2/N_{\mathrm{dp}}$ slightly above unity. The Pb momentum fraction ($x_{\mathrm{Pb}}$) regions probed by the data can be estimated using LO kinematics for the process, see Eq.~(\ref{eq-dy-variables2}). To exhibit them in the context of our NNLO calculation, we plot in Fig. \ref{fig-WZboson-data} the normalized NNLO cross section as a function of $x_{\mathrm{Pb}}$, integrated over the rapidity ranges relevant for ATLAS and CMS. We notice that the ATLAS data cover a somewhat broader range in $x_{\mathrm{Pb}}$ due to broader acceptance in $y_Z$ and larger fiducial phase-space. In case of CMS the decay leptons are accepted only if $p_{\mathrm{T}}^l > 20~\text{GeV}$ and $|\eta_{\mathrm{lab}}^l| < 2.4$ and the studied mass window for the accepted dileptons is $60 < m_Z < 120~\text{GeV}$. In the ATLAS case the data have been corrected for the limited acceptance in the lepton reconstruction, so the only remaining kinematical cut is the mass window of the dilepton pair, $66 < m_Z < 116~\text{GeV}$, explaining also the larger cross section compared to the CMS data. Using Eq.~(\ref{eq-dy-variables2}) and turning back to Fig.~\ref{fig-Zboson}, we see that at forward rapidities with respect to the proton beam ($\eta > 0$, corresponding to $x_{\mathrm{Pb}} \lesssim 0.02$), both data sets clearly favor a suppression of the nPDFs relative to the proton PDFs that is well captured by the fits. In the backward direction ($\eta < 0$, $x_{\mathrm{Pb}} \gtrsim 0.02$) both data sets may have a slight tendency toward an enhancement (anti-shadowing) of p+Pb over p+p. Overall, the NNLO results are slightly higher than the NLO ones, and the shifts by the systematic uncertainties needed to obtain optimal agreement with the data are small and partly compensate for the differences between NLO and NNLO.

\begin{figure*}[bt!]
     \begin{center}                           
          \subfigure{\includegraphics[width=0.4\textwidth]{ATLAS_Z_pPb_nlo.pdf}} 
          \subfigure{\includegraphics[width=0.4\textwidth]{CMS_Z_pPb_nlo.pdf}}
          \subfigure{\includegraphics[width=0.4\textwidth]{ATLAS_Z_pPb_nnlo.pdf}}
          \subfigure{\includegraphics[width=0.4\textwidth]{CMS_Z_pPb_nnlo.pdf}}
          \end{center} 
    \caption{Comparison to $Z$-boson production data from the ATLAS \cite{Aad:2015gta} and CMS experiments \cite{Khachatryan:2015pzs} in p+Pb collisions in LHC Run~I (solid and dotted). For the solid lines, we just use our fitted nuclear PDFs to compute the cross section, while for the brown dash-dotted lines with hatched uncertainty band we also include the normalization shift of the theoretical cross section that we obtained during the fitting procedure and that resulted in the optimal ${\large{\chi}}^2$ value. For illustration, we also show the cross sections computed by using free-proton PDFs for the lead nucleus (dashed lines).}
\label{fig-Zboson}
\end{figure*}

\begin{figure*}[bt!]
     \begin{center}                           
          \subfigure{\includegraphics[width=0.4\textwidth]{CMS_Wp_pPb_nlo.pdf}} 
          \subfigure{\includegraphics[width=0.4\textwidth]{CMS_Wm_pPb_nlo.pdf}}
          \subfigure{\includegraphics[width=0.4\textwidth]{CMS_Wp_pPb_nnlo.pdf}} 
          \subfigure{\includegraphics[width=0.4\textwidth]{CMS_Wm_pPb_nnlo.pdf}}
          \end{center}
    \caption{Comparison to $W^{\pm}$-boson production data obtained by the CMS experiment \cite{CMS:2019leu} in p+Pb collisions in LHC Run~II. The lines are as in Fig.~\ref{fig-Zboson}.}
\label{fig-Wboson}    
    \end{figure*}

The CMS data for $W^{\pm}$ boson production from the p+Pb run at $\sqrt{s_{\mathrm{NN}}} = 8.16~\text{TeV}$ offer a significantly increased precision compared to the earlier data from Run~I. As the measured cross section is now given as a function of the pseudorapidity of the charged decay muon, $\eta^{\mu}$, the LO kinematics does not provide an immediate estimate for the kinematic reach of the data. Nevertheless, in Fig. \ref{fig-WZboson-data} we plot the kinematic reach in $x_{\mathrm{Pb}}$ for $W^+$ and $W^-$ production. One can see that there is sensitivity down to $x_{\mathrm{Pb}}\sim 10^{-3}$. Figure \ref{fig-Wboson} shows the comparisons of our NLO and NNLO $W^+$ and $W^-$ cross section to the ATLAS and CMS data. Again we also compute the cross sections using free-proton PDFs instead of the ones for lead. We find very clear suppression with respect to this proton baseline at $\eta^{\mu} > 0$, which is well in line with the shadowing observed in the resulting nuclear PDFs. At $\eta^{\mu} < 0$ we find only modest effects from our nuclear PDFs, which partly follows from rather modest anti-shadowing and partly since the cross section can get contributions both from regions of $x_{\mathrm{Pb}}$ where nPDF effects provide enhancement or suppression. The shapes of the resulting cross sections are well in line with the data for $W^+$ and $W^-$. This holds for both, NLO and NNLO, but in case of NNLO somewhat larger shifts ($\sim 7\%$) are required to match the data than at NLO ($\sim 3\%$), both being below or around the quoted correlated uncertainties. In all cases a reasonably good agreement, ${\large{\chi}}^2/N_{\mathrm{dp}} \sim 1.7$, is observed.
    
\begin{figure*}[bt!]
     \begin{center}                           
          \subfigure{\includegraphics[width=0.4\textwidth]{xdist_Z_pPb.pdf}} 
          \subfigure{\includegraphics[width=0.4\textwidth]{xdist_W_pPb.pdf}}
      \end{center} 
    \caption{Distribution of Pb momentum fraction ($x_{\mathrm{Pb}}$) values probed by the EW-boson production data in p+Pb collisions in 
    Run I for $Z$ bosons and Run II for $W$ bosons, obtained for our NNLO calculation.}
\label{fig-WZboson-data}    
\end{figure*}

In our previous TUJU19 analysis which included only DIS data for the proton baseline and the nuclear PDFs, we found that the resulting ${\large{\chi}}^2/N_{\mathrm{dp}}$ values were rather similar for the NLO and the NNLO analysis. Interestingly, when we add the EW boson data into the analysis, the resulting ${\large{\chi}}^2/N_{\mathrm{dp}}$ values are clearly smaller at NNLO than at NLO. For the proton baseline fit we obtain ${\large{\chi}}^2/N_{\mathrm{dp}} = 1.30$ at NLO and 1.24 at NNLO. In case of the nuclear PDF analysis, the resulting values are 0.94 (NLO) and 0.84 (NNLO). This demonstrates that even though there are still significant uncertainties in the nuclear effects in large parts of the relevant kinematical regions, it will be necessary to go beyond NLO in such analyses when more high-precision LHC data are included. It is, however, still possible that using a different parameterization or employing only normalized nuclear cross sections in the fit could result in more similar ${\large{\chi}}^2/N_{\mathrm{dp}}$ values at NLO and NNLO.

\subsection{Comparisons to other recent nPDF sets}

We next compare the nuclear PDFs obtained in this work to those of other nPDF analyses. Figure \ref{fig-nPDF-comparisonQ2-NLO-ALL} presents the nuclear modification ratios for $g,\bar{u},u_{\mathrm{v}},d_{\mathrm{v}},\mathrm{{\mathrm{v}}}$ (total valence) and $\Sigma\equiv u+\bar{u}+d+\bar{d}+s+\bar{s}$ at $Q^2 = 100~\text{GeV}^2$ for bound protons in a lead nucleus from EPPS16~\cite{Eskola:2016oht}, nCTEQ15wz~\cite{Kusina:2020lyz}, nNNPDF2.0~\cite{AbdulKhalek:2020yuc} and from this work, TUJU21. For our comparisons we stick to the more recent sets of nPDFs; similar comparisons could be carried out to earlier sets such as EPS09~\cite{Eskola:2009uj} or DSSZ~\cite{deFlorian:2011fp}. Notice that the uncertainty bands we provide in the figure do not include the uncertainty of the respective proton baselines. In the case of nNNPDF2.0 the proton PDF uncertainty would even affect the central result of the fit as the latter is obtained as the median of the ratios of replica sets. All analyses shown include some data from p+Pb collisions at the LHC, including EW boson production and also dijet production in case of EPPS16. Overall, some level of qualitative agreement can be observed within the given uncertainties, but notable differences persist even now that data from the LHC have been included. For example, a tendency for opposite valence quark modifications for $u$ and $d$ are found for nCTEQ15wz and TUJU21, whereas the wide uncertainty bands in nNNPDF2.0 and EPPS16 cover all possibilities, reflecting the fact that individual valence flavors are not well constrained by the available data as only certain combinations of them can be accessed in nuclear collisions. By contrast, the total valence nuclear modification is rather well under control for $0.01 < x < 0.5$. All the considered analyses support small-$x$ gluon shadowing, though with varying uncertainty. Also gluon anti-shadowing, enhancement somewhere around $0.05 < x < 0.5$, is present in all analyses, but the precise location in $x$ varies among the sets, with TUJU21 and nNNPDF2.0 being more in line. The largest differences arise for the $\bar{u}$ distribution. TUJU21 and nCTEQ15wz show a strong small-$x$ shadowing, whereas EPPS16 and nNNPDF2.0 have only modest suppression at $x<0.01$. The reason for these differences is likely associated with the different treatment of flavor dependence of sea quarks in the various analyses. We note that when considering the nuclear effects for the singlet distribution, $\Sigma$, all analyses are very well in line with each other, and also the uncertainties are well matched for $0.01 < x < 0.5$. 
\begin{figure*}[bt!]
     \begin{center}                           
          \subfigure{\includegraphics[width=0.32\textwidth]{pdfs_NLO_R_g_100_ALL.pdf}} 
          \subfigure{\includegraphics[width=0.32\textwidth]{pdfs_NLO_R_ubar_100_ALL.pdf}}
          \subfigure{\includegraphics[width=0.32\textwidth]{pdfs_NLO_R_Sigma_100_ALL.pdf}} 
          \subfigure{\includegraphics[width=0.32\textwidth]{pdfs_NLO_R_uv_100_ALL.pdf}}
          \subfigure{\includegraphics[width=0.32\textwidth]{pdfs_NLO_R_dv_100_ALL.pdf}}
          \subfigure{\includegraphics[width=0.32\textwidth]{pdfs_NLO_R_V_100_ALL.pdf}} 
          \end{center} 
\caption{Ratios $R_i^{\mathrm{p/Pb}}$ of PDFs in a proton bound in a lead nucleus compared to the PDFs in a free proton for TUJU21, nCTEQ15wz \cite{Kusina:2020lyz}, EPPS16 \cite{Eskola:2016oht} and nNNPDF2.0 \cite{AbdulKhalek:2020yuc}, shown at $Q^2=100\,\mathrm{GeV}^2$.}
\label{fig-nPDF-comparisonQ2-NLO-ALL}
\end{figure*}

The only other recent nPDF fit performed at NNLO (apart from nNNPDF1.0 which is now superseded by nNNPDF2.0), is the KSASG20 analysis \cite{Khanpour:2020zyu}. However, no LHC data have been included for this set. As the provided LHAPDF grids include only limited information about the uncertainties of the individual distributions, we show only the comparisons to the full-nucleus PDFs. These are provided in Fig. \ref{fig-nPDF-comparison-NNLO} for $g,s,u,d$, where $s$ represents a generic sea quark distribution since a full flavor separation for the sea quarks is not available in either of these analyses. Overall, we find reasonable agreement, although generally the gluons tend to be higher and the different quark PDFs lower in TUJU21 than in KSASG20. Notably there is some discrepancy for the sea quarks at $x\approx 0.1$, where the KSASG20 fit shows a particularly large sea-quark antishadowing.
\begin{figure*}[tb!]
     \begin{center}
          \subfigure{\includegraphics[width=0.237\textwidth]{pdfs_NNLO_KSASG20_f_g_100.pdf}}
          \subfigure{\includegraphics[width=0.237\textwidth]{pdfs_NNLO_KSASG20_f_ubar_100.pdf}}
          \subfigure{\includegraphics[width=0.237\textwidth]{pdfs_NNLO_KSASG20_f_u_100.pdf}}
          \subfigure{\includegraphics[width=0.237\textwidth]{pdfs_NNLO_KSASG20_f_d_100.pdf}}
          \end{center} 
    \caption{Nuclear parton distribution functions TUJU21 in lead at NNLO compared to the KSASG20 results \cite{Khanpour:2020zyu}, shown at $Q^2=100\,\mathrm{GeV}^2$. The comparison is presented for the distribution functions $xf_i(x,Q^2)$ with $i=g,\,s=\bar{s}=\bar{u}=\bar{d},\,u,\,d$ for a proton bound in a lead nucleus.}
\label{fig-nPDF-comparison-NNLO}    
    \end{figure*}

\section{Applications to EW-boson production in nuclear collisions}

\label{sec:applications}

\subsection{EW-boson production in Pb+Pb}

Since the $W^{\pm}$ and $Z$ bosons do not experience the strong interactions, it is expected that once they are produced in a high-energy heavy-ion collision they will be able to exit from the collision without much attenuation, even in the presence of a quark-gluon plasma that has been formed. Therefore such processes can be applied to study initial-state effects and to test the Glauber model required to normalize the measured centrality-dependent yields and to convert minimum-bias results into cross sections. No deviations from the theoretical predictions with nuclear PDFs were observed for the Run-I LHC Pb+Pb data at $\sqrt{s_{\mathrm{NN}}} = 2.76~\text{TeV}$ \cite{CMS:2011zfr, CMS:2012fgk, ATLAS:2012qdj, ATLAS:2014sic} but uncertainties in these data were fairly sizable. However, the more recent high-precision Run II data at $\sqrt{s_{\mathrm{NN}}} = 5.02~\text{TeV}$ by ATLAS \cite{ATLAS:2019maq, ATLAS:2019ibd} do show some difference in normalization when compared to NNLO calculations with NNPDF3.1 NNLO proton PDFs and EPPS16 NLO nuclear modifications \cite{Eskola:2020lee}. Similarly the recent Run-II CMS data for $Z$ boson production in Pb+Pb collisions seem to sit at the upper edge of nPDF-based predictions based on an NLO computation matched to a parton shower \cite{CMS:2021kvd}, the latter providing an approximation to leading-logarithmic resummation. Unlike ATLAS, CMS has not relied on the Glauber model to convert the measured minimum bias yield to a cross section, but utilized the measured luminosity instead. A surprising feature is that the centrality dependence of these two measurements is opposite: CMS finds a decreasing trend for the normalized yield towards more peripheral collisions, whereas in the ATLAS data it increases with centrality. The former has been explained with different possible biases in centrality classification in high-scale processes \cite{Loizides:2017sqq}, and the latter could be due to nuclear shadowing for $\sigma^{\mathrm{inel}}_{\mathrm{NN}}$ \cite{Eskola:2020lee} or anchor-point bias \cite{Jonas:2021xju}. Here we, however, confront the minimum bias data with our nPDFs constrained by the $W^{\pm}$ and $Z$ boson production in p+Pb at NLO and NNLO, in order to investigate whether these data are compatible with our nPDFs, and whether they are mutually consistent.

To calculate the EW boson production cross section at NLO and NNLO in QCD we use the \textsc{MCFM} code version 10.1 which includes an improved calculation of PDF-uncertainties from Ref.~\cite{Campbell:2019dru}. Before applying the setup to heavy-ion collisions, we validate our computation against the p+p data taken by ATLAS \cite{ATLAS:2018pyl} at the same energy, $\sqrt{s} = 5.02~\text{TeV}$, using our proton baseline PDFs that actually did not include these data sets. The comparisons for $Z$, $W^+$ and $W^-$ bosons are separately shown in Fig. \ref{fig:applicationspp} at NLO and NNLO. We find a very good agreement in all cases, which can be expected since similar data at other energies were used in our proton baseline fit.
\begin{figure*}[tb!]
     \begin{center}
          \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Z_pp_nlo.pdf}}
           \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Wp_pp_nlo.pdf}}
           \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Wm_pp_nlo.pdf}}          
      %
          \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Z_pp_nnlo.pdf}} 
           \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Wp_pp_nnlo.pdf}}
           \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Wm_pp_nnlo.pdf}}  
          \end{center} 
\caption{EW boson production cross sections in p+p collisions at $\sqrt{s} = 5.02~\text{TeV}$ for $Z$ (left), $W^+$ (center) and $W^-$ (right), calculated at NLO (upper panels) and NNLO (lower panels) and compared to data from ATLAS \cite{ATLAS:2018pyl}.}
\label{fig:applicationspp}    
\end{figure*}

In Fig.~\ref{fig:applicationsPbPb-W} we present a comparison of our NLO and NNLO calculations of $W^{\pm}$ production in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}} = 5.02~\text{TeV}$ to the measurement by ATLAS \cite{ATLAS:2019ibd}. Here we find that for both $W^+$ and $W^-$ our calculations with nuclear PDFs tend to undershoot the data, both at NLO and NNLO. Accounting for all the uncertainties, including also the normalization uncertainty, there is rough qualitative agreement, but the overall trend is that the data is consistently above the calculation. This is in line with the observation made in \cite{Eskola:2020lee}. Interestingly the calculation with only proton PDFs (but accounting for isospin effects) shows a better agreement with the data. This is unexpected when keeping in mind the strong preference for nPDF effects that we found in our analysis of p+Pb collisions. To quantify the impact of the NNLO corrections we also plot ratios between the NNLO and NLO results and compare them with the data in Fig.~\ref{fig:applicationsPbPb-W}. For this observable the NNLO corrections are only of the order few percent and do not significantly improve the agreement with the data.
\begin{figure*}[tb!]
     \begin{center}
           \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Wp_PbPb_nlo.pdf}}
          \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Wp_PbPb_nnlo.pdf}}
          \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Wp_PbPb_nlo_vs_nnlo.pdf}}
           \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Wm_PbPb_nlo.pdf}}          
           \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Wm_PbPb_nnlo.pdf}}
          \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Wm_PbPb_nlo_vs_nnlo.pdf}}
          \end{center} 
\caption{Comparison of $W^+$ (top) and $W^-$ (bottom) boson production in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}} = 5.02~\text{TeV}$ at NLO (left) and NNLO (center) with (solid with uncertainty band) and without (dashed) nuclear PDF modifications to the ATLAS data \cite{ATLAS:2019ibd}. In the right part we plot the ratios of the NNLO (red with uncertainty) and NLO (dot-dashed brown with hatched uncertainty) together with the data.}
\label{fig:applicationsPbPb-W}    
\end{figure*}

The same trend is visible in Fig.~\ref{fig:applicationsPbPb-Z} for the $Z$ boson production data in Pb+Pb collisions when comparing the NLO and NNLO calculations to the ATLAS data \cite{ATLAS:2019maq}. However, the recent CMS data for the same observable and the same collision energy \cite{CMS:2021kvd}, also shown in the figure, is well in line with our calculations with nuclear PDFs. At NNLO, it even appears that the calculated cross section is somewhat above the data, contrary to the ATLAS comparison. The differences are well visible also in the plots showing the ratio between the NNLO and NLO results together with the data in Fig.~\ref{fig:applicationsPbPb-Z}. There, the NNLO corrections grow with $y_Z$ and are of the order 10\% at the largest rapidities. The ATLAS data seem to agree with the NNLO result, whereas the CMS results seem to fall a bit below the NNLO calculation at larger rapidities and are better in line with the NLO result. Therefore our results point to possible tensions between the two data sets. That said, one should keep in mind that there are some differences in the experimental analyses: In case of ATLAS, the Glauber model was used to calculate the normalization, whereas CMS applied the measured luminosity. Also, ATLAS provides the result only in the fiducial phase-space region, while the CMS data has been corrected to include also the phase space removed by cuts on the final-state leptons. These features make direct comparisons of the two data sets difficult.
\begin{figure*}[tb!]
     \begin{center}
          \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Z_PbPb_nlo.pdf}}
          \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Z_PbPb_nnlo.pdf}}
          \subfigure{\includegraphics[width=0.325\textwidth]{ATLAS_Z_PbPb_nlo_vs_nnlo.pdf}}
          \subfigure{\includegraphics[width=0.325\textwidth]{CMS_Z_PbPb_nlo.pdf}}
          \subfigure{\includegraphics[width=0.325\textwidth]{CMS_Z_PbPb_nnlo.pdf}}
          \subfigure{\includegraphics[width=0.325\textwidth]{CMS_Z_PbPb_nlo_vs_nnlo.pdf}}
          \end{center} 
\caption{Comparison of $Z$ boson production in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}} = 5.02~\text{TeV}$ at NLO (left) and NNLO (center) with (solid with uncertainty band) and without (dashed) nuclear PDF modifications to ATLAS \cite{ATLAS:2019maq} (upper panels) and CMS \cite{CMS:2021kvd} (lower panels) data. In the right part we plot the ratios of the NNLO (red with uncertainty) and NLO (dot-dashed brown with hatched uncertainty) together with the data.}
\label{fig:applicationsPbPb-Z}
\end{figure*}


\subsection{DY production in p+Pb}
\label{sec:DYpPb}

A recent dataset that has proved difficult to include in an nPDF analysis at the NLO is the CMS DY production in p+Pb collisions \cite{CMS:2021ynu}. It has been anticipated that for the lower-mass bin ($15 < M < 60$ GeV) the NNLO corrections could be significant \cite{Khalek:2022zqe} and for the higher-mass bin ($60 < M < 120$ GeV) it has been noted that due to large fluctuations at the mid-rapidity it is difficult to have acceptable $\chi^2$ values with any PDF-based calculation \cite{Eskola:2021nhw}. Here we quantify the impact of the NNLO corrections on these data to study whether these could explain the observed differences in the low-mass bin.

Comparison with this CMS data is presented in Fig.~\ref{fig:applicationspPb-DY} for both mass windows as a function of the rapidity $y$ for the dilepton pair including also the ratios between the NNLO and NLO results. The comparisons are made for the fiducial cross section that has not been corrected for the limited acceptance. Here we notice that the NNLO corrections are rather mild, around 5\% for the high-mass bin but become significant for the low-mass bin, reaching 20\% at the largest (absolute) rapidities. To further quantify this effect we have calculated the $\chi^2/N_{\mathrm{dp}}$ values for these data at NLO and NNLO, shown in Table~\ref{tab:chi2DY} separately for the low- and high-mass bins and for the combination of these two. The common luminosity uncertainty is not included in the data uncertainties for $\chi^2$ calculation but has been accounted for by finding a common normalization factor that minimizes the combined $\chi^2/N_{\mathrm{dp}}$. In both cases the scaling factor is consistent with the quoted luminosity uncertainty of $3.5\%$. For the high-mass bin the data actually seem to be better described by the NLO calculation at negative rapidites whereas for positive rapidities it is in good agreement with the NNLO result. For the lower mass bin it seems clear that the NNLO corrections are needed to have a good agreement with this data and also the combined $\chi^2/N_{\mathrm{dp}}$ is significantly smaller at NNLO (1.554) than at NLO (2.261).
\begin{center}
\renewcommand{\arraystretch}{1.25}
\begin{table}[tb!]
\caption{The values of $\chi^2/N_{\mathrm{dp}}$ for the CMS DY data in Fig.~\ref{fig:applicationspPb-DY}. The data points have been scaled by a factor that minimizes $\chi^2$ to account for the correlated luminosity uncertainty ($3.5\%$).}
\label{tab:chi2DY}
\begin{tabular}{l|c|c}
 & $\chi^2/N_{\mathrm{dp}}$(NLO) & $\chi^2/N_{\mathrm{dp}}$(NNLO) \\ 
\hline
$15 < M < 60 $ GeV & 3.002 & 0.735 \\ 
$60 < M < 120 $ GeV & 1.894 & 2.009 \\ 
\hline
combined & 2.261 & 1.554 \\
\hline
\hline
Scaling factor & 0.989 & 1.037 \\ 
\end{tabular}

\end{table}
\end{center}

\begin{figure*}[tb!]
     \begin{center}
          \subfigure{\includegraphics[width=0.325\textwidth]{CMS_DYlow_pPb_nlo.pdf}}
          \subfigure{\includegraphics[width=0.325\textwidth]{CMS_DYlow_pPb_nnlo.pdf}}
          \subfigure{\includegraphics[width=0.325\textwidth]{CMS_DYlow_pPb_nlo_vs_nnlo.pdf}}
          \subfigure{\includegraphics[width=0.325\textwidth]{CMS_DYhigh_pPb_nlo.pdf}}
          \subfigure{\includegraphics[width=0.325\textwidth]{CMS_DYhigh_pPb_nnlo.pdf}}
          \subfigure{\includegraphics[width=0.325\textwidth]{CMS_DYhigh_pPb_nlo_vs_nnlo.pdf}}
          \end{center} 
\caption{Comparison of DY production in p+Pb collisions at $\sqrt{s_{\mathrm{NN}}} = 8.16~\text{TeV}$ at NLO (left) and NNLO (center) results with (solid with uncertainty band) and without (dashed) nuclear PDF modifications in two invariant mass bins, $15 < M < 60$ GeV (upper panels) and $60 < M < 120$ GeV (lower panels) to CMS data \cite{CMS:2021ynu}. In the right part we plot the ratios of the NNLO (red with uncertainty) and NLO (dot-dashed brown with hatched uncertainty) together with the data.}
\label{fig:applicationspPb-DY}
\end{figure*}


\section{Summary and outlook}
\label{sec:summary}

We have presented new analyses of nuclear PDFs at NLO and NNLO, TUJU21. We have adopted the same framework as in our previous TUJU19 analysis, but in addition to neutral-current DIS and charged-current neutrino DIS data we have now also included new electroweak-boson production data from the LHC, both for our proton baseline fit and for the nuclear modifications. The resulting nPDFs provide a fully consistent setup for cross section calculations at NNLO in nuclear collisions, for the first time incorporating the LHC data in a full NNLO analysis for nuclear PDFs. The comparisons to the existing nPDF sets demonstrate a reasonable agreement within the error bands, although some discrepancies in flavor dependence were observed. We do point out, however, that the adopted parametrization is rather restrictive, likely resulting in uncertainties that are underestimated in the region $x < 0.001$ where no data have been included. The resulting cross sections show very good agreement with the included experimental data, as confirmed by the total ${\large{\chi}}^2/N_{\mathrm{dp}}<1.0$ for the nuclear part of the analysis. In the presented framework, the fit performed at NNLO was found to have a significantly lower ${\large{\chi}}^2/N_{\mathrm{dp}}$ value than the NLO one, 0.84 instead of 0.94. The resulting PDFs will become available in LHAPDF6 format from the LHAPDF home page\footnote{https://lhapdf.hepforge.org/} or by request from the authors.

As an application, we have studied EW-boson production in Pb+Pb collisions at the LHC and have compared the results to recent ATLAS and CMS data. We found that for ATLAS both the NLO and the NNLO computation with the fitted nuclear PDFs tends to be below the data, even though the cross sections were well reproduced for p+Pb collisions. We find better agreement when comparing to the very recent CMS data for $Z$ boson production \cite{CMS:2021kvd}, which hints at a possible tension between the two experimental data sets. We compare our results also to the recent CMS data for DY dilepton production in p+Pb collisions \cite{CMS:2021ynu} that were not included in the presented analysis. Here we find that the NNLO corrections are significant, especially for the lower mass bin, and necessary to have a good description of the data. This demonstrates that for some observables the NNLO corrections can be larger than the uncertainties in nuclear PDF analyses and that these corrections should be taken into account when considering such data.

A possible future improvement would be to analyze the $W^{\pm}$-boson production data removing the constraint $\bar{u}=\bar{d}=\bar{s}$. However, a release of that constraint will increase the number of free parameters and have an impact on the convergence of the fit, likely requiring an extension of the analysis. Another avenue for a step forward will be a combined analysis of proton and nuclear PDFs. The ensuing doubling of the number of fit parameters will demand a full re-thinking of the minimization and PDF determination procedure.

Clearly, our results -- especially the comparisons to other sets of nPDFs -- show that we still have a long way to go until we can be confident to have a good understanding of nuclear modifications of parton distributions. Despite the still large uncertainties, we are encouraged by the improvement that the inclusion of NNLO corrections appears to provide. Ultimately, we hope that the future Electron Ion Collider (EIC) will provide precision constraints on the nuclear PDFs. On the time scale of the EIC, we also expect new analysis technologies to become available that offer new methods for extending the possibilities of theoretical investigations. Among them could be physics simulations on a quantum computer. As an example, a recent study \cite{Lamm:2019uyc} presents an algorithm for computing predictions of parton distribution functions and the hadronic tensor, where it is also speculated that ``quantum supremacy'' could possibly be demonstrated in this area of research as such a task has proven very difficult with classical computers. As another example, a recent study has presented a proof-of-principle demonstration of the determination of proton PDFs with a quantum computer \cite{Perez-Salinas:2020nem}, paving the way for further applications to, for example, protons bound in a nucleus. These new methods and new technologies will be most relevant when aiming at the highest possible perturbative precision, and we are confident that our NNLO analysis will provide a solid reference for such future studies.

\section*{Acknowledgments}

The authors acknowledge support by the state of Baden-W\"urttemberg through bwHPC providing the possibility to run the computational calculations on the high-performance cluster and also wish to acknowledge CSC – IT Center for Science, Finland, for additional computational resources through project jyy2580. Also the support from the Academy of Finland, project numbers 308301 and 331545, and through the CoE in Quark Matter, is acknowledged (I.H.). Furthermore, the authors thank HQS Quantum Simulations GmbH for the supportive working environment allowing the finalization of the publication of the presented results. This work was also supported in part by the Bundesministerium f\"{u}r Bildung und Forschung (BMBF), grant 05P18VTCA1. We thank the \textsc{xFitter} team for the support and I. Novikov for the help with the \textsc{APPLgrid} interfacing. Thanks also to H. Paukkunen for many useful discussions.

