# Impact of heavy quark and quarkonium data on nuclear gluon PDFs

P. Duwentäster<EMAIL_ADDRESS>Institut für Theoretische Physik,
Westfälische Wilhelms-Universität Münster, Wilhelm-Klemm-Straße 9, D-48149
Münster, Germany  T. Ježo<EMAIL_ADDRESS>Institut für
Theoretische Physik, Westfälische Wilhelms-Universität Münster, Wilhelm-Klemm-
Straße 9, D-48149 Münster, Germany  M. Klasen  Institut für Theoretische
Physik, Westfälische Wilhelms-Universität Münster, Wilhelm-Klemm-Straße 9,
D-48149 Münster, Germany  K. Kovařík  Institut für Theoretische Physik,
Westfälische Wilhelms-Universität Münster, Wilhelm-Klemm-Straße 9, D-48149
Münster, Germany  A. Kusina  Institute of Nuclear Physics Polish Academy of
Sciences, PL-31342 Krakow, Poland  K. F. Muzakka  Institut für Theoretische
Physik, Westfälische Wilhelms-Universität Münster, Wilhelm-Klemm-Straße 9,
D-48149 Münster, Germany  F. I. Olness  Southern Methodist University, Dallas,
TX 75275, USA  R. Ruiz  Institute of Nuclear Physics Polish Academy of
Sciences, PL-31342 Krakow, Poland  I. Schienbein  Laboratoire de Physique
Subatomique et de Cosmologie, Université Grenoble-Alpes, CNRS/IN2P3, 53 avenue
des Martyrs, 38026 Grenoble, France  J. Y. Yu  Laboratoire de Physique
Subatomique et de Cosmologie, Université Grenoble-Alpes, CNRS/IN2P3, 53 avenue
des Martyrs, 38026 Grenoble, France

###### Abstract

A clear understanding of nuclear parton distribution functions (nPDFs) plays a
crucial role in the interpretation of collider data taken at the Relativistic
Heavy Ion Collider (RHIC), the Large Hadron Collider (LHC) and in the near
future at the Electron-Ion Collider (EIC). Even with the recent inclusions of
vector boson and light meson production data, the uncertainty of the gluon PDF
remains substantial and limits the interpretation of heavy ion collision data.
To obtain new constraints on the nuclear gluon PDF, we extend our recent
nCTEQ15WZ+SIH analysis to inclusive quarkonium and open heavy-flavor meson
production data from the LHC. This vast new data set covers a wide kinematic
range and puts strong constraints on the nuclear gluon PDF down to $x\lesssim
10^{-5}.$ The theoretical predictions for these data sets are obtained from a
data-driven approach, where proton-proton data are used to determine effective
scattering matrix elements. This approach is validated with detailed
comparisons to existing next-to-leading order (NLO) calculations in non-
relativistic QCD (NRQCD) for quarkonia and in the general-mass variable-
flavor-number scheme (GMVFNS) for the open heavy-flavored mesons. In addition,
the uncertainties from the data-driven approach are determined using the
Hessian method and accounted for in the PDF fits. This extension of our
previous analyses represents an important step toward the next generation of
PDFs not only by including new data sets, but also by exploring new methods
for future analyses.

††preprint:  MS-TP-22-10, IFJPAN-IV-2022-5

###### Contents

  1. I Introduction
  2. II Theoretical approach
    1. II.1 The nCTEQ framework
    2. II.2 The data-driven approach
  3. III Proton-Proton baseline
    1. III.1 Cuts and excluded data
    2. III.2 Baseline fit
    3. III.3 Comparison with $J/\psi$ production in NRQCD
    4. III.4 Comparison with $D^{0}$ production in the GMVFNS
  4. IV Impact of HQ data on nPDF fits
    1. IV.1 Data selection
    2. IV.2 Resulting PDFs
    3. IV.3 Fit quality
  5. V Conclusions
  6. A Data comparisons

##  I Introduction

Parton distribution functions (PDFs) are fundamental quantities required to
calculate predictions for processes involving hadronic initial states. The
underlying theoretical framework is based on factorization theorems, which
have been proven from first principles of QCD for a number of collider
processes in $ep$ and $pp$ collisions Collins:1989gx . This formalism provides
both a field theoretical definition of the PDFs, and the definition of the
short distance hard scattering cross sections at the partonic level.
Additionally, it includes a statement about the error of the collinear
factorization formula, which is inversely proportional to some power of the
hard scale of the process. The predictive power of this formalism lies in the
fact that the PDFs are universal, i.e. process independent, whereas the
process dependent short distance cross sections can be systematically
calculated in perturbation theory. This approach has been widely used in
analyses of proton PDFs, which have been constrained with great precision
Hou:2019efy ; Ball:2017nwa ; Khalek:2018mdn ; Gao:2017yyd ; Kovarik:2019xvh ;
Alekhin:2017olj ; Nadolsky:2008zw ; Sato:2019yez ; Harland-Lang:2014zoa ;
Thorne:2019mpt ; Ball:2009mk ; Lin:2017snn ; Lin:2020rut . Assuming that the
twist-2 collinear factorization remains valid also in the case of $eA$ and
$pA$ collisions, nuclear parton distribution functions (nPDFs) have been
determined Kovarik:2015cma ; Eskola:2016oht ; AbdulKhalek:2019mzd ;
AbdulKhalek:2020yuc ; Ethier:2020way ; Guzey:2019kik ; Klasen:2017kwb ;
Klasen:2018gtb ; Kovarik:2019xvh ; Armesto:2015lrg ; Eskola:2021nhw  as well,
albeit with significantly larger uncertainties compared to the case of proton
PDFs. In particular, the poorly constrained nuclear gluon PDF has been the
focus of recent nCTEQ studies Kusina:2020lyz ; Duwentaster:2021ioo , but
unfortunately the uncertainties have remained substantial. Through DGLAP
evolution, the gluons also produce significant uncertainties on other flavors.
Even the vector boson production and single inclusive hadron production data
sets, that were included in our recent nCTEQ15WZ(+SIH) studies Kusina:2020lyz
; Duwentaster:2021ioo , do not constrain the gluon below $x\approx 10^{-3}$.

In this paper, we perform a global analysis of nuclear PDFs in the nCTEQ
framework, including open heavy flavor meson and heavy quarkonium data. Heavy
quarkonium and open heavy-flavor meson production data (in the following
collectively called heavy quark production, HQ) has the potential to yield new
constraints on the gluon PDF, because the gluon-gluon channel contributes the
dominant part to the overall cross section of these processes. This was first
shown in the case of proton PDFs in Ref. PROSA:2015yid , before the
reweighting study presented in Ref. Kusina:2017gkz  was the first to
demonstrate that the LHC heavy-quark data is also useful to constrain the
nuclear gluon PDF at small $x$. The $x$-dependence of the data can be
estimated at leading order as

$\displaystyle x\approx\frac{2p_{T}}{\sqrt{s}}\exp(-|y|),$ (1)

where $y$ is the rapidity in the centre-of-mass frame111All rapidities
mentioned in the following are in the centre-of-mass frame and are therefore
denoted simply by $y$.. Fig. 1 shows the kinematic coverage of the available
data with contours for the corresponding $x$ dependence according to the
estimate above. This shows the data to be sensitive to $x$ values below
$10^{-5}$ for the most forward or backward rapidity at low $p_{T}$.

Figure 1: Coverage of the kinematic $(p_{T},y_{\mathrm{cms}})$-plane of the
quarkonium and open heavy quark production data sets from proton-lead
collisions. ALICE data is shown in red, ATLAS in blue, CMS in orange and LHCb
in green. The dashed and solid contours show the estimated $x$-dependence for
$\sqrt{s}=5$ and 8 TeV, respectively.

$D$-meson production data in particular has been used in recent nuclear PDF
analyses to reduce the uncertainty on the gluon PDF Eskola:2021nhw ;
Khalek:2022zqe . In Ref. Eskola:2021nhw , the ratio of double differential
cross sections between proton-proton and proton-lead is compared to pQCD
predictions in the GMVFNS Helenius:2018uul ; Kniehl:2004fy ; Kniehl:2005mk
and the data is included directly in the global analysis. In contrast, Ref.
Khalek:2022zqe  compares the same data to predictions in a fixed-flavor number
scheme with POWHEG and Pythia8 and uses this to perform a Bayesian reweighting
of the PDFs.

In this investigation we will study single inclusive222The meaning of the word
“inclusive” can be somewhat ambiguous in this context. Here, we use it to
denote processes, where the hadron is produced together with an arbitrary
number of other particles, as opposed to exclusive production. This is the
case for all processes studied in this paper. Later we will therefore use the
word “inclusive” to denote the sum of prompt and non-prompt production of a
particle. production of open heavy-flavor mesons and heavy quarkonia in
proton-proton and proton-lead collisions. Including these processes in a
global PDF fit is not as straight-forward as others, since there is no
universally accepted theoretical model for quarkonia. Therefore we introduce a
data-driven approach, that relies only on the following assumptions:

  * •

the gluon-gluon channel contributes the dominant fraction of the total cross
section,

  * •

it is sufficient to focus on subprocesses with 2 $\rightarrow$ 2 kinematics
and subprocesses with more than 2 hard final state particles can be neglected.

If these assumptions are valid, the data-driven approach has considerable
advantages over the other available calculations in perturbative QCD. Firstly,
the studied processes can be accurately described across a large kinematic
region with well controlled uncertainties. The second advantage is the speed
of the calculation, which is significantly faster than the available pQCD
calculations Kniehl:2009ar . Perturbative calculations at NLO would require
very efficient gridding to be included in a fit, where the predictions for
hundreds of data points need to be evaluated thousands of times.

The focus of this study will be on incorporating this process into a new
global analysis, while accounting for the uncertainty of the data-driven
theory, to determine the nuclear gluon PDF with greater accuracy than
previously possible. A precise knowledge of the nuclear PDFs is important for
the following reasons: Firstly, they provide a description of the hadronic or
nuclear structure in terms of quark and gluon degrees of freedom. In the
context of the standard pQCD formalism the PDFs are universal and therefore
required to make predictions for a wide range of collider observables.
Finally, they provide a starting point for comparisons with microscopic models
predicting the nuclear modifications (at the twist-2 level) in different $x$
regions. This includes microscopic models for nuclear effects on PDFs in the
shadowing region Armesto:2006ph ; Frankfurt:2011cs ; Kopeliovich:2012kw ;
Kulagin:2004ie , the antishadowing region Brodsky:1989qz ; Brodsky:2004qa ;
Kulagin:2004ie , or the EMC effect Geesaman:1995yd ; Norton:2003cb ;
Hen:2013oha ; Malace:2014uea ; Hen:2016kwk ; Kulagin:2004ie . At small-$x$ and
moderately hard scales, the density of gluons becomes very large such that the
assumptions underlying collinear factorization are expected to break down.
This kinematic region is described by the theory of Color Glass Condensates
Iancu:2000hn ; Gelis:2010nm  and there are also promising unified approaches
which interpolate between the CGC at small $x$ and collinear factorization at
large $x$, see, e.g. Ref. Jalilian-Marian:2020xvg  and references therein.
Nevertheless, it is fair to say that for now there is no unambiguous
microscopic picture of the inner workings of heavier nuclei.

It should be stressed again, that throughout this paper, our main underlying
assumption is that the twist-2 collinear factorization remains valid also in
the case of $eA$ and $pA$ collisions for the same observables. As it has been
discussed in Refs. Qiu:2003cg ; Accardi:2004be  this is reasonable, even if
higher twist terms may be enhanced in the nuclear case up to higher hard
scales $(\propto A^{1/3})$. We impose kinematic cuts on the data to
effectively reduce the impact of these higher twist effects and confirm
phenomenologically that all remaining data is well described. In the future,
such higher twist effects could be modelled to extend the reach towards data
with lower hard scales. One example is the effects due to fully coherent
energy loss Arleo:2012hn ; Arleo:2012rs ; Arleo:2021bpv . These contributions
are formally higher twist (twist 3), but have been shown to be relevant for
hard process data up to moderately large transverse momenta $p_{T}\approx
10\,\mathrm{GeV}$. It could therefore be interesting to include such effects
in future global analyses, however more work would be needed both on the
conceptual and the phenomenological side.

The next section provides an overview of the nCTEQ framework and the
integration of the new data-driven approach. Following that, we perform and
evaluate the fit of the proton-proton baseline for the theory in Sec. III. In
Sec. IV we present the fits obtained using the HQ data and evaluate the
compatibility between the new and old data. Finally, in Sec. V we summarize
our findings and give an outlook for future work.

##  II Theoretical approach

###  II.1 The nCTEQ framework

The nCTEQ project expands upon the foundation of the proton PDF global fitting
analysis by including the nuclear dimension. In early proton PDF analyses
(e.g. Ref. Olness:2003wz ), the nuclear data was used to calculate correction
factors which were then applied to the proton PDF fit without any
uncertainties. In contrast, the nCTEQ framework enables full communication
between nuclear and proton data, which means that observed tensions between
data sets can be investigated through the lens of nuclear corrections.

The details of the nCTEQ15 nPDFs are presented in Ref. Kovarik:2015cma . The
current analysis, along with the other recent nCTEQ analyses, such as
nCTEQ15WZ Kusina:2020lyz , nCTEQ15HIX Segarra:2020gtj  and nCTEQ15WZ+SIH
Duwentaster:2021ioo , is performed with a new C++-based code nCTEQ++. This
allows us to easily interface external programs such as HOPPET Salam:2008qg ,
APPLgrid Carli:2010rw , and INCNLO INCNLO . In particular, we work at leading
twist and next-to-leading order (NLO) of QCD for both the PDF and FF evolution
equations as well as the hard scattering coefficients. The calculation code
for the quarkonia and open heavy quarks is a partial C++ adaption of HELAC-
Onia 2.0 Shao:2015vga  and uses the data-driven approach explained in Sec.
II.2 instead of a pQCD calculation.

For the fits in this investigation, we use the same 19 parameters as for the
nCTEQ15WZ(+SIH) sets. These 19 parameters include the 16 free parameters of
the nCTEQ15 analysis, with an additional 3 open parameters for the strange
distribution. For the nCTEQ15 set, the strange PDF was constrained by the
relation $s=\bar{s}={(\kappa/2)(\bar{u}{+}\bar{d})}$ at the initial scale
$Q_{0}=1.3$ GeV, which forces it into the same form as the other sea quarks.

Our PDFs are parameterized at the initial scale $Q_{0}=1.3$ GeV as

$\displaystyle
xf_{i}^{p/A}(x,Q_{0})=c_{0}x^{c_{1}}(1-x)^{c_{2}}e^{c_{3}x}(1+e^{c_{4}}x)^{c_{5}}\
,$ (2)

and the nuclear $A$ dependence is encoded in the coefficients as

$\displaystyle c_{k}\longrightarrow c_{k}(A)\equiv p_{k}+a_{k}(1-A^{-b_{k}})\
,$ (3)

where $k=\\{1,...,5\\}$. The 16 free parameters used for the nCTEQ15 set
describe the $x$-dependence of the {$g,u_{v},d_{v},\bar{d}{+}\bar{u}$} PDF
combinations, and we do not vary the $\bar{d}/\bar{u}$ parameters; see Ref.
Kovarik:2015cma  for details. As in the nCTEQ15WZ(+SIH) analysis, we have
added three strange PDF parameters:
$\\{a_{0}^{s+\bar{s}},a_{1}^{s+\bar{s}},a_{2}^{s+\bar{s}}\\}$; these
parameters correspond to the nuclear modification of the overall
normalization, the low-$x$ exponent and the large-$x$ exponent of the strange
quark distribution, respectively.

In total, the 19 open parameters are:

$\displaystyle\\{a^{u_{v}}_{1},\ a^{u_{v}}_{2},\ a^{u_{v}}_{4},\
a^{u_{v}}_{5},\ a^{d_{v}}_{1},\ a^{d_{v}}_{2},\ a^{d_{v}}_{5},\
a^{\bar{u}+\bar{d}}_{1},\ a^{\bar{u}+\bar{d}}_{5},$ $\displaystyle a^{g}_{1},\
a^{g}_{4},\ a^{g}_{5},\ b^{g}_{0},\ b^{g}_{1},\ b^{g}_{4},\ b^{g}_{5},\
\bm{a^{s+\bar{s}}_{0},\ a^{s+\bar{s}}_{1},\ a^{s+\bar{s}}_{2}}\\}.$

All the fixed parameters are kept as they were in nCTEQ15.

###  II.2 The data-driven approach

Instead of performing the cross section calculations of the heavy mesons in
perturbative QCD, we take the data-driven approach outlined initially in Ref.
Kom:2011bd  and used for a reweighting study in Refs. Kusina:2017gkz ;
Kusina:2020dki . In this approach, the cross section for two nuclei $A$ and
$B$ scattering and producing a quarkonium or open heavy-flavor meson
$\mathcal{Q}$ is calculated as the convolution integral of the two initial
state gluon PDFs $f_{1,g}(x_{1},\mu),f_{2,g}(x_{2},\mu)$ and a fitted
effective scattering matrix element
$\overline{\left|\mathcal{A}_{gg\rightarrow\mathcal{Q}+X}\right|^{2}}$ over
the $AB\rightarrow\mathcal{Q}$ phase space

$\displaystyle\sigma(AB\rightarrow\mathcal{Q}+X)=$
$\displaystyle\int\mathrm{d}x_{1}\mathrm{~{}d}x_{2}f_{1,g}\left(x_{1},\mu\right)f_{2,g}\left(x_{2},\mu\right)\frac{1}{2\hat{s}}\overline{\left|\mathcal{A}_{gg\rightarrow\mathcal{Q}+X}\right|^{2}}\mathrm{dPS}.$

The effective scattering matrix element is parameterized with the Crystal Ball
function

$\displaystyle\begin{split}&\overline{\left|\mathcal{A}_{gg\rightarrow\mathcal{Q}+X}\right|^{2}}=\frac{\lambda^{2}\kappa\hat{s}}{M_{\mathcal{Q}}^{2}}e^{a|y|}\\\
&\times\begin{cases}e^{-\kappa\frac{p_{T}^{2}}{M_{\mathcal{Q}}^{2}}}&\text{ if
}p_{T}\leq\left\langle p_{T}\right\rangle\\\ e^{-\kappa\frac{\left\langle
p_{T}\right\rangle^{2}}{M_{\mathcal{Q}}^{2}}}\left(1+\frac{\kappa}{n}\frac{p_{T}^{2}-\left\langle
p_{T}\right\rangle^{2}}{M_{\mathcal{Q}}^{2}}\right)^{-n}&\text{ if
}p_{T}>\left\langle p_{T}\right\rangle\end{cases},\end{split}$ (4)

where the five parameters333The parameter name $``\left<p_{T}\right>"$ is
somewhat misleading. The parameterization was initially invented for a
different purpose, where this parameter did have the physical meaning of the
particle’s average transverse momentum, but this interpretation is lost in the
current context. However, we decided to keep the name to keep consistency with
previous works. $\lambda$, $\kappa$, $\left<p_{T}\right>$, $n$ and $a$ are
then fitted for each final state $\mathcal{Q}$. We have introduced the fifth
parameter $a$, which was not present in the original parameterization
Gaiser:1982 , to allow for a more accurate reproduction of the rapidity
dependence ellis2003qcd . The parameters are then fitted to
$pp\rightarrow\mathcal{Q}+X$ data. Once the optimal parameters are found, we
can also determine the uncertainty of our Crystal Ball fit via the same
Hessian method used to calculate our PDF uncertainties. We can then account
for these uncertainties by adding them in quadrature to the systematic
uncertainties of the $pPb\rightarrow\mathcal{Q}+X$ data. The included final
states in this analysis are $D^{0}$, $J/\psi$, $\Upsilon(1S)$ and $\psi(2S)$
mesons. Note, however, that prompt and non-prompt production of the same
particle need to be considered as two different final states. Inclusive
production is generally not fitted separately, but calculated as the sum of
the other two. The exception to this is $\Upsilon(1S)$, where all available
data is for inclusive production. Other final states, like $D^{\pm}$ or higher
excitations of $\Upsilon$, are excluded due to insufficient data, but in
principle the approach can be used for any measurement of single inclusive
hadrons, as long as sufficient proton-proton data is available for the
baseline and the cross section is dominated by the gluon-gluon channel.

The default scales $\mu=\mu_{0}$ that enter the PDFs are chosen as they were
in the previous reweighting study, as shown in Tab. 1. The entire procedure is
also repeated once with the scale doubled and once with the scale halved to
obtain an estimate of the impact that this choice has on the final result.

Table 1: Scale choices for the different particles. | $D^{0}$ | $J/\psi$ | $B\rightarrow J/\psi$ | $\Upsilon(1S)$ | $\psi(2S)$ | $B\rightarrow\psi(2S)$
---|---|---|---|---|---|---
$\mu_{0}^{2}$ | 4$M_{D}^{2}+p_{T,D}^{2}$ | $M_{J/\psi}^{2}+p_{T,J/\psi}^{2}$ | 4$M_{B}^{2}+\frac{M_{B}^{2}}{M_{J/\psi}^{2}}p_{T,J/\psi}^{2}$ | $M_{\Upsilon(1S)}^{2}+p_{T,\Upsilon(1S)}^{2}$ | $M_{\psi(2S)}^{2}+p_{T,\psi(2S)}^{2}$ | 4$M_{B}^{2}+\frac{M_{B}^{2}}{M_{\psi(2S)}^{2}}p_{T,\psi(2S)}^{2}$

##  III Proton-Proton baseline

In the first step of the analysis, we use the proton-proton data sets listed
in Tabs. 6 \- 6 to determine the parameters of the Crystal Ball function. The
largest portion of the proton-proton data comes from $J/\psi$ production,
followed by $\Upsilon(1S)$, then $\psi(2S)$ and finally $D^{0}$.

###  III.1 Cuts and excluded data

Both in the baseline fit and the PDF fit, we cut all heavy-quark data with
$p_{T}<3\,$GeV and outside of the rapidity range $-4<y_{cms}<4$. These cuts
allow predictions with reasonable $\chi^{2}$ for all remaining data points.
Relaxing the cuts introduces data points with
$\chi^{2}/N_{dof}=\mathcal{O}(10)$. We also exclude the lowest $p_{T}$ bin of
the 2011 CMS $J/\psi$ data sets, because the corresponding paper mentions that
there may be large acceptance effects in this bin that are not included in the
uncertainties. The lowest rapidity bin of the 7 TeV $\Upsilon(1S)$ production
data set from LHCb is included in the fit, but the normalization is determined
separately from the rest of the data set. This $y$-bin contains 20 $p_{T}$
bins (17 after cuts), which qualitatively agree with the remaining fit, but
the normalization is off by 25%. Therefore, keeping the normalization of this
bin the same as the remainder of the data set causes the $\chi^{2}/N_{dof}$
for the entire $\Upsilon(1S)$ fit to increase from 0.92 to 1.6. This
particular bin is also not described by other models based on the color-octet
mechanism LHCb:2015log .

We also limit the study to LHC data only, as the available RHIC data taken at
very different energies is not guaranteed to work with the same fit parameters
and would also likely not provide additional strong constraints.

###  III.2 Baseline fit

A comparison of the fits with the data is shown in Figs. 11 \- 16. The
uncertainties of the fit are determined via the same Hessian method that is
used in our PDF fits. All predictions show very close agreement with the data
across the region included in the fit. The large number of data points
produces very small uncertainties for $J/\psi$ and $\Upsilon(1S)$ (not visible
on the logarithmic scale). While the uncertainties of $D^{0}$ and $\psi(2S)$
are somewhat larger, due to the lower number of data points, they are still
small compared to the experimental uncertainties.

The obtained parameters and values of $\chi^{2}/N_{dof}{}$ for each process
are given in Tab. 2. The parameters for $J/\psi$ and $\psi(2S)$ are each
obtained in a combined fit with prompt, non-prompt and inclusive data, which
results in only one value of $\chi^{2}/N_{dof}$ for each. There is overall
very good agreement between the data and the fitted theory with
$\chi^{2}/N_{dof}$ values slightly below one for $J/\psi$, $\psi(2S)$ and
$\Upsilon(1S)$ and a particularly low $\chi^{2}/N_{dof}$ value of 0.25
obtained for $D^{0}$ production.

Alternative baseline fits have been performed with a variety of
parameterizations including more degrees of freedom, like the extended Crystal
Ball function (CB2) SignalExtractionALICE , which has a polynomial tail on
both sides of the Gaussian for a total of seven parameters. In those fits,
there is no significant improvement in the description of the currently
included data and relaxing the cuts still leads to unreasonable $\chi^{2}$
values for the data points that would be introduced. Therefore, we keep the
parameterization as written in Eq. (4) as a good compromise between a
reasonable number of parameters and accuracy of the data description.

Table 2: Crystal Ball parameters and $\chi^{2}/d.o.f.$ values for the Crystal Ball function for the different processes. | $D^{0}$ | $J/\psi$ | $B\rightarrow J/\psi$ | $\Upsilon(1S)$ | $\psi(2S)$ | $B\rightarrow\psi(2S)$
---|---|---|---|---|---|---
$\kappa$ | 0.33457 | 0.47892 | 0.15488 | 0.94524 | 0.21589 | 0.45273
$\lambda$ | 1.82596 | 0.30379 | 0.12137 | 0.06562 | 0.07528 | 0.13852
$\left<p_{T}\right>$ | 2.40097 | 5.29310 | -7.65026 | 8.63780 | 8.98819 | 7.80526
$n$ | 2.00076 | 2.17366 | 1.55538 | 1.93239 | 1.07203 | 1.64797
$a$ | -0.03295 | 0.02816 | -0.08083 | 0.22389 | -0.10614 | 0.06179
$N_{\mathrm{points}}$ | 34 | 501 | 375 | 55
$\chi^{2}/N_{dof}$ | 0.25 | 0.88 | 0.92 | 0.77
Table 3: Overview of the available $pp\ \longrightarrow\ D^{0}+X$ production
data sets and their number of data points.

Group | Year | Ref. | ID | Type | Points after / before cuts
---|---|---|---|---|---
ALICE | 2012 | ALICE:2012mhc  | 3008 | Prompt | 7 / 9
LHCb | 2013 | LHCb:2013xam  | 3014 | Prompt | 22 / 38
ALICE | 2016 | ALICE:2016yta  | 3021 | Prompt | 7 / 10

Table 4: Overview of the available $pp\ \longrightarrow\ J/\psi+X$ production
data sets and their number of data points. Note that the ALICE 2017 data set
is split into two for technical reasons (one part is taken at 5 TeV and one at
13 TeV).

Group | Year | Ref. | ID | Type | Points after / before cuts
---|---|---|---|---|---
ATLAS | 2011 | ATLAS:2011aqv  | 3003 | Prompt | 63 / 64
CMS | 2011 | CMS:2011rxs  | 3001 | Prompt | 39 / 44
LHCb | 2011 | LHCb:2011zfl  | 3002 | Prompt | 50 / 66
CMS | 2017 | CMS:2017exb  | 3020 | Prompt | 50 / 52
ATLAS | 2018 | ATLAS:2017prf  | 3017 | Prompt | 33 / 33
CMS | 2010 | CMS:2010nis  | 3007 | Non-prompt | 11 / 14
ATLAS | 2011 | ATLAS:2011aqv  | 3006 | Non-prompt | 63 / 64
CMS | 2011 | CMS:2011rxs  | 3005 | Non-prompt | 39 / 44
LHCb | 2011 | LHCb:2011zfl  | 3004 | Non-prompt | 50 / 66
CMS | 2017 | CMS:2017exb  | 3019 | Non-prompt | 50 / 52
ATLAS | 2018 | ATLAS:2017prf  | 3018 | Non-prompt | 33 / 33
ALICE | 2015 | ALICE:2015pgg  | 3022 | Inclusive | 10 / 13
ALICE | 2017 | ALICE:2017leg  | 3023 | Inclusive | 8 / 11
ALICE | 2017 | ALICE:2017leg  | 3024 | Inclusive | 14 / 18
ALICE | 2019 | ALICE:2019pid  | 3016 | Inclusive | 4 / 7

Table 5: Overview of the available $pp\ \longrightarrow\ \Upsilon(1S)+X$
production data sets and their number of data points. Note that the LHCb 2015
data set is split into two for technical reasons (one half is taken at 7 TeV
and one at 8 TeV).

Group | Year | Ref. | ID | Type | Points after / before cuts
---|---|---|---|---|---
ATLAS | 2012 | ATLAS:2012lmu  | 3012 | Inclusive | 88 / 100
LHCb | 2012 | LHCb:2012aa  | 3025 | Inclusive | 55 / 75
CMS | 2013 | CMS:2013qur  | 3013 | Inclusive | 30 / 42
ALICE | 2014 | ALICE:2014uja  | 3009 | Inclusive | 3 / 5
LHCb | 2015 | LHCb:2015log  | 3011 | Inclusive | 89 / 109
LHCb | 2015 | LHCb:2015log  | 3015 | Inclusive | 89 / 109
ALICE | 2015 | ALICE:2015pgg  | 3031 | Inclusive | 3 / 5
ATLAS | 2017 | ATLAS:2017prf  | 3031 | Inclusive | 18 / 24

Table 6: Overview of the available $pp\ \longrightarrow\ \psi(2S)+X$
production data sets and their number of data points.

Group | Year | Ref. | ID | Type | Points after / before cuts
---|---|---|---|---|---
ATLAS | 2017 | ATLAS:2017prf  | 3026 | Prompt | 15 / 15
ATLAS | 2017 | ATLAS:2017prf  | 3027 | Non-prompt | 15 / 15
ALICE | 2015 | ALICE:2015pgg  | 3028 | Inclusive | 6 / 9
CMS | 2018 | CMS:2018gbb  | 3029 | Inclusive | 9 / 10
ALICE | 2017 | ALICE:2017leg  | 3030 | Inclusive | 9 / 12

###  III.3 Comparison with $J/\psi$ production in NRQCD

The predictions for $J/\psi$ production from the data-driven method can also
be compared to predictions made using non-relativistic QCD (NRQCD). The NRQCD
framework assumes the speed of the heavy quarks to be slow compared to the
speed of light to derive a factorization theorem where the physics is
separated into short- and long-distance factors Bodwin:1994jh . The short-
distance physics, i.e. the production of the heavy quark-antiquark pair, can
then be calculated perturbatively, while the non-perturbative long distance
matrix elements parameterizing the formation of the bound state need to be
determined empirically. This has been done by various independent groups
including Ma et al. Ma:2010yw ; Ma:2010jj  and Butenschoen et al.
Butenschoen:2010rq ; Butenschoen:2011yh . The latter group has provided us
with predictions for prompt $J/\psi$ production using their NLO calculation
and our nCTEQ15 proton PDF, which is shown in Fig. 2. The uncertainties shown
for NRQCD are calculated by varying all scales together by a factor of two
around their default values $\mu_{r,0}=\mu_{f,0}$ =
$\sqrt{p_{T}^{2}+4m_{c}^{2}}$ and $m_{\mathrm{NRQCD},0}=m_{c}=1.5$ GeV. There
is a very good agreement between the two methods across the entire kinematic
range, but the Crystal Ball method produces significantly smaller
uncertainties. Varying the scales in the Crystal Ball fit produces only
negligible differences, as the scale dependence is mostly absorbed into the
parameters of the effective matrix element.

Figure 2: Comparison between prompt $J/\psi$ production in $pp$ collisions for
LHCbLHCb:2017ygo , ALICEALICE:2018mml  and ATLASATLAS:2015mpz  kinematics as
predicted by NRQCD and with the data-driven approach. The uncertainties of the
NRQCD predictions come from scale variation
$1/2<\mu_{r}/\mu_{r,0}=\mu_{f}/\mu_{f,0}=\mu_{\mathrm{NRQCD}}/\mu_{\mathrm{NRQCD},0}<2$
around the base scale $\mu_{r,0}=\mu_{f,0}$ = $\sqrt{p_{T}^{2}+4m_{c}^{2}}$
and $m_{\mathrm{NRQCD},0}=m_{c}$. Different rapidity bins are separated by
multiplying the cross sections by powers of ten for visual clarity.

###  III.4 Comparison with $D^{0}$ production in the GMVFNS

The predictions for $D^{0}$ production can also be compared with perturbative
calculations. These calculations can be carried out using the General-Mass
Variable-Flavor-Number-Scheme (GMVFNS) implementation of heavy quark
production at NLO QCD by Kniehl et al. Kniehl:2004fy ; Kniehl:2005mk . Fig. 3
shows a comparison of the predictions obtained from the GMVFNS code, with
those from our Crystal Ball fit for the data sets used in the fit for all the
$pp\to D^{0}+X$ data used in the baseline fit. In the input of GMVFNS we use
the nCTEQ15 proton PDF set, we set the $c$ quark mass to $m_{c}=1.3$ GeV, and
the renormalization and the initial/final factorization scales to
$\mu_{r}=\mu_{i}=\mu_{f}=\sqrt{p_{T}^{2}+4m_{c}^{2}}$. As a fragmentation
function we use the one with identifier 712 from the KKKS08 set of
fragmentation functions Kneesch:2007ey  which was obtained in a global fit to
Belle, CLEO, ALEPH and OPAL data. The uncertainties of the GMVFNS predictions
are obtained by varying the three scales individually by a factor of two, such
that there is never a factor four between two scales. These uncertainties are
similar in size as the data uncertainty, except for the low-$p_{T}$ region,
where they are somewhat larger. Overall the central prediction of the GMVFNS
calculation slightly overshoots the data. This can perhaps be attributed to
the contribution from largely unconstrained gluon component of the
fragmentation function, which contributes at almost 50%. However, there is
always overlap between the data and GMVFNS theory uncertainty. The uncertainty
of the Crystal Ball fit is similar in size as the GMVFNS one for large
$p_{T}$, but contrary to the latter it decreases for lower $p_{T}$ values. The
central values are very close to the data points, as indicated by the low
$\chi^{2}/N_{dof}$ value seen in Tab. 2. Overall the two methods are in very
good agreement with only minor discrepancies seen in the highest $p_{T}$ bins.

Figure 3: Comparison between prompt $D^{0}$ production as predicted in the
GMVFNS (red) and with the data-driven approach (blue). The uncertainties of
the GMVFNS predictions come from varying the scales individually by a factor
of 2, such that there is never a factor 4 between two scales. Different
rapidity bins are separated by multiplying the cross sections by powers of ten
for visual clarity.

We also compared our Crystall Ball fit and GMVFNS predictions against more
recent data of $D^{0}$ production in $pp$ collisions from ALICE ALICE:2017olh
; ALICE:2019nxm  and LHCb LHCb:2016ikn , which we have not been used in the
present analysis. We do not show the comparisons here, but we report that both
the Crystall Ball fit as well as GMVFNS reproduce the ALICE and LHCb data
well. This data could provide further constraints on the $D^{0}$ Crystall Ball
parameters.

##  IV Impact of HQ data on nPDF fits

Using the Crystal Ball parameters determined in the previous section we can
now perform a new global nPDF fit using the available heavy-quark data. The
new fits are using the same framework as nCTEQ15WZ+SIH, including all settings
like open parameters, scales and cuts for the previously included data. We do
not include the changes made for nCTEQ15HIX Segarra:2020gtj  and nCTEQ15$\nu$
nCTEQ15nu  as these developments are mostly orthogonal to those made in this
study and do not affect the low-$x$ gluon PDF. One minor change from the
previous analyses is the treatment of normalizations. Previously,
$\chi^{2}$-penalties were assigned individually for each affected data set,
whereas now they are applied only once per normalization parameter.

###  IV.1 Data selection

We add the heavy-quark data sets shown in Tabs. 10 \- 10 to the new PDF fit
for a total of 1484 (548 new, 936 old) data points. Similar to the
fragmentation function uncertainties of the SIH data in Ref.
Duwentaster:2021ioo , we can compensate for the theoretical uncertainty of the
data-driven approach by adding the uncertainty from the Crystal Ball fit as a
systematic uncertainty to all new data sets.

Table 7: Overview of the available $pPb\ \longrightarrow\ D^{0}+X$ production
data sets and their number of data points.

Group | Year | Ref. | ID | Type | Points after / before cuts
---|---|---|---|---|---
ALICE | 2014 | ALICE:2014xjz  | 3101 | Prompt | 8 / 10
ALICE | 2016 | ALICE:2016yta  | 3123 | Prompt | 8 / 11
LHCb | 2017 | LHCb:2017yua  | 3102 | Prompt | 53 / 92
ALICE | 2019 | ALICE:2019fhe  | 3122 | Prompt | 13 / 21

Table 8: Overview of the available $pPb\ \longrightarrow\ J/\psi+X$ production
data sets and their number of data points.

Group | Year | Ref. | ID | Type | Points after / before cuts
---|---|---|---|---|---
LHCb | 2013 | LHCb:2013gmv  | 3108 | Prompt | 25 / 40
ATLAS | 2015 | ATLAS:2015mpz  | 3118 | Prompt | 10 / 10
LHCb | 2017 | LHCb:2017ygo  | 3105 | Prompt | 88 / 140
CMS | 2017 | CMS:2017exb  | 3120 | Prompt | 51 / 53
ATLAS | 2018 | ATLAS:2017prf  | 3117 | Prompt | 8 / 8
LHCb | 2013 | LHCb:2013gmv  | 3107 | Non-prompt | 25 / 40
ATLAS | 2015 | ATLAS:2015mpz  | 3119 | Non-prompt | 10 / 10
LHCb | 2017 | LHCb:2017ygo  | 3106 | Non-prompt | 88 / 140
CMS | 2017 | CMS:2017exb  | 3121 | Non-prompt | 51 / 53
ATLAS | 2018 | ATLAS:2017prf  | 3116 | Non-prompt | 8 / 8
ALICE | 2013 | ALICE:2013snh  | 3103 | Inclusive | 0 / 12
ALICE | 2015 | ALICE:2015sru  | 3104 | Inclusive | 10 / 25
ALICE | 2018 | ALICE:2018mml  | 3112 | Inclusive | 9 / 24

Table 9: Overview of the available $pPb\ \longrightarrow\ \Upsilon(1S)+X$
production data sets and their number of data points.

Group | Year | Ref. | ID | Type | Points after / before cuts
---|---|---|---|---|---
ALICE | 2014 | ALICE:2014ict  | 3110 | Inclusive | 0 / 4
LHCb | 2014 | LHCb:2014rku  | 3111 | Inclusive | 0 / 2
ATLAS | 2018 | ATLAS:2017prf  | 3109 | Inclusive | 6 / 8
LHCb | 2018 | LHCb:2018psc  | 3113 | Inclusive | 36 / 66
ALICE | 2019 | ALICE:2019qie  | 3114 | Inclusive | 3 / 10

Table 10: Overview of the available $pPb\ \longrightarrow\ \psi(2S)+X$
production data sets and their number of data points.

Group | Year | Ref. | ID | Type | Points after / before cuts
---|---|---|---|---|---
ALICE | 2014 | ALICE:2014cgk  | 3127 | Inclusive | 2 / 8
ALICE | 2020 | ALICE:2020vjy  | 3126 | Inclusive | 3 / 10
ATLAS | 2017 | ATLAS:2017prf  | 3124 | Prompt | 8 / 8
CMS | 2018 | CMS:2018gbb  | 3115 | Prompt | 17 / 17
ATLAS | 2017 | ATLAS:2017prf  | 3125 | Non-prompt | 8 / 8

For the new HQ data, we use the same cuts as in the proton-proton baseline and
additionally exclude $D^{0}$ data points with $p_{T}>15$ GeV, because there is
no baseline data. Furthermore, we remove two individual points from the 2018
LHCb $\Upsilon(1S)$ data set that are described very poorly with $\chi^{2}$
values of 66 and 26, respectively. Both points are at the high-$p_{T}$ edge of
the experiment’s kinematic range, which makes systematic errors a likely
explanation, since the remaining 36 data points of the set can be well
described.

To estimate the impact of the scale choice on the final PDFs, we repeat the
entire procedure, including the proton-proton baseline fit, with the scale of
the heavy quark production processes set to $\frac{1}{2}$ or 2 times their
regular value.

###  IV.2 Resulting PDFs

Fig. 4 shows a comparison of the new fit, labelled nCTEQ15HQ444Since
nCTEQ15WZ+SIH+HQ would be too long, we break this naming convention and
shorten it to nCTEQ15HQ. This does not imply removal of the WZ and SIH data.,
with the PDFs obtained in previous nCTEQ analyses at a scale of $Q=2$ GeV. The
same comparison is shown in Fig. 5 in terms of nuclear modification factors.
The central value of the gluon PDF retains a similar shape as the
nCTEQ15WZ+SIH fit in the $x>10^{-3}$ region. Below this point, the
nCTEQ15WZ+SIH fit starts diverging towards higher values, while the new fit
moves towards similar values as is nCTEQ15WZ. From $x\approx 0.2$ downwards
the uncertainties of the new gluon PDF are reduced significantly compared to
the previous analyses, particularly below $x<10^{-4}$, where the previous
uncertainties start to increase rapidly. The ratio plots underline that there
are no significant changes to the up- and down-quark distributions, but there
is a minor reduction in $u_{\mathrm{valence}}$ around $x\approx 0.1$. The
uncertainties of the new fit still include the central values of all previous
fits. The gluon ratio on the other hand shows a reduction in the nuclear
modification for $x>0.008$. In this region, the gluon ratio is mostly
compatible with unity, but the central value goes from 0.8 at the lowest end
of this region to 1.2 at its peak and then down to 0.6 at very high $x$. At
lower $x$-values, however, the reduced uncertainties clearly show a
suppression of the gluon in lead between 20% and 40% compared to the proton
case. The modification of the strange quark is similar to the other quark
flavours in the new fit, but the uncertainties are still larger than the
modification for all $x<0.3$. The central values of the new up and down quark
PDFs are mostly very close to those of nCTEQ15WZ+SIH, with a minor downward
shift at low-$x$ and reduced uncertainties in the same region due to the
better constraints on the gluon. The strange quark PDF stays similar to
previous fits in the $x>0.01$ region, but reverts to a low value similar to
the original nCTEQ15, while nCTEQ15WZ and particularly nCTEQ15WZ+SIH start
increasing. It is important to note however, that the relative uncertainty of
the strange quark is still larger than 100% at low $x$, so that it is
impossible to draw any strong conclusions from this analysis. We refer the
reader to our recent neutrino analysis nCTEQ15nu  for complementary
information.

The impact of different scale choices is shown in Fig. 7, where each flavor
$i$ is shown as a ratio over the corresponding flavor $i_{\mathrm{central}}$
from the fit with the central scale choice. The central values of the up, down
and gluon PDFs with modified scales lie close together with differences no
larger than a few percent. The only exception to this is the high-$x$ region,
where the gluon PDF vanishes and therefore the relative uncertainties blow up.
The gluon PDFs cross each other at $x=0.02$ and $x=0.15$, with the lower scale
being above the central scale in the high- and low-$x$ regions and the upper
scale showing the opposite behaviour. The PDFs extracted with the modified
scale stay well within the uncertainties of the regular fit and show slightly
larger uncertainties themselves, which are likely due to a slightly worse fit.
The strange quark sees variations of up to 15%, but given the very large
uncertainties, this is still well within expectations. This conclusion is
somewhat different from the one reached in the previous reweighting analysis
Kusina:2017gkz ; Kusina:2020dki , where the scale choice resulted in
significant differences. Two updates in the methodology are responsible for
this change: Firstly, the fact that normalizations of the data sets are now
included as nuisance parameters in the fit, mitigates the scale dependence due
to the large normalization uncertainties on many data sets. Secondly, the more
restrictive kinematic cuts remove the data that is most sensitive to the scale
choice.

A final set of alternative fits is shown in Fig. 7, where the previous three
fits are compared to equivalent fits where the heavy-quark data sets keep
their bare uncertainties, i.e. the uncertainties from the Crystal Ball fit are
not added to the systematic uncertainties. For the up- and down-quark flavours
this leads to shifts on a similar order as the scale variation, while the
changes in the gluon are even smaller. The changes in the strange quark PDF
are somewhat larger due to it’s weakly constrained nature. The difference
between the dashed and solid lines gives an estimate of the potential
enhancement of the impact from the heavy-quark data if further proton-proton
data is added to the baseline to reduce the Crystal Ball uncertainty.

Figure 4: Lead PDFs from different nCTEQ15 versions. The baseline nCTEQ15 fit
is shown in black, nCTEQ15WZ in blue, nCTEQ15WZSIH in green, and the new fit
in red. Figure 5: Ratio of lead and proton PDF from different nCTEQ15
versions. The baseline nCTEQ15 fit is shown in black, nCTEQ15WZ in blue,
nCTEQ15WZSIH in green, and the new fit in red.

Figure 6: Comparisons of fits, where the scale for heavy quark production is
varied by a factor two around the central value $\mu_{0}$.

Figure 7: Comparisons of PDF central values with Crystal Ball uncertainties
included in the fit (solid lines) and without theory uncertainties (dashed
lines).

###  IV.3 Fit quality

Figs. 18 \- 21 show a comparison of the proton-lead data with predictions from
the new nCTEQ15HQ PDFs including PDF uncertainties. For all four particle
species, there is very close agreement across the entire kinematic range and
the PDF uncertainties are barely visible. Note however, that the uncertainties
from the Crystal Ball fit, which are not shown here, would still be of the
same size as in the baseline predictions in Fig. 11 \- 16.

For a more quantitative evaluation of the fit quality, one can take a look at
Fig. 8, which shows the $\chi^{2}/N_{dof}$ for each data set of the previous
nCTEQ15WZ+SIH and the new nCTEQ15HQ fits in the upper and lower panels,
respectively. The comparison shows no significant rise in $\chi^{2}$ for any
of the established data sets. The $\chi^{2}$ of the outlier 6215 (ATLAS Run I,
$Z$ production) among the vector boson production data sets is even reduced
due to the new normalization treatment. This means that there are no
incompatibilities between the previous and the new data sets. The new heavy-
quark data sets themselves mostly show $\chi^{2}/N_{dof}$ values below or
around one with a distribution similar to that of the DIS data sets. The only
major outlier among the new sets is the $\psi(2S)$ production from a 2014
measurement by ALICE with just two data points.

Figure 8: $\chi^{2}/N_{dof}{}$ values for each data set in the previous
nCTEQ15WZ+SIH fit (upper panel) and the new nCTEQ15HQ fit (lower panel).

Tab. 11 shows a comparison of the $\chi^{2}/N_{dof}$ values for each process
of the new fit and the preceding three nCTEQ15 generations. The original
nCTEQ15 is the only set of PDFs that does not give a good description of the
heavy-quark data, especially for $J/\psi$ mesons. nCTEQ15WZ already displays
significant improvement for $D^{0}$ and $J/\psi$, thus giving a reasonable
$\chi^{2}/N_{dof}$ of 0.92 for the new data sets. nCTEQ15WZ+SIH retains a very
similar overall $\chi^{2}$, but does worse than nCTEQ15WZ for $D^{0}$ and
$\Upsilon(1S)$, while improving the description of the $J/\psi$ data. Finally,
the new nCTEQ15HQ fit gives a good $\chi^{2}$ for all processes with only
minor increases for DIS and WZ production compared to both nCTEQ15WZ and
nCTEQ15WZ+SIH, which are outweighed by a significant improvement in the large
sample of new heavy-quark data. Interestingly, the $\chi^{2}/N_{dof}$ value
for $\psi(2S)$ production barely changes between the four fits, but since it
is the smallest of the new data sets with large theory uncertainties from the
baseline fit, the small impact is to be expected.

Table 11: $\chi^{2}/N_{dof}{}$ values for the individual heavy-quark final
states, the individual processes DIS, DY, WZ, SIH, HQ, and the total. The
shown $\chi^{2}$ is the sum of regular $\chi^{2}$ and normalization penalty.
Excluded processes are shown in parentheses. Note that both nCTEQ15 AND
nCTEQ15WZ included the neutral pions from STAR and PHENIX.
| $D^{0}$ | $J/\psi$ | $\Upsilon(1S)$ | $\psi(2S)$ | DIS | DY | WZ | SIH | HQ | Total
---|---|---|---|---|---|---|---|---|---|---
nCTEQ15 | (0.56) | (2.50) | (0.82) | (1.06) | 0.86 | 0.78 | (2.19) | (0.78) | (1.96) | 1.23
nCTEQ15WZ | (0.32) | (1.04) | (0.76) | (1.02) | 0.91 | 0.77 | 0.63 | (0.47) | (0.92) | 0.90
nCTEQ15WZ+SIH | (0.46) | (0.84) | (0.90) | (1.07) | 0.91 | 0.77 | 0.72 | 0.40 | (0.93) | 0.92
nCTEQ15HQ | 0.35 | 0.79 | 0.79 | 1.06 | 0.93 | 0.77 | 0.78 | 0.40 | 0.77 | 0.86

To estimate the impact of the new data on the PDF parameters, it is
instructive to look at the $\Delta\chi^{2}$ profiles along specific
parameters. Fig. 9 shows parameter scans for all of the seven gluon parameters
included in the fit with the $\Delta\chi^{2}$ split into processes. It is
immediately evident that the new heavy-quark data is the dominant constraint
on most of the parameters with the WZ production always constraining one side,
because its minimum is somewhat off-center. This indicates some slight
tensions with the WZ production data, but looking back at Tab. 11 this tension
largely comes from the SIH data and is still well within expectations. The
only parameter that is not constrained mostly by vector boson and heavy quark
production is $b_{0}^{g}$. This parameter determines the $A$ dependence of the
gluon normalization and minor shifts can therefore be compensated through the
normalization parameters of the data sets. Additionally, since parameters with
$b$ indices are related to the nuclear $A$ dependence, they may cancel out
when looking only at a singular $A$ as is the case in the LHC data, which is
only taken in proton-lead collisions.

Figure 9: Scans along the seven open gluon parameters, divided by experiment
types with the total shown in black.

##  V Conclusions

In conclusion, we have incorporated a large new data set of heavy quarkonium
and open heavy-flavor production in the nCTEQ++ framework and extended our PDF
analysis of the gluon to significantly lower $x$ values than was previously
possible.

We employed a data-driven approach to determine the theory and investigated
the advantages and limitations of this approach in detail. In particular, we
determined the kinematic range where the approach is applicable and verified
the predictions with those from rigorous pQCD calculations for $J/\psi$ and
$D^{0}$ production.

We obtained good $\chi^{2}/N_{dof}$ values for the new data sets without
compromising those of the established sets. The new data has a tremendous
impact on the gluon PDF, especially in the region $x<0.01$, and lowers the
uncertainty considerably. Through the DGLAP evolution, this impact on the
gluon also has consequences for the quark PDFs at low $x$, which remain at
similar central values, but also experience a reduction of their
uncertainties.

The data-driven approach has proven to be both reliable with a reasonable
estimation of its uncertainties, and also a powerful tool for the difficult
task of constraining the nuclear gluon PDF. Therefore, the presented analysis
and the PDFs obtained therein showcase a considerable advancement on the path
towards a precise understanding of nuclear structure.

## Acknowledgements

The work of P.D., T.J., M.K. and K.K. was funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) – project-id
273811115 – SFB 1225. We thank Peter Braun-Munzinger and Johanna Stachel for
stimulating this research and all our SFB collaborators for very useful
discussions. P.D., T.J., M.K., K.K. and K.F.M. also acknowledge support of the
DFG through the Research Training Group GRK 2149. The work of I.S. was
supported in part by the French National Centre for Scientific Research CNRS
through IN2P3 Project GLUE@NLO. F.O. acknowledges support through US DOE grant
DE-SC0010129. A.K. acknowledges the support of Narodowe Centrum Nauki under
Sonata Bis Grant No. 2019/34/E/ST2/00186.

##  Appendix A Data comparisons

This appendix contains comparisons between data and predictions made with the
Crystal Ball fit. Figs. 10-16 show the predictions for the fitted proton-
proton data with the uncertainties of the Crystal Ball fit, while Figs. 17-21
show those for the proton-lead data with the uncertainties of the nCTEQ15HQ
nuclear PDFs. Note that the proton-lead data is scaled by a factor
$\frac{1}{208}$ to give the cross section per nucleon, which is used in the
fit.

Figure 10: Predictions for $D^{0}$ production in proton-proton collisions with
uncertainties from the Crystal Ball fit. Different rapidity bins are separated
by multiplying the cross sections by powers of ten for visual clarity.

Figure 11: Predictions for $J/\psi$ production in proton-proton collisions
with uncertainties from the Crystal Ball fit. Different rapidity bins are
separated by multiplying the cross sections by powers of ten for visual
clarity.

Figure 12: Predictions for $J/\psi$ production in proton-proton collisions
with uncertainties from the Crystal Ball fit. Different rapidity bins are
separated by multiplying the cross sections by powers of ten for visual
clarity.

Figure 13: Predictions for $J/\psi$ production in proton-proton collisions
with uncertainties from the Crystal Ball fit. Different rapidity bins are
separated by multiplying the cross sections by powers of ten for visual
clarity.

Figure 14: Predictions for $\Upsilon(1S)$ production in proton-proton
collisions with uncertainties from the Crystal Ball fit. Different rapidity
bins are separated by multiplying the cross sections by powers of ten for
visual clarity.

Figure 15: Predictions for $\Upsilon(1S)$ production in proton-proton
collisions with uncertainties from the Crystal Ball fit. Different rapidity
bins are separated by multiplying the cross sections by powers of ten for
visual clarity.

Figure 16: Predictions for $\psi(2S)$ production in proton-proton collisions
with uncertainties from the Crystal Ball fit. Different rapidity bins are
separated by multiplying the cross sections by powers of ten for visual
clarity.

Figure 17: Predictions for $D^{0}$ production in proton-lead collisions with
PDF uncertainties of the nCTEQ15HQ fit. Different rapidity bins are separated
by multiplying the cross sections by powers of ten for visual clarity.

Figure 18: Predictions for $J/\psi$ production in proton-lead collisions with
PDF uncertainties of the nCTEQ15HQ fit. Different rapidity bins are separated
by multiplying the cross sections by powers of ten for visual clarity.

Figure 19: Predictions for $J/\psi$ production in proton-lead collisions with
PDF uncertainties of the nCTEQ15HQ fit. Different rapidity bins are separated
by multiplying the cross sections by powers of ten for visual clarity.

Figure 20: Predictions for $\Upsilon(1S)$ production in proton-lead collisions
with PDF uncertainties of the nCTEQ15HQ fit. Different rapidity bins are
separated by multiplying the cross sections by powers of ten for visual
clarity.

Figure 21: Predictions for $\psi(2S)$ production in proton-lead collisions
with PDF uncertainties of the nCTEQ15HQ fit. Different rapidity bins are
separated by multiplying the cross sections by powers of ten for visual
clarity.

## References

  * (1) J. C. Collins, D. E. Soper, and G. F. Sterman, “Factorization of Hard Processes in QCD,” Adv. Ser. Direct. High Energy Phys. 5 (1989) 1–91, hep-ph/0409313.
  * (2) T.-J. Hou et al., “New CTEQ global analysis of quantum chromodynamics with high-precision data from the LHC,” 1912.10053.
  * (3) NNPDF Collaboration, R. D. Ball et al., “Parton distributions from high-precision collider data,” Eur. Phys. J. C77 (2017), no. 10, 663, 1706.00428.
  * (4) R. Abdul Khalek, S. Bailey, J. Gao, L. Harland-Lang, and J. Rojo, “Towards Ultimate Parton Distributions at the High-Luminosity LHC,” Eur. Phys. J. C 78 (2018), no. 11, 962, 1810.03639.
  * (5) J. Gao, L. Harland-Lang, and J. Rojo, “The Structure of the Proton in the LHC Precision Era,” Phys. Rept. 742 (2018) 1–121, 1709.04922.
  * (6) K. Kovarik, P. M. Nadolsky, and D. E. Soper, “Hadron structure in high-energy collisions,” Rev.Mod.Phys. 92 (2020) 045003, 1905.06957.
  * (7) S. Alekhin, J. Blümlein, and S. Moch, “Strange sea determination from collider data,” Phys. Lett. B 777 (2018) 134–140, 1708.01067.
  * (8) P. M. Nadolsky, H.-L. Lai, Q.-H. Cao, J. Huston, J. Pumplin, D. Stump, W.-K. Tung, and C. P. Yuan, “Implications of CTEQ global analysis for collider observables,” Phys. Rev. D 78 (2008) 013004, 0802.0007.
  * (9) JAM Collaboration, N. Sato, C. Andres, J. J. Ethier, and W. Melnitchouk, “Strange quark suppression from a simultaneous Monte Carlo analysis of parton distributions and fragmentation functions,” Phys. Rev. D 101 (2020), no. 7, 074020, 1905.03788.
  * (10) L. Harland-Lang, A. Martin, P. Motylinski, and R. Thorne, “Parton distributions in the LHC era: MMHT 2014 PDFs,” Eur. Phys. J. C 75 (2015), no. 5, 204, 1412.3989.
  * (11) R. S. Thorne, S. Bailey, T. Cridge, L. A. Harland-Lang, A. Martin, and R. Nathvani, “Updates of PDFs using the MMHT framework,” PoS DIS2019 (2019) 036, 1907.08147.
  * (12) NNPDF Collaboration, R. D. Ball, L. Del Debbio, S. Forte, A. Guffanti, J. I. Latorre, A. Piccione, J. Rojo, and M. Ubiali, “Precision determination of electroweak parameters and the strange content of the proton from neutrino deep-inelastic scattering,” Nucl. Phys. B 823 (2009) 195–233, 0906.1958.
  * (13) H.-W. Lin et al., “Parton distributions and lattice QCD calculations: a community white paper,” Prog. Part. Nucl. Phys. 100 (2018) 107–160, 1711.07916.
  * (14) H.-W. Lin et al., “Parton distributions and lattice QCD calculations: toward 3D structure,” 2006.08636.
  * (15) K. Kovarik et al., “nCTEQ15 - Global analysis of nuclear parton distributions with uncertainties in the CTEQ framework,” Phys. Rev. D 93 (2016), no. 8, 085037, 1509.00792.
  * (16) K. J. Eskola, P. Paakkinen, H. Paukkunen, and C. A. Salgado, “EPPS16: Nuclear parton distributions with LHC data,” Eur. Phys. J. C 77 (2017), no. 3, 163, 1612.05741.
  * (17) NNPDF Collaboration, R. Abdul Khalek, J. J. Ethier, and J. Rojo, “Nuclear parton distributions from lepton-nucleus scattering and the impact of an electron-ion collider,” Eur. Phys. J. C79 (2019), no. 6, 471, 1904.00018.
  * (18) R. Abdul Khalek, J. J. Ethier, J. Rojo, and G. van Weelden, “nNNPDF2.0: Quark Flavor Separation in Nuclei from LHC Data,” 2006.14629.
  * (19) J. J. Ethier and E. R. Nocera, “Parton Distributions in Nucleons and Nuclei,” Ann. Rev. Nucl. Part. Sci. (2020), no. 70, 1–34, 2001.07722.
  * (20) V. Guzey and M. Klasen, “Constraints on nuclear parton distributions from dijet photoproduction at the LHC,” Eur. Phys. J. C 79 (2019), no. 5, 396, 1902.05126.
  * (21) M. Klasen, K. Kovarik, and J. Potthoff, “Nuclear parton density functions from jet production in DIS at an EIC,” Phys. Rev. D 95 (2017), no. 9, 094013, 1703.02864.
  * (22) M. Klasen and K. Kovařík, “Nuclear parton density functions from dijet photoproduction at the EIC,” Phys. Rev. D 97 (2018), no. 11, 114013, 1803.10985.
  * (23) N. Armesto, H. Paukkunen, J. M. Penín, C. A. Salgado, and P. Zurita, “An analysis of the impact of LHC Run I proton–lead data on nuclear parton densities,” Eur. Phys. J. C 76 (2016), no. 4, 218, 1512.01528.
  * (24) K. J. Eskola, P. Paakkinen, H. Paukkunen, and C. A. Salgado, “EPPS21: A global QCD analysis of nuclear PDFs,” 2112.12462.
  * (25) A. Kusina et al., “Impact of LHC vector boson production in heavy ion collisions on strange PDFs,” Eur. Phys. J. C 80 (2020), no. 10, 968, 2007.09100.
  * (26) nCTEQ Collaboration, P. Duwentäster, L. A. Husová, T. Ježo, M. Klasen, K. Kovařík, A. Kusina, K. F. Muzakka, F. I. Olness, I. Schienbein, and J. Y. Yu, “Impact of inclusive hadron production data on nuclear gluon PDFs,” 2105.09873.
  * (27) PROSA Collaboration, O. Zenaiev et al., “Impact of heavy-flavour production cross sections measured by the LHCb experiment on parton distribution functions at low x,” Eur. Phys. J. C 75 (2015), no. 8, 396, 1503.04581.
  * (28) A. Kusina, J.-P. Lansberg, I. Schienbein, and H.-S. Shao, “Gluon Shadowing in Heavy-Flavor Production at the LHC,” Phys. Rev. Lett. 121 (2018), no. 5, 052004, 1712.07024.
  * (29) R. A. Khalek, R. Gauld, T. Giani, E. R. Nocera, T. R. Rabemananjara, and J. Rojo, “nNNPDF3.0: Evidence for a modified partonic structure in heavy nuclei,” 2201.12363.
  * (30) I. Helenius and H. Paukkunen, “Revisiting the D-meson hadroproduction in general-mass variable flavour number scheme,” JHEP 05 (2018) 196, 1804.03557.
  * (31) B. A. Kniehl, G. Kramer, I. Schienbein, and H. Spiesberger, “Inclusive D*+- production in p anti-p collisions with massive charm quarks,” Phys. Rev. D 71 (2005) 014018, hep-ph/0410289.
  * (32) B. A. Kniehl, G. Kramer, I. Schienbein, and H. Spiesberger, “Collinear subtractions in hadroproduction of heavy quarks,” Eur. Phys. J. C 41 (2005) 199–212, hep-ph/0502194.
  * (33) B. A. Kniehl, G. Kramer, I. Schienbein, and H. Spiesberger, “Open charm hadroproduction and the charm content of the proton,” Phys. Rev. D 79 (2009) 094009, 0901.4130.
  * (34) N. Armesto, “Nuclear shadowing,” J. Phys. G 32 (2006) R367–R394, hep-ph/0604108.
  * (35) L. Frankfurt, V. Guzey, and M. Strikman, “Leading Twist Nuclear Shadowing Phenomena in Hard Processes with Nuclei,” Phys. Rept. 512 (2012) 255–393, 1106.2091.
  * (36) B. Z. Kopeliovich, J. G. Morfin, and I. Schmidt, “Nuclear Shadowing in Electro-Weak Interactions,” Prog. Part. Nucl. Phys. 68 (2013) 314–372, 1208.6541.
  * (37) S. A. Kulagin and R. Petti, “Global study of nuclear structure functions,” Nucl. Phys. A 765 (2006) 126–187, hep-ph/0412425.
  * (38) S. J. Brodsky and H. J. Lu, “Shadowing and Antishadowing of Nuclear Structure Functions,” Phys. Rev. Lett. 64 (1990) 1342.
  * (39) S. J. Brodsky, I. Schmidt, and J.-J. Yang, “Nuclear antishadowing in neutrino deep inelastic scattering,” Phys. Rev. D 70 (2004) 116003, hep-ph/0409279.
  * (40) D. F. Geesaman, K. Saito, and A. W. Thomas, “The nuclear EMC effect,” Ann. Rev. Nucl. Part. Sci. 45 (1995) 337–390.
  * (41) P. R. Norton, “The EMC effect,” Rept. Prog. Phys. 66 (2003) 1253–1297.
  * (42) O. Hen, D. W. Higinbotham, G. A. Miller, E. Piasetzky, and L. B. Weinstein, “The EMC Effect and High Momentum Nucleons in Nuclei,” Int. J. Mod. Phys. E 22 (2013) 1330017, 1304.2813.
  * (43) S. Malace, D. Gaskell, D. W. Higinbotham, and I. Cloet, “The Challenge of the EMC Effect: existing data and future directions,” Int. J. Mod. Phys. E 23 (2014), no. 08, 1430013, 1405.1270.
  * (44) O. Hen, G. A. Miller, E. Piasetzky, and L. B. Weinstein, “Nucleon-Nucleon Correlations, Short-lived Excitations, and the Quarks Within,” Rev. Mod. Phys. 89 (2017), no. 4, 045002, 1611.09748.
  * (45) E. Iancu, A. Leonidov, and L. D. McLerran, “Nonlinear gluon evolution in the color glass condensate. 1.,” Nucl. Phys. A 692 (2001) 583–645, hep-ph/0011241.
  * (46) F. Gelis, E. Iancu, J. Jalilian-Marian, and R. Venugopalan, “The Color Glass Condensate,” Ann. Rev. Nucl. Part. Sci. 60 (2010) 463–489, 1002.0333.
  * (47) J. Jalilian-Marian, “From small to large $x$: toward a unified description of high energy collisions,” J. Phys. Conf. Ser. 1602 (2020), no. 1, 012024.
  * (48) J.-w. Qiu, “QCD factorization and rescattering in proton nucleus collisions,” hep-ph/0305161.
  * (49) A. Accardi et al., “Hard probes in heavy ion collisions at the lhc: pdfs, shadowing and pa collisions,” in 3rd Workshop on Hard Probes in Heavy Ion Collisions: 3rd Plenary Meeting.  10, 2004.  hep-ph/0308248.
  * (50) F. Arleo and S. Peigne, “J/$\psi$ suppression in p-A collisions from parton energy loss in cold QCD matter,” Phys. Rev. Lett. 109 (2012) 122301, 1204.4609.
  * (51) F. Arleo and S. Peigne, “Heavy-quarkonium suppression in p-A collisions from parton energy loss in cold QCD matter,” JHEP 03 (2013) 122, 1212.0434.
  * (52) F. Arleo, G. Jackson, and S. Peigné, “Impact of fully coherent energy loss on heavy meson production in pA collisions,” JHEP 01 (2022) 164, 2107.05871.
  * (53) F. Olness, J. Pumplin, D. Stump, J. Huston, P. M. Nadolsky, H. L. Lai, S. Kretzer, J. F. Owens, and W. K. Tung, “Neutrino dimuon production and the strangeness asymmetry of the nucleon,” Eur. Phys. J. C 40 (2005) 145–156, hep-ph/0312323.
  * (54) E. P. Segarra et al., “nCTEQ15HIX – Extending nPDF Analyses into the High-$x$, Low $Q^{2}$ Region,” 2012.11566.
  * (55) G. P. Salam and J. Rojo, “A Higher Order Perturbative Parton Evolution Toolkit (HOPPET),” Comput. Phys. Commun. 180 (2009) 120–156, 0804.3755.
  * (56) T. Carli, D. Clements, A. Cooper-Sarkar, C. Gwenlan, G. P. Salam, F. Siegert, P. Starovoitov, and M. Sutton, “A posteriori inclusion of parton density functions in NLO QCD final-state calculations at hadron colliders: The APPLGRID Project,” Eur. Phys. J. C 66 (2010) 503–524, 0911.2985.
  * (57) M. Werlen, “INCNLO-direct photon and inclusive hadron production code website.” Version 1.4. http://lapth.cnrs.fr/PHOX_FAMILY.
  * (58) H.-S. Shao, “HELAC-Onia 2.0: an upgraded matrix-element and event generator for heavy quarkonium physics,” Comput. Phys. Commun. 198 (2016) 238–259, 1507.03435.
  * (59) C. H. Kom, A. Kulesza, and W. J. Stirling, “Pair Production of J/psi as a Probe of Double Parton Scattering at LHCb,” Phys. Rev. Lett. 107 (2011) 082002, 1105.4186.
  * (60) A. Kusina, J.-P. Lansberg, I. Schienbein, and H.-S. Shao, “Reweighted nuclear PDFs using heavy-flavor production data at the LHC,” Phys. Rev. D 104 (2021), no. 1, 014010, 2012.11462.
  * (61) J. Gaiser, “Charmonium spectroscopy from radiative decays of the $j/\psi$ and $\psi(2s)$,”. https://www.slac.stanford.edu/cgi-bin/getdoc/slac-r-255.pdf.
  * (62) R. Ellis, W. Stirling, and B. Webber, QCD and Collider Physics.  Cambridge Monographs on Particle Physics, Nuclear Physics and Cosmology. Cambridge University Press, 2003.
  * (63) LHCb Collaboration, R. Aaij et al., “Forward production of $\Upsilon$ mesons in $pp$ collisions at $\sqrt{s}=7$ and 8TeV,” JHEP 11 (2015) 103, 1509.02372.
  * (64) ALICE Collaboration, “Quarkonium signal extraction in alice.” https://cds.cern.ch/record/2060096/ files/QuarkoniumSignalExtraction_v5.pdf.
  * (65) ALICE Collaboration, B. Abelev et al., “Measurement of charm production at central rapidity in proton-proton collisions at $\sqrt{s}=7$ TeV,” JHEP 01 (2012) 128, 1111.1553.
  * (66) LHCb Collaboration, R. Aaij et al., “Prompt charm production in pp collisions at sqrt(s)=7 TeV,” Nucl. Phys. B 871 (2013) 1–20, 1302.2864.
  * (67) ALICE Collaboration, J. Adam et al., “$D$-meson production in $p$-Pb collisions at $\sqrt{s_{\rm NN}}=$5.02 TeV and in pp collisions at $\sqrt{s}=$7 TeV,” Phys. Rev. C 94 (2016), no. 5, 054908, 1605.07569.
  * (68) ATLAS Collaboration, G. Aad et al., “Measurement of the differential cross-sections of inclusive, prompt and non-prompt $J/\psi$ production in proton-proton collisions at $\sqrt{s}=7$ TeV,” Nucl. Phys. B 850 (2011) 387–444, 1104.3038.
  * (69) CMS Collaboration, S. Chatrchyan et al., “$J/\psi$ and $\psi_{2S}$ production in $pp$ collisions at $\sqrt{s}=7$ TeV,” JHEP 02 (2012) 011, 1111.1557.
  * (70) LHCb Collaboration, R. Aaij et al., “Measurement of $J/\psi$ production in $pp$ collisions at $\sqrt{s}=7~{}\rm{TeV}$,” Eur. Phys. J. C 71 (2011) 1645, 1103.0423.
  * (71) CMS Collaboration, A. M. Sirunyan et al., “Measurement of prompt and nonprompt $\mathrm{J}/{\psi}$ production in $\mathrm{p}\mathrm{p}$ and $\mathrm{p}\mathrm{Pb}$ collisions at $\sqrt{s_{\mathrm{NN}}}=5.02\,\text{TeV}$,” Eur. Phys. J. C 77 (2017), no. 4, 269, 1702.01462.
  * (72) ATLAS Collaboration, M. Aaboud et al., “Measurement of quarkonium production in proton–lead and proton–proton collisions at $5.02~{}\mathrm{TeV}$ with the ATLAS detector,” Eur. Phys. J. C 78 (2018), no. 3, 171, 1709.03089.
  * (73) CMS Collaboration, V. Khachatryan et al., “Prompt and Non-Prompt $J/\psi$ Production in $pp$ Collisions at $\sqrt{s}=7$ TeV,” Eur. Phys. J. C 71 (2011) 1575, 1011.4193.
  * (74) ALICE Collaboration, J. Adam et al., “Inclusive quarkonium production at forward rapidity in pp collisions at $\sqrt{s}=8$ TeV,” Eur. Phys. J. C 76 (2016), no. 4, 184, 1509.08258.
  * (75) ALICE Collaboration, S. Acharya et al., “Energy dependence of forward-rapidity $\mathrm{J}/\psi$ and $\psi\mathrm{(2S)}$ production in pp collisions at the LHC,” Eur. Phys. J. C 77 (2017), no. 6, 392, 1702.00557.
  * (76) ALICE Collaboration, S. Acharya et al., “Inclusive J/$\psi$ production at mid-rapidity in pp collisions at $\sqrt{s}$ = 5.02 TeV,” JHEP 10 (2019) 084, 1905.07211.
  * (77) ATLAS Collaboration, G. Aad et al., “Measurement of Upsilon production in 7 TeV pp collisions at ATLAS,” Phys. Rev. D 87 (2013), no. 5, 052004, 1211.7255.
  * (78) LHCb Collaboration, R. Aaij et al., “Measurement of Upsilon production in pp collisions at $\sqrt{s}$ = 7 TeV,” Eur. Phys. J. C 72 (2012) 2025, 1202.6579.
  * (79) CMS Collaboration, S. Chatrchyan et al., “Measurement of the $\Upsilon(1S),\Upsilon(2S)$, and $\Upsilon(3S)$ Cross Sections in $pp$ Collisions at $\sqrt{s}$ = 7 TeV,” Phys. Lett. B 727 (2013) 101–125, 1303.5900.
  * (80) ALICE Collaboration, B. B. Abelev et al., “Measurement of quarkonium production at forward rapidity in $pp$ collisions at $\sqrt{s}=7$ TeV,” Eur. Phys. J. C 74 (2014), no. 8, 2974, 1403.3648.
  * (81) CMS Collaboration, A. M. Sirunyan et al., “Measurement of prompt $\psi$(2S) production cross sections in proton-lead and proton-proton collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}=$ 5.02 TeV,” Phys. Lett. B 790 (2019) 509–532, 1805.02248.
  * (82) G. T. Bodwin, E. Braaten, and G. P. Lepage, “Rigorous QCD analysis of inclusive annihilation and production of heavy quarkonium,” Phys. Rev. D 51 (1995) 1125–1171, hep-ph/9407339. [Erratum: Phys.Rev.D 55, 5853 (1997)].
  * (83) Y.-Q. Ma, K. Wang, and K.-T. Chao, “$J/\psi(\psi^{\prime})$ production at the Tevatron and LHC at ${\cal O}(\alpha_{s}^{4}v^{4})$ in nonrelativistic QCD,” Phys. Rev. Lett. 106 (2011) 042002, 1009.3655.
  * (84) Y.-Q. Ma, K. Wang, and K.-T. Chao, “A complete NLO calculation of the $J/\psi$ and $\psi^{\prime}$ production at hadron colliders,” Phys. Rev. D 84 (2011) 114001, 1012.1030.
  * (85) M. Butenschoen and B. A. Kniehl, “Reconciling $J/\psi$ production at HERA, RHIC, Tevatron, and LHC with NRQCD factorization at next-to-leading order,” Phys. Rev. Lett. 106 (2011) 022003, 1009.5662.
  * (86) M. Butenschoen and B. A. Kniehl, “World data of J/psi production consolidate NRQCD factorization at NLO,” Phys. Rev. D 84 (2011) 051501, 1105.0820.
  * (87) LHCb Collaboration, R. Aaij et al., “Prompt and nonprompt J/$\psi$ production and nuclear modification in $p$Pb collisions at $\sqrt{s_{\text{NN}}}=8.16$ TeV,” Phys. Lett. B 774 (2017) 159–178, 1706.07122.
  * (88) ALICE Collaboration, S. Acharya et al., “Inclusive J/$\psi$ production at forward and backward rapidity in p-Pb collisions at $\sqrt{s_{\rm NN}}$ = 8.16 TeV,” JHEP 07 (2018) 160, 1805.04381.
  * (89) ATLAS Collaboration, G. Aad et al., “Measurement of differential $J/\psi$ production cross sections and forward-backward ratios in p + Pb collisions with the ATLAS detector,” Phys. Rev. C 92 (2015), no. 3, 034904, 1505.08141.
  * (90) T. Kneesch, B. A. Kniehl, G. Kramer, and I. Schienbein, “Charmed-meson fragmentation functions with finite-mass corrections,” Nucl. Phys. B 799 (2008) 34–59, 0712.0481.
  * (91) ALICE Collaboration, S. Acharya et al., “Measurement of D-meson production at mid-rapidity in pp collisions at ${\sqrt{s}=7}$ TeV,” Eur. Phys. J. C 77 (2017), no. 8, 550, 1702.00766.
  * (92) ALICE Collaboration, S. Acharya et al., “Measurement of ${{\mathrm{D}}^{0}}$ , ${{\mathrm{D}}^{+}}$ , ${{\mathrm{D}}^{*+}}$ and ${{\mathrm{D}}^{+}_{\mathrm{s}}}$ production in pp collisions at ${\sqrt{{\textit{s}}}~{}=~{}5.02~{}{\text{TeV}}}$ with ALICE,” Eur. Phys. J. C 79 (2019), no. 5, 388, 1901.07979.
  * (93) LHCb Collaboration, R. Aaij et al., “Measurements of prompt charm production cross-sections in pp collisions at $\sqrt{s}=5$ TeV,” JHEP 06 (2017) 147, 1610.02230.
  * (94) nCTEQ Collaboration, K. F. Muzakka, P. Duwentäster, T. Hobbs, T. Ježo, M. Klasen, K. Kovařík, A. Kusina, J. Morfin, F. I. Olness, R. Ruiz, I. Schienbein, and J. Y. Yu, “Compatibility of neutrino DIS data and its impact on nuclear parton distribution functions,”. In preparation.
  * (95) ALICE Collaboration, B. B. Abelev et al., “Measurement of prompt $D$-meson production in $p-Pb$ collisions at $\sqrt{s_{NN}}$ = 5.02 TeV,” Phys. Rev. Lett. 113 (2014), no. 23, 232301, 1405.3452.
  * (96) LHCb Collaboration, R. Aaij et al., “Study of prompt D0 meson production in $p$Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5$ TeV,” JHEP 10 (2017) 090, 1707.02750.
  * (97) ALICE Collaboration, S. Acharya et al., “Measurement of prompt D0, D+, D∗+, and ${\mathrm{D}}_{\mathrm{S}}^{+}$ production in p–Pb collisions at $\sqrt{{\mathrm{s}}_{\mathrm{NN}}}$ = 5.02 TeV,” JHEP 12 (2019) 092, 1906.03425.
  * (98) LHCb Collaboration, R. Aaij et al., “Study of $J/\psi$ production and cold nuclear matter effects in $pPb$ collisions at $\sqrt{s_{NN}}=5$ TeV,” JHEP 02 (2014) 072, 1308.6729.
  * (99) ALICE Collaboration, B. B. Abelev et al., “$J/\psi$ production and nuclear effects in p-Pb collisions at $\sqrt{S_{NN}}$ = 5.02 TeV,” JHEP 02 (2014) 073, 1308.6726.
  * (100) ALICE Collaboration, J. Adam et al., “Rapidity and transverse-momentum dependence of the inclusive J/$\psi$ nuclear modification factor in p-Pb collisions at $\sqrt{s_{N\ N}}=$ 5.02 TeV,” JHEP 06 (2015) 055, 1503.07179.
  * (101) ALICE Collaboration, B. B. Abelev et al., “Production of inclusive $\Upsilon$(1S) and $\Upsilon$(2S) in p-Pb collisions at $\mathbf{\sqrt{s_{{\rm NN}}}=5.02}$ TeV,” Phys. Lett. B 740 (2015) 105–117, 1410.2234.
  * (102) LHCb Collaboration, R. Aaij et al., “Study of $\Upsilon$ production and cold nuclear matter effects in $p$Pb collisions at $\sqrt{s_{NN}}$=5 TeV,” JHEP 07 (2014) 094, 1405.5152.
  * (103) LHCb Collaboration, R. Aaij et al., “Study of $\Upsilon$ production in $p$Pb collisions at $\sqrt{s_{NN}}=8.16$ TeV,” JHEP 11 (2018) 194, 1810.07655. [Erratum: JHEP 02, 093 (2020)].
  * (104) ALICE Collaboration, S. Acharya et al., “$\Upsilon$ production in p–Pb collisions at $\sqrt{s_{NN}}$=8.16 TeV,” Phys. Lett. B 806 (2020) 135486, 1910.14405.
  * (105) ALICE Collaboration, B. B. Abelev et al., “Suppression of $\psi$(2S) production in p-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV,” JHEP 12 (2014) 073, 1405.3796.
  * (106) ALICE Collaboration, S. Acharya et al., “Measurement of nuclear effects on $\psi\rm{(2S)}$ production in p-Pb collisions at $\sqrt{\textit{s}_{\rm NN}}=8.16$ TeV,” JHEP 07 (2020) 237, 2003.06053.