NNLO nuclear parton distribution functions with electroweak-boson production data from the LHC

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ NNLO nuclear parton distribution functions with electroweak-boson production data from the LHC © Authors, 2022 Published version Helenius, Ilkka; Walt, Marina; Vogelsang, Werner Helenius, I., Walt, M., & Vogelsang, W. (2022). NNLO nuclear parton distribution functions with electroweak-boson production data from the LHC. Physical Review D, 105(9), Article 094031. https://doi.org/10.1103/PhysRevD.105.094031 2022 NNLO nuclear parton distribution functions with electroweak-boson production data from the LHC Ilkka Helenius ,1,2,* Marina Walt ,3,‡and Werner Vogelsang3,† 1University of Jyvaskyla, Department of Physics, P.O. Box 35, FI-40014 University of Jyvaskyla, Finland 2Helsinki Institute of Physics, P.O. Box 64, FI-00014 University of Helsinki, Finland 3Institute for Theoretical Physics, University of Tübingen, Auf der Morgenstelle 14, 72076 Tübingen, Germany (Received 14 January 2022; accepted 11 April 2022; published 24 May 2022) We present new sets of nuclear parton distribution functions (nPDFs) at next-to-leading order and next-to-next-to-leading order in perturbative QCD. Our analyses are based on deeply inelastic scattering data with charged-lepton and neutrino beams on nuclear targets, and experimental data from measurements of W;Zboson production in p þPb collisions at the LHC. In addition, a set of proton baseline PDFs is fitted within the same framework and with the same theoretical assumptions. The results of our global QCD analysis are compared to existing nPDF sets and to the previous nPDF set TUJU19, which was based on deep-inelastic scattering data only. Our work is performed using an open-source tool, X F ITTER , and the required extensions of the code are discussed as well. We find good agreement with the data included in the fit and a lower value for χ2=Ndp when performing the fit at next-to-next-to-leading order. We apply the resulting nuclear PDFs to electroweak boson production in Pb þPb collisions at the LHC and compare the results to the most recent data from ATLAS and CMS. DOI: 10.1103/PhysRevD.105.094031 I. INTRODUCTION In the framework of collinear factorization [1], 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 [2] and references therein). One may also apply the same factorization to processes involving high-energy nuclei, such as e þAorpþ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 fixedtarget 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 D-meson production at forward rapidities measured by LHCb [3] which is particularly sensitive to small-xgluon nPDFs [4], and Wand Zboson production data from ATLAS [5] and CMS [6,7], which are primarily sensitive to different quark-flavor combinations at moderate values of xand provide potential constraints for flavor dependence [8] and also for strange quarks and gluons [9]. Furthermore, dijet production [10] can offer further constraints on gluon nPDFs in the intermediate xregion [11]. Other possible observables include inclusive direct-photon and pion production in *[email protected] †werner[email protected] ‡Present address: HQS Quantum Simulations GmbH, Haidund-Neu-Straße 7, 76131 Karlsruhe, Germany. Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP3. PHYSICAL REVIEW D 105, 094031 (2022) 2470-0010=2022=105(9)=094031(19) 094031-1 Published by the American Physical Society pþPb collisions at the LHC. The former has been recently measured by ATLAS [12] and for the latter there are already several datasets available from ALICE [13–15]. 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 [16] and DSSZ [17] 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 [18] was the first to include also data from p þPb collisions at the LHC by analyzing run I data for dijets and Wand Zboson production. The original nCTEQ15 analysis [19] did not include any LHC data, but it has recently been extended to include, e.g., run II Wand Z data from p þPb collisions (nCTEQ15wz analysis [9]). Similar datasets have been considered also by the NNPDF collaboration in their most recent analysis nNNPDF2.0 [20]. The nCTEQ15 analysis has been recently updated to contain also single-inclusive hadron production [21] and also other active groups have recently prepared new analyses. The EPPS21 [22] and nNNPDF3.0 [23] analyses both include a significant amount of new LHC data including, e.g., dijet, Wand Zboson and inclusive Dmeson 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-toleading order (NNLO): nNNPDF1.0 [24], TUJU19 [25], and KSASG20 [26]. 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 neutralcurrent 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 þAcollisions 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 [27–29] which, contrary to expectations, have previously been found to be difficult to describe within a factorization-based NNLO calculation [30]. Our paper is organized as follows. We describe the theoretical framework for our analysis in Sec. II, then discuss the analysis procedure in Sec. III and the selection of experimental data in Sec. III C. The results of the analysis are presented in Sec. IV. We discuss the application to Pb þPb collisions at the LHC in Sec. V. Finally, we summarize our work in Sec. VI, where also an outlook toward future developments is given. II. THEORETICAL FRAMEWORK The theoretical framework is very similar to that adopted in our earlier analysis, TUJU19 [25], which includes also a detailed description of neutral-current (NC) and chargedcurrent (CC) DIS processes. Here we give an overview of the calculational framework that we have set up to use the new LHC data. A. Drell-Yan and W,Zboson production processes As new constraints we include LHC data for inclusive hadroproduction of electroweak (EW) bosons, Wand Z=γ. Often such processes are referred to as DY processes, originally relating to reactions where two hadrons collide and form a highly virtual photon that decays to a lepton pair [31]. The signature of such an event is a lepton pair (l¯ lwith l¼e,μ,τ) 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, q¯ q→γ→l¯ l: ð1Þ At high enough energies, similar scattering may form also other EW bosons, namely Zand W. 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¯ d initial state and the leptonic decay channel to include production of a charged lepton and a corresponding neutrino, q¯ q→Z→l¯ l; and qi¯ qj→W→lνl:ð2Þ 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 [32–34]. 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 [35,36]: d2σ dM2dy¼4πα2τ 3NCM4X q e2 q½qðx1Þ¯ qðx2Þþ¯ qðx1Þqðx2Þ;ð3Þ where NCis the number of colors, αthe fine structure constant, τ¼M2=s (with ffiffiffi s pthe center-of-mass energy), and we have introduced the variables HELENIUS, WALT, and VOGELSANG PHYS. REV. D 105, 094031 (2022) 094031-2 M2¼ðp1þp2Þ2≡ˆ s¼x1x2sðinvariant massÞ; y¼1 2logx1 x2ðrapidityÞ:ð4Þ Here p1;2are momenta of the incoming partons carrying momentum fraction x1;2of the total hadron momentum. In terms of M2and ywe have (again, to LO): x1¼M ffiffiffi s pexpðyÞand x2¼M ffiffiffi s pexpð−yÞ:ð5Þ 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 ηis used instead of yas the former does not require information of the mass of the given system. In Eq. (3),qðx1Þand ¯ qðx2Þ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¯ qð0Þ 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 →EW boson þXby ωab, we have the perturbative expansion ωab ¼ωð0Þ ab þαs πωð1Þ ab þαs π2 ωð2Þ ab þOðα3 sÞ;ð6Þ where αsis the strong coupling evaluated at some renormalization scale μR. The NLO coefficient functions ωð1Þ ab have been known for a long time [37]. The NNLO corrections ωð2Þ ab are also fully available [38–43], and their inclusion is a new feature of our TUJU21 analysis. B. Input parametrization 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 [44–47], 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 nonperturbative 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 [25] xfp iðx; Q2 0Þ¼c0xc1ð1−xÞc2ð1þc3xþc4x2Þ;ð7Þ where the index iruns over parton flavor i¼ g; dv;u v;¯ u; ¯ d; ¯ s, with the subscript “v”referring to the up or down valence-quark contribution. We have kept the flavor dependence in the valence sector, but for the sea quarks we had to assume ¯ u¼¯ d¼¯ s¼sas 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 parametrization at the charm mass threshold, Q2 0¼m2 c¼1.69 GeV2.Wedo not include any intrinsic charm content at this scale and apply the FONLL-A(C) general-mass variable-flavornumber-scheme [48,49] to calculate the heavy-quark cross sections at Q>Q 0in our NLO (NNLO) analysis. To obtain the PDFs for protons bound in a nucleus, we add dependence on the nuclear mass number Ato the parameters ciby reparametrizing them as ck→ckðAÞ¼ck;0þck;1ð1−A−ck;2Þ;ð8Þ where k¼0;…;4. A similar form has been used also in the nCTEQ15 analysis [9,19]. It is worth pointing out that with A¼1the A-dependent right-hand part of Eq. (8) vanishes and the free proton PDFs are recovered when indentifying ck;0¼ck. 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 un=A ¼dp=A and dn=A ¼up=A, and likewise for the light antiquarks. Then, an average PDF for a nucleon bound in a nucleus with Zprotons and A−Zneutrons is obtained as fN=A iðx; Q2Þ¼Z·fp=A iþðA−ZÞ·fn=A i A:ð9Þ III. ANALYSIS PROCEDURE A. 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 [25], this is done by minimizing χ2defined as χ2ðm;bÞ¼X i ½mi−Σαγi αμibα−μi2 ðδi;stat ffiffiffiffiffiffiffiffiffi μimi pÞ2þðδi;uncorrmiÞ2 þX α b2 α:ð10Þ NNLO NUCLEAR PARTON DISTRIBUTION FUNCTIONS WITH …PHYS. REV. D 105, 094031 (2022) 094031-3 Here, μiis the value of the measured data point for a given observable, whereas miis the actual theoretical value calculated using DGLAP-evolved PDFs with given parameters fckg. The uncertainties are represented by the δi;stat and δi;uncorr, which are the relative statistical and uncorrelated systematic uncertainties, respectively, as well as by the γi αμi, which are the correlated errors. Furthermore, the bαare the so-called nuisance parameters that are determined during the fitting. To account for the correlated systematic uncertainties for each dataset, we allow for shifts of the calculated cross section within the quoted uncertainty, penalizing such shifts by the additional b2 αcontributions to χ2. The minimization of χ2, Eq. (10), 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 [50,51] is used for the analysis of the uncertainties. The Hessian error analysis is performed assuming a quadratic expansion of the function χ2¼χ2 0þΔχ2around its global minimum. Here, χ2 0is the value of the function at the global minimum (with the bestfit parameters fk0g) and Δχ2is the displacement from the minimum [50,51]. During the performed error analysis Δχ2 defines the tolerance criterion determining the allowed growth of χ2. As argued in our previous work [25], for the proton baseline with 13 free fit parameters we select Δχ2¼20 and for the nuclear PDF error analysis we choose Δχ2¼50 for our 16 free parameters. FIG. 1. Schematic view of the high-level X F ITTER functionalities as relevant for Drell-Yan and W,Zboson production processes, with the newly implemented convolution routine in X F ITTER as used for the nPDF fit TUJU21. X F ITTER logo credited from [60]. HELENIUS, WALT, and VOGELSANG PHYS. REV. D 105, 094031 (2022) 094031-4 B. The fitting framework Our global analyses of the baseline proton and nuclear PDFs are performed with the X F ITTER [52,53] tool. The main goal of the X F ITTER project is to provide an opensource 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. [25]. 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 APPL GRID [54], can be linked with X F ITTER , facilitating fast cross section calculations using the precomputed 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 X F ITTER requires grids in a ROOT format [55]. It is also possible to use Kfactors when going from NLO to NNLO cross sections to reduce the computing effort. We have used such Kfactors in our proton baseline fit in cases where they were available in X F ITTER for a given dataset. The settings and PDFs used to calculate the Kfactors can be found from the references provided in Table II. The schematic overview of the fitting routine and the required tool set ( X F ITTER and all other modules) as applied for the nuclear PDF fit is shown in Fig. 1. For this, a new convolution routine had to be implemented in X F ITTER 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 [56]. In this work, MCFM8.3 [57–59] 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 MCFM up to order NNLO, without relying on approximative Kfactors. During the fitting procedure the actual theoretical predictions were obtained by convoluting the fitted nPDFs based on updated TABLE I. Summary of experimental DIS data used to determine our proton PDF baseline. In the last two columns the χ2 values at NLO and NNLO obtained in our analysis are provided. Exp. Dataset Year Ref. Ndp χ2NLO χ2NNLO BCDMS F2p 100 GeV 1996 [65] 83 96.30 93.69 F2p 120 GeV 90 70.54 68.70 F2p 200 GeV 79 91.81 86.32 F2p 280 GeV 75 67.52 69.71 HERA 1þ2NCep 920 2015 [66] 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 NMC-97 NCep 1997 [67] 100 117.72 111.31 In total: 1559 TABLE II. Summary of experimental DYand W;Zboson production data used to determine our proton PDFs. In the last two columns the χ2values at NLO and NNLO obtained in our analysis are provided. Exp. Dataset Year Ref. Ndp χ2NLO χ2NNLO ATLAS DY high mass DY 2013 [61] 13 10.91 11.54 low mass DY 2014 [62] 8 22.86 8.68 ATLAS W;Z Wþlepton η2012 [63] 11 15.77 14.00 W−lepton η11 7.98 8.54 Zy 8 4.10 2.71 Wþlepton y2016 [64] 11 19.57 10.71 W−lepton y11 11.11 11.82 high mass CC Zy 6 7.61 6.05 high mass CF Zy 6 3.92 3.95 low mass Zy 6 41.32 23.84 peak CC Zy 12 45.76 14.40 peak CF Zy 9 21.85 7.55 CMS WWþσ8TeV 2016 [68] 11 10.64 4.81 W−σ8 TeV 11 8.14 9.27 In total: 134 NNLO NUCLEAR PARTON DISTRIBUTION FUNCTIONS WITH …PHYS. REV. D 105, 094031 (2022) 094031-5 parameters with the precalculated 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. C. Experimental data 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 I, also showing the resulting χ2 values at NLO and NNLO obtained in the analysis. The newly included experimental data from DY and Wand Zboson production for the proton baseline are listed in Table II. For the experimental proton data used, the fast interpolation grids were publicly available and the details of grid generation for each proton PDF dataset can be found from the references provided in Table I. These data include highand low-mass DY data from ATLAS at ffiffiffi s p¼7TeV [61,62], EW boson production data from ATLAS at ffiffiffi s p¼7TeV [63], and increased-luminosity data for these from ATLAS at ffiffiffi s p¼7TeV [64]. In addition, also W production cross sections from CMS at ffiffiffi s p¼8TeV have been included. In total the newly included data sets consist of 134 data points. Table III provides a list of nuclear-DIS data as also used in the TUJU19 analysis, but with the χ2values obtained in this analysis. The input data files and the fast interpolation grids used for the fitting procedure laid out in Fig. 1for nPDFs have been collected and prepared as part of this analysis for the newly included data points summarized in Table IV. The added data include run I measurements of Z boson production in p þPb collisions at the LHC by ATLAS [5] and CMS [7] at ffiffiffiffiffiffiffiffi sNN p¼5.02 TeV and the more recent run II measurement of Wboson production in pþPb collisions at ffiffiffiffiffiffiffiffi sNN p¼8.16 TeV by CMS [69].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 [70] and LHCb [71] 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production from CMS [72] 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 [73], E866 [74], and CLAS [75] 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 Sec. VB we consider the recent run II CMS measurement for Z=γproduction for which it has proven difficult to obtain χ2=Ndp values close to unity in a NLO QCD analysis [22,23]. As we have discussed in the context of the TUJU19 nPDF analysis, some authors have found tension between the neutrinoand charged-lepton DIS data. To check for potential tension with the new LHC data we have also performed fits without any neutrino-DIS TABLE III. Summary of experimental DIS data used to determine the nuclear PDFs. In the last two columns the χ2 values at NLO and NNLO obtained in our analysis are provided. Nucleus Exp. Year Ref. Ndp χ2NLO χ2NNLO D NMC 97 1996 [67] 120 151.61 121.52 EMC 90 1989 [76] 21 24.31 22.89 He=D HERMES 2002 [77] 7 6.79 8.92 NMC 95, re. 1995 [78] 13 10.67 10.42 SLAC E139 1994 [79] 11 6.47 4.42 Li=D NMC 95 1995 [80] 12 9.10 9.00 Be=D SLAC E139 1994 [79] 10 11.58 11.51 Be=C NMC 96 1996 [81] 14 13.56 16.06 C EMC 90 1989 [76] 17 13.41 13.44 C=D FNAL E665 1995 [82] 3 2.00 1.83 SLAC E139 1994 [79] 6 20.69 13.86 EMC 88 1988 [83] 9 3.70 4.22 NMC 95, re. 1995 [78] 13 34.96 19.49 C=Li NMC 95, re. 1995 [78] 10 7.77 10.18 N=D HERMES 2002 [77] 1 0.95 2.08 Al=D SLAC E139 1994 [79] 10 18.49 9.49 Al=C NMC 96 1996 [81] 14 7.29 6.23 Ca EMC 90 1989 [76] 19 13.41 13.44 Ca=D NMC 95, re. 1995 [78] 12 34.75 21.82 FNAL E665 1995 [82] 3 1.84 2.41 SLAC E139 1994 [79] 6 16.74 9.01 Ca=Li NMC 95, re. 1995 [78] 10 1.45 1.33 Ca=C NMC 95, re. 1995 [78] 10 9.35 8.00 NMC 96 1996 [81] 14 8.45 6.42 Fe SLAC E140 1993 [84] 2 0.14 0.04 Fe=D SLAC E139 1994 [79] 14 44.53 32.07 Fe=C NMC 96 1996 [81] 14 11.17 9.93 νFe CDHSW 1991 [85] 464 404.26 358.19 ¯ νFe CDHSW 1991 [85] 462 439.53 395.99 Cu=D EMC 88 1988 [83] 9 8.38 5.71 EMC 93 1993 [86] 19 26.38 12.58 Kr=D HERMES 2002 [77] 1 2.02 2.02 Ag=D SLAC E139 1994 [79] 6 21.37 18.80 Sn=D EMC 88 1988 [83] 8 14.37 13.98 Sn=C NMC 96 1996 [81] 14 6.48 8.52 NMC 96, Q2dep. 1996 [87] 134 76.13 75.03 Xe=D FNAL E665 1992 [88] 3 1.64 1.34 Au=D SLAC E139 1994 [79] 11 16.89 18.66 Pb=D FNAL E665 1995 [82] 2 8.28 7.72 Pb=C NMC 96 1996 [81] 14 8.32 5.42 νPb CHORUS 2005 [89] 405 259.48 237.85 ¯ νPb CHORUS 2005 [89] 405 356.01 352.09 In total: 2336 HELENIUS, WALT, and VOGELSANG PHYS. REV. D 105, 094031 (2022) 094031-6 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 datasets. IV. RESULTS A. Proton baseline 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;Z boson production processes taken by the ATLAS and CMS collaborations at the LHC; see Sec. III C for details. The comparisons are presented in Fig. 2for NLO and in Fig. 3for 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;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 the Appendix. B. Nuclear PDFs The resulting nuclear PDFs, referred to as TUJU21, are presented in Figs. 4at NLO and 5at NNLO for the boundproton-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-xsuppression is significantly milder than in TUJU19, along with a TABLE IV. Summary of experimental W;Zboson 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 χ2values at NLO and NNLO obtained in our analysis are provided. Nucleus Proc. Exp. Year Ref. Ndp χ2NLO χ2NNLO pPb ZLHC run I ATLAS 2015 [5] 14 16.40 12.82 ZLHC run I CMS 2015 [7] 12 8.76 7.30 W−LHC run II CMS 2019 [69] 24 39.58 42.83 WþLHC run II CMS 24 41.08 39.07 In total: 74 FIG. 2. Proton baseline PDFs in TUJU21 at NLO compared to the previous TUJU19 results, shown at the initial scale Q2 0¼ 1.69 GeV2(upper panels) and at Q2¼100 GeV2(lower panels) after DGLAP evolution. NNLO NUCLEAR PARTON DISTRIBUTION FUNCTIONS WITH …PHYS. REV. D 105, 094031 (2022) 094031-7 slightly reduced uncertainty. Also, the strong antishadowing enhancement at intermediate values of xwe found previously is now tamed to a more moderate ∼10% effect. At the initial scale of the fit the sea quark nPDFs are now slightly lower in the small-xregion, 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 parametrization. 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. 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 uvand dvis 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 the Appendix. The error sets, FIG. 3. Same as for Fig. 2, but at NNLO. FIG. 4. NLO nuclear parton distribution functions in TUJU21 for a lead nucleus, compared to the previous TUJU19 results, shown at the initial scale Q2 0¼1.69 GeV2(upper panels) and at Q2¼100 GeV2after DGLAP evolution (center panels). The lower panels show the corresponding ratios of PDFs for a proton bound in lead over the free proton PDFs. HELENIUS, WALT, and VOGELSANG PHYS. REV. D 105, 094031 (2022) 094031-8 FIG. 14. Comparison of Zboson production in Pb þPb collisions at ffiffiffiffiffiffiffiffi sNN p¼5.02 TeV at NLO (left) and NNLO (center) with (solid with uncertainty band) and without (dashed) nuclear PDF modifications to ATLAS [28] (upper panels) and CMS [29] (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. FIG. 15. Comparison of DY production in p þPb collisions at ffiffiffiffiffiffiffiffi sNN p¼8.16 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 [98]. 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. NNLO NUCLEAR PARTON DISTRIBUTION FUNCTIONS WITH …PHYS. REV. D 105, 094031 (2022) 094031-15 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 datasets. 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 datasets difficult. B. DY production in p + Pb 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 [98]. It has been anticipated that for the lower-mass bin (15 <M<60 GeV) the NNLO corrections could be significant [23] and for the higher-mass bin (60 <M<120 GeV) it has been noted that due to large fluctuations at the midrapidity it is difficult to have acceptable χ2values with any PDF-based calculation [22]. 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. 15 for both mass windows as a function of the rapidity yfor 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 χ2=Ndp values for these data at NLO and NNLO, shown in Table Vseparately for the lowand high-mass bins and for the combination of these two. The common luminosity uncertainty is not included in the data uncertainties for χ2calculation but has been accounted for by finding a common normalization factor that minimizes the combined χ2=Ndp. 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 rapidities 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 χ2=Ndp is significantly smaller at NNLO (1.554) than at NLO (2.261). VI. SUMMARY AND OUTLOOK 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 electroweakboson 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 χ2=Ndp <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 χ2=Ndp value than the NLO one, 0.84 instead of 0.94. The resulting PDFs will become available in LHAPDF6 format from the LHAPDF home page1or 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 Zboson production [29], which hints at a possible tension between the two experimental datasets. We compare our results also to the recent CMS data for DY dilepton production in p þPb collisions [98] 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 TABLE V. The values of χ2=Ndp for the CMS DY data in Fig. 15. The data points have been scaled by a factor that minimizes χ2to account for the correlated luminosity uncertainty (3.5%). χ2=NdpðNLOÞχ2=NdpðNNLOÞ 15 <M<60 GeV 3.002 0.735 60 <M<120 GeV 1.894 2.009 Combined 2.261 1.554 Scaling factor 0.989 1.037 1https://lhapdf.hepforge.org/. HELENIUS, WALT, and VOGELSANG PHYS. REV. D 105, 094031 (2022) 094031-16 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-boson production data removing the constraint ¯ u¼¯ d¼¯ 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 rethinking 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 will provide precision constraints on the nuclear PDFs. On the time scale of the electron ion collider, 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 [99] 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 [100], paving the way for further applications such as, 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. ACKNOWLEDGMENTS The authors acknowledge support by the state of Baden-Württemberg through bwHPC providing the possibility to run the computational calculations on the highperformance cluster and also wish to acknowledge CSC-IT Center for Science, Finland, for additional computational resources through Project No. jyy2580. Also the support from the Academy of Finland, Projects No. 308301 and No. 331545, and through the Center of Excellence 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ür Bildung und Forschung (BMBF), Grant No. 05P18VTCA1. We thank the X F ITTER team for the support and I. Novikov for the help with the APPL GRID interfacing. Thanks also to H. Paukkunen for many useful discussions. APPENDIX: PDF PARAMETERS Here we collect the input parameters obtained for the proton and nuclear parton distribution functions presented in Sec. IV. The naming convention corresponds to the PDF parametrization given in Eqs. (7) and (8). Table VI provides the NLO parameters, while Table VII presents the NNLO ones. Some of the parameters were manually set to zero if the data that were used did not provide enough sensitivity to constrain them without large uncertainty. The Adependence was implemented for a subset of parameters, again selected such that the data provided enough sensitivity to result in a converged fit. TABLE VI. Values of the NLO fit parameters at the initial scale, Q2 0¼1.69 GeV2. (SR) means that the normalization for that particular parton is fixed by the momentum and valence number sum rules. A dash indicates that this parameter was excluded from the fit. Parameter values for the sea quarks, apart from ¯ u, were derived from the applied constraints ¯ s¼s¼¯ d¼¯ u. gValue uvValue dvValue ¯ uValue cg 0;08.9596 cuv 0;0(SR) cdv 0;0(SR) c¯ u 0;0(SR) cg 1;00.3270 cuv 1;00.7121 cdv 1;00.7629 c¯ u 1;0−0.1815 cg 2;013.438 cuv 2;03.4290 cdv 2;02.0996 c¯ u 2;05.2593 cg 3;06.4371 cuv 3;01.4506 cdv 3;0−1.4391 c¯ u 3;02.4151 cg 4;0 cuv 4;0 cdv 4;0 c¯ u 4;0 cg 1;1−5.4728 cuv 1;1−0.0462 cdv 1;1−19.16 c¯ u 1;1251.91 cg 1;2−0.0013 cuv 1;20.3411 cdv 1;2−0.0026 c¯ u 1;20.0002 cg 2;1−2.000 cuv 2;14.2325 cdv 2;11.2264 c¯ u 2;1−276.53 cg 2;20.3695 cuv 2;20.0025 cdv 2;20.4273 c¯ u 2;2−0.0017 NNLO NUCLEAR PARTON DISTRIBUTION FUNCTIONS WITH …PHYS. REV. D 105, 094031 (2022) 094031-17 [1] J. C. Collins, D. E. Soper, and G. F. Sterman, Adv. Ser. Dir. High Energy Phys. 5, 1 (1989). [2] J. Gao, L. Harland-Lang, and J. Rojo, Phys. Rep. 742,1 (2018). [3] R. Aaij et al. (LHCb Collaboration), J. High Energy Phys. 10 (2017) 090. [4] K. J. Eskola, I. Helenius, P. Paakkinen, and H. Paukkunen, J. High Energy Phys. 05 (2020) 037. [5] G. Aad et al. (ATLAS Collaboration), Phys. Rev. C 92, 044915 (2015). [6] V. Khachatryan et al. (CMS Collaboration), Phys. Lett. B 750, 565 (2015). [7] V. Khachatryan et al. (CMS Collaboration), Phys. Lett. B 759, 36 (2016). [8] H. Paukkunen and C. A. Salgado, J. High Energy Phys. 03 (2011) 071. [9] A. Kusina et al.,Eur. Phys. J. C 80, 968 (2020). [10] A. M. Sirunyan et al. (CMS Collaboration), Phys. Rev. Lett. 121, 062002 (2018). [11] K. J. Eskola, P. Paakkinen, and H. Paukkunen, Eur. Phys. J. C 79, 511 (2019). [12] M. Aaboud et al. (ATLAS Collaboration), Phys. Lett. B 796, 230 (2019). [13] S. Acharya et al. (ALICE Collaboration), Eur. Phys. J. C 78, 624 (2018). [14] J. Adam et al. (ALICE Collaboration), Phys. Lett. B 760, 720 (2016). [15] S. Acharya et al. (ALICE Collaboration), Phys. Lett. B 827, 136943 (2022). [16] K. J. Eskola, H. Paukkunen, and C. A. Salgado, J. High Energy Phys. 04 (2009) 065. [17] D. de Florian, R. Sassot, P. Zurita, and M. Stratmann, Phys. Rev. D 85, 074028 (2012). [18] K. J. Eskola, P. Paakkinen, H. Paukkunen, and C. A. Salgado, Eur. Phys. J. C 77, 163 (2017). [19] K. Kovarik et al.,Phys. Rev. D 93, 085037 (2016). [20] R. Abdul Khalek, J. J. Ethier, J. Rojo, and G. van Weelden, J. High Energy Phys. 09 (2020) 183. [21] P. Duwentäster, L. A. Husová, T. Ježo, M. Klasen, K.Kovaˇ rík, A. Kusina, K. F. Muzakka, F. I. Olness, I. Schienbein, and J. Y. Yu, Phys.Rev.D104, 094005 (2021). [22] K. J. Eskola, P. Paakkinen, H. Paukkunen, and C. A. Salgado, arXiv:2112.12462. [23] R. A. Khalek, R. Gauld, T. Giani, E. R. Nocera, T. R. Rabemananjara, and J. Rojo, arXiv:2201.12363. [24] R. Abdul Khalek, J. J. Ethier, and J. Rojo (NNPDF Collaboration), Eur. Phys. J. C 79, 471 (2019). [25] M. Walt, I. Helenius, and W. Vogelsang, Phys. Rev. D 100, 096015 (2019). [26] H. Khanpour, M. Soleymaninia, S. Atashbar Tehrani, H. Spiesberger, and V. Guzey, Phys. Rev. D 104, 034010 (2021). [27] G. Aad et al. (ATLAS Collaboration), Eur. Phys. J. C 79, 935 (2019). [28] G. Aad et al. (ATLAS Collaboration), Phys. Lett. B 802, 135262 (2020). [29] A. M. Sirunyan et al. (CMS Collaboration), Phys. Rev. Lett. 127, 102002 (2021). [30] K. J. Eskola, I. Helenius, M. Kuha, and H. Paukkunen, Phys. Rev. Lett. 125, 212301 (2020). [31] S. D. Drell and T.-M. Yan, Phys. Rev. Lett. 25, 316 (1970); 25, 902(E) (1970). [32] G. T. Bodwin, Phys. Rev. D 31, 2616 (1985);34, 3932(E) (1986). [33] J. C. Collins, D. E. Soper, and G. F. Sterman, Nucl. Phys. B261, 104 (1985). [34] J. C. Collins, D. E. Soper, and G. F. Sterman, Nucl. Phys. B308, 833 (1988). [35] R. K. Ellis, W. J. Stirling, and B. R. Webber, QCD and Collider Physics (Cambridge University Press, Cambridge, England, 1996). [36] G. Sterman et al.,Handbook of Perturbative QCD (CTEQ, College Park, Maryland, 2001). [37] G. Altarelli, R. K. Ellis, and G. Martinelli, Nucl. Phys. B157, 461 (1979). [38] R. Hamberg, W. L. van Neerven, and T. Matsuura, Nucl. Phys. B359, 343 (1991);B644, 403(E) (2002). TABLE VII. Same as Table VI, but at NNLO. gValue uvValue dvValue ¯ uValue cg 0;06.4747 cuv 0;0(SR) cdv 0;0(SR) c¯ u 0;0(SR) cg 1;00.2858 cuv 1;00.7157 cdv 1;00.9101 c¯ u 1;0−0.1197 cg 2;07.6890 cuv 2;03.6964 cdv 2;03.8936 c¯u 2;08.0188 cg 3;0−0.0413 cuv 3;02.5811 cdv 3;0−0.5844 c¯u 3;0 cg 4;0 cuv 4;0 cdv 4;0 c¯u 4;011.960 cg 1;12.9882 cuv 1;1−0.0235 cdv 1;1−0.6681 c¯ u 1;1−85.228 cg 1;20.0003 cuv 1;20.6564 cdv 1;2−0.0376 c¯ u 1;2−0.0005 cg 2;1−0.6166 cuv 2;115.614 cdv 2;11.2905 c¯ u 2;1−0.1323 cg 2;20.4518 cuv 2;2−0.0011 cdv 2;20.3396 c¯ u 2;2−0.4051 HELENIUS, WALT, and VOGELSANG PHYS. REV. D 105, 094031 (2022) 094031-18 [39] R. V. Harlander and W. B. Kilgore, Phys. Rev. Lett. 88, 201801 (2002). [40] C. Anastasiou, L. J. Dixon, K. Melnikov, and F. Petriello, Phys. Rev. Lett. 91, 182002 (2003). [41] C. Anastasiou, L. J. Dixon, K. Melnikov, and F. Petriello, Phys. Rev. D 69, 094008 (2004). [42] S. Catani, L. Cieri, G. Ferrera, D. de Florian, and M. Grazzini, Phys. Rev. Lett. 103, 082001 (2009). [43] R. Gavin, Y. Li, F. Petriello, and S. Quackenbush, Comput. Phys. Commun. 182, 2388 (2011). [44] V. N. Gribov and L. N. Lipatov, Yad. Fiz. 15, 781 (1972) [Sov. J. Nucl. Phys. 15, 438 (1972)]. [45] L. N. Lipatov, Yad. Fiz. 20, 181 (1974) [Sov. J. Nucl. Phys. 20, 94 (1975)]. [46] G. Altarelli and G. Parisi, Nucl. Phys. B126, 298 (1977). [47] Y. L. Dokshitzer, Zh. Eksp. Teor. Fiz. 73, 1216 (1977) [Sov. Phys. JETP 46, 641 (1977)]. [48] M. Cacciari, M. Greco, and P. Nason, J. High Energy Phys. 05 (1998) 007. [49] S. Forte, E. Laenen, P. Nason, and J. Rojo, Nucl. Phys. B834, 116 (2010). [50] J. Pumplin, D. Stump, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai, and W. K. Tung, Phys. Rev. D 65, 014013 (2001). [51] J. Pumplin, D. R. Stump, and W. K. Tung, Phys. Rev. D 65, 014011 (2001). [52] O. Zenaiev ( X F ITTER team Collaboration), Proc. Sci., DIS2016 (2016) 033. [53] V. Bertone et al. ( X F ITTER Developers’Team), Proc. Sci., DIS2017 203 (2018). [54] T. Carli, D. Clements, A. Cooper-Sarkar, C. Gwenlan, G. P. Salam, F. Siegert, P. Starovoitov, and M. Sutton, Eur. Phys. J. C 66, 503 (2010). [55] R. Brun and F. Rademakers, Nucl. Instrum. Methods Phys. Res., Sect. A 389, 81 (1997). [56] A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr, and G. Watt, Eur. Phys. J. C 75, 132 (2015). [57] A. Falkowski, M. L. Mangano, A. Martin, G. Perez, and J. Winter, Phys. Rev. D 87, 034039 (2013). [58] J. M. Campbell, R. K. Ellis, and W. T. Giele, Eur. Phys. J. C75, 246 (2015). [59] R. Boughezal, J. M. Campbell, R. K. Ellis, C. Focke, W. Giele, X. Liu, F. Petriello, and C. Williams, Eur. Phys. J. C 77, 7 (2017). [60] X F ITTER :https://www.xfitter.org/xfitter. [61] G. Aad et al. (ATLAS Collaboration), Phys. Lett. B 725, 223 (2013). [62] G. Aad et al. (ATLAS Collaboration), J. High Energy Phys. 06 (2014) 112. [63] G. Aad et al. (ATLAS Collaboration), Phys. Rev. Lett. 109, 012001 (2012). [64] M. Aaboud et al. (ATLAS Collaboration), Eur. Phys. J. C 77, 367 (2017). [65] A. C. Benvenuti et al. (BCDMS Collaboration), Phys. Lett. B223, 485 (1989). [66] H. Abramowicz et al. (H1, ZEUS Collaborations), Eur. Phys. J. C 75, 580 (2015). [67] M. Arneodo et al. (New Muon Collaboration), Nucl. Phys. B483, 3 (1997). [68] V. Khachatryan et al. (CMS Collaboration), Eur. Phys. J. C 76, 469 (2016). [69] A. M. Sirunyan et al. (CMS Collaboration), Phys. Lett. B 800, 135048 (2020). [70] J. Adam et al. (ALICE Collaboration), J. High Energy Phys. 02 (2017) 077. [71] R. Aaij et al. (LHCb Collaboration), J. High Energy Phys. 09 (2014) 030. [72] V. Khachatryan et al. (CMS Collaboration), Phys. Lett. B 750, 565 (2015). [73] D. M. Alde et al.,Phys. Rev. Lett. 64, 2479 (1990). [74] M. A. Vasilev et al. (NuSea Collaboration), Phys. Rev. Lett. 83, 2304 (1999). [75] B. Schmookler et al. (CLAS Collaboration), Nature (London) 566, 354 (2019). [76] M. Arneodo et al. (European Muon Collaboration), Nucl. Phys. B333, 1 (1990). [77] A. Airapetian et al. (HERMES Collaboration), arXiv:hepex/0210068. [78] P. Amaudruz et al. (New Muon Collaboration), Nucl. Phys. B441, 3 (1995). [79] J. Gomez et al.,Phys. Rev. D 49, 4348 (1994). [80] M. Arneodo et al. (New Muon Collaboration), Nucl. Phys. B441, 12 (1995). [81] M. Arneodo et al. (New Muon Collaboration), Nucl. Phys. B481, 3 (1996). [82] M. R. Adams et al. (E665 Collaboration), Z. Phys. C 67, 403 (1995). [83] J. Ashman et al. (European Muon Collaboration), Phys. Lett. B 202, 603 (1988). [84] S. Dasu et al.,Phys. Rev. D 49, 5641 (1994). [85] J. P. Berge et al.,Z. Phys. C 49, 187 (1991). [86] J. Ashman et al. (European Muon Collaboration), Z. Phys. C57, 211 (1993). [87] M. Arneodo et al. (New Muon Collaboration), Nucl. Phys. B481, 23 (1996). [88] M. R. Adams et al. (E665 Collaboration), Phys. Rev. Lett. 68, 3266 (1992). [89] G. Onengut et al. (CHORUS Collaboration), Phys. Lett. B 632, 65 (2006). [90] S. Chatrchyan et al. (CMS Collaboration), Phys. Rev. Lett. 106, 212301 (2011). [91] S. Chatrchyan et al. (CMS Collaboration), Phys. Lett. B 715, 66 (2012). [92] G. Aad et al. (ATLAS Collaboration), Phys. Rev. Lett. 110, 022301 (2013). [93] G. Aad et al. (ATLAS Collaboration), Eur. Phys. J. C 75, 23 (2015). [94] C. Loizides and A. Morsch, Phys. Lett. B 773, 408 (2017). [95] F. Jonas and C. Loizides, Phys. Rev. C 104, 044905 (2021). [96] J. Campbell and T. Neumann, J. High Energy Phys. 12 (2019) 034. [97] M. Aaboud et al. (ATLAS Collaboration), Eur. Phys. J. C 79, 128 (2019);79, 374(E) (2019). [98] A. M. Sirunyan et al. (CMS Collaboration), J. High Energy Phys. 05 (2021) 182. [99] H. Lamm, S. Lawrence, and Y. Yamauchi (NuQS Collaboration), Phys. Rev. Research 2, 013272 (2020). [100] A. P´erez-Salinas, J. Cruz-Martinez, A. A. Alhajri, and S. Carrazza, Phys. Rev. D 103, 034027 (2021). NNLO NUCLEAR PARTON DISTRIBUTION FUNCTIONS WITH …PHYS. REV. D 105, 094031 (2022) 094031-19