# Fixed-target charmonium production and pion parton distributions

Wen-Chen Chang  Institute of Physics, Academia Sinica, Taipei 11529, Taiwan
Jen-Chieh Peng  Department of Physics, University of Illinois at Urbana-
Champaign, Urbana, Illinois 61801, USA  Stephane Platchkov  IRFU, CEA,
Université Paris-Saclay, 91191 Gif-sur-Yvette, France  Takahiro Sawada  Nambu
Yoichiro Institute of Theoretical and Experimental Physics, Osaka Metropolitan
University, Osaka 558-8585, Japan

(August 27, 2024)

###### Abstract

We investigate how charmonium hadroproduction at fixed-target energies can be
used to constrain the gluon distribution in pion. Using non-relativistic QCD
(NRQCD) formulation, the $J/\psi$ and $\psi(2S)$ cross sections as a function
of longitudinal momentum fraction $x_{F}$ from pions and protons colliding
with light targets, as well as the $\psi(2S)$ to $J/\psi$ cross section
ratios, are included in the analysis. The color-octet long-distance matrix
elements are found to have a pronounced dependence on the pion parton
distribution functions (PDFs). This study shows that the $x_{F}$ differential
cross sections of pion-induced charmonium production impose strong constraints
on the pion’s quark and gluon PDFs. In particular, the pion PDFs with larger
gluon densities provide a significantly better description of the data. It is
also found that the production of the $\psi(2S)$ state is associated with a
larger quark-antiquark contribution, compared with $J/\psi$.

##  I Introduction

The pion, as the lightest QCD bound state, plays an essential role in the
nucleon-nucleon interactions over nuclear-size distances Horn and Roberts
(2016). Theoretically, its partonic structure is easier to construct than that
of the nucleon. Pion distribution amplitudes and parton distribution functions
(PDFs) have been predicted by a number of recent calculations based on the
chiral-quark model Nam (2012); Watanabe  _et al._ (2016, 2018), Nambu-Jona-
Lasinio model Hutauruk  _et al._ (2016), light-front Hamiltonian Lan  _et al._
(2019, 2020, 2022), holographic QCD de Teramond  _et al._ (2018); Watanabe
_et al._ (2020); Lan and Mondal (2020), maximum entropy method Han  _et al._
(2020, 2021), Dyson-Schwinger equations (DSE) Chang  _et al._ (2014); Chang
and Thomas (2015); Chen  _et al._ (2016); Shi  _et al._ (2018); Bednar  _et
al._ (2020); Ding  _et al._ (2020); Cui  _et al._ (2020); Freese  _et al._
(2021); Chang and Roberts (2021); Cui  _et al._ (2022a, b), and lattice QCD
Zhang  _et al._ (2019); Sufian  _et al._ (2019); Izubuchi  _et al._ (2019);
Joó  _et al._ (2019); Sufian  _et al._ (2020); Chen  _et al._ (2020); Gao  _et
al._ (2020); Alexandrou  _et al._ (2021a, b); Fan and Lin (2021); Detmold  _et
al._ (2022); Gao  _et al._ (2022). In contrast, the partonic structure of pion
is much less explored experimentally, due to the absence of a pion target. The
present knowledge on the pion PDFs comes primarily from fixed-target pion-
induced Drell-Yan (DY) measurements Chang and Dutta (2013). However, the DY
data are mainly sensitive to the valence-quark distributions, leaving the sea
and gluon distributions essentially unknown. The sea-quark contributions can
in principle be extracted by comparing measurements with the positive and
negative pion beams Londergan  _et al._ (1995), although the existing
measurements are scarce and of insufficient statistical accuracy.

The gluon distribution in the pion can be accessed through processes such as
prompt-photon production Bonesini  _et al._ (1988), leading-neutron deep-
inelastic scattering (DIS) Khoze  _et al._ (2006); McKenney  _et al._ (2016)
or heavy quarkonia production Gluck  _et al._ (1978); Barger  _et al._ (1980).
Each of these processes has its own advantages and limitations. With the
exception of Ref. Owens (1984), the pion-induced $J/\psi$ and $\psi(2S)$
production data were not included in the global analysis, possibly reflecting
the concern that the production mechanism for charmonium production was not
well understood. Significant progress in understanding the $J/\psi$ production
mechanism has been made in recent decades, and it is timely to investigate how
the charmonium production data can provide useful constraints on the pion
PDFs.

The color evaporation model (CEM) Einhorn and Ellis (1975); Fritzsch (1977);
Halzen (1977), although successful for some observables, fails to explain some
others observables in charmonium production Bodwin  _et al._ (2005). Within
the more rigorous nonrelativistic QCD (NRQCD) Bodwin  _et al._ (1995)
framework, the production of the heavy quark pair is treated perturbatively,
whereas its hadronization to a bound state is described in terms of a set of
long-distance matrix elements (LDMEs), extracted phenomenologically from the
data.

From the experimental perspective, charmonium production has one important
advantage: the cross sections are large, between one to two orders of
magnitude higher than the DY ones, depending on the experimental conditions. A
large number of fixed-target charmonium production experiments have been
performed in the past, including experiments with pion beams Schuler (1994);
Vogt (1999). These data, collected mostly at CERN or at Fermilab, provide a
wealth of additional information on the pion structure, and are expected to
shed new light on its gluon distribution.

In this paper we investigate how charmonium production could help to
differentiate between the available pion PDFs by imposing further constraints
on the gluon distribution function Chang  _et al._ (2020); Hsieh  _et al._
(2021). In the fixed-target energy domain, charmonium production is dominated
by the quark-antiquark annihilation ($q\bar{q}$) and gluon-gluon fusion ($GG$)
partonic subprocesses. The longitudinal momentum $x_{F}$-differential cross
sections is sensitive to the quark and gluon parton distributions of the
colliding hadrons. Since the nucleon PDFs are known with good accuracy, these
differential cross sections should provide additional constraints on the
pion’s quark and gluon PDFs.

Figure 1: Momentum density distributions $[xf(x)]$ of valence quarks, sea
quarks and gluons of SMRS, GRV, xFitter and JAM pion PDFs and their ratios to
the SMRS PDFs, at the scale of $J/\psi$ mass ($Q^{2}$= 9.6 GeV2). The
uncertainty bands associated with JAM and xFitter PDFs are also shown.

To perform this study, we employ the NRQCD framework, along the lines
developed in Ref. Beneke and Rothstein (1996). Although limited to leading
order (LO), this approach provides an adequate description of the fixed-target
$J/\psi$ and $\psi(2S)$ production data and can be used as a tool for
accessing the pion PDFs. Our primary goal is to obtain a good phenomenological
description of both pion and proton-induced data, and to explore the
sensitivity of the results to the pion quark-gluon structure. Assuming that
the LDMEs are independent of the beam species, the proton-induced cross
sections are also included in this analysis. Since the proton PDFs are well
known, the proton-induced data should help constraining the values of LDMEs
common to both the proton and pion data.

Results of an earlier study limited to total cross sections of charmonium
production were recently reported Hsieh  _et al._ (2021). A new set of color-
octet LDMEs, leading to a good agreement between the charmonium production
data and the NRQCD fit was obtained. Here, we extend the study by including
the $x_{F}$-dependent cross sections for $J/\psi$ and $\psi(2S)$ production as
well as their ratios for both pion and proton beams in the global fit. In
order to minimize nuclear matter effects that are not well understood, the
present analysis is limited to data taken with the lightest targets available:
hydrogen, lithium and beryllium. Data with heavier targets were considered
only for the $J/\psi$ to $\psi(2S)$ cross section ratios, assuming nuclear
effects are largely independent of the charmonium states.

This paper is organized as follows. In Sec. II, we describes distinctive
features of parton densities in four pion PDFs. The NRQCD formalism used for
this study is introduced in Sec. III. Section IV briefly describes the
$J/\psi$ and $\psi(2S)$ data sets used in the global fit. We present the
results of NRQCD calculations using various pion PDFs and the comparison with
the charmonium data in Sec. V. Finally we comment on the fit results in Sec.
VI and conclude in Sec. VII

##  II Pion PDFs

As mentioned before, pion-induced Drell-Yan data are used in all global
analyses for constraining the valence-quark distribution of the pion PDFs.
Without data from other processes, the sea and gluon distributions can only be
inferred through the momentum sum rule and valence-quark sum rule. The two
most recent global analyses dedicated to the extraction of the pion PDFs are
JAM Barry  _et al._ (2018); Cao  _et al._ (2021); Barry  _et al._ (2021) and
xFitter Novikov  _et al._ (2020). The two groups consider the same DY data,
but differ in the choice of the additional processes. The xFitter group makes
use of the pion-induced prompt-photon production data, whereas the JAM
collaboration includes the leading-neutron DIS cross section measurements
instead. The Sutton-Martin-Roberts-Stirling (SMRS) global fit Sutton  _et al._
(1992) also incorporates the prompt-photon data, but instead of calculating
the fit uncertainties, it considers three different options for the gluon and
sea contents. Another widely used parametrization is the fit of Gluck-Reya-
Vogt (GRV) Gluck  _et al._ (1992), in which the gluon and sea distributions
are dynamically generated from the QCD evolution.

We utilize the LHAPDF framework Whalley  _et al._ (2005); Buckley  _et al._
(2015) to access these four pion PDFs for our study. The corresponding pion
PDF sets are “SMRSPI.LHgrid”, “GRVPI1”, “JAM21PionPDFnlo”, and
“xFitterPI_NLO_EIG”, respectively. Out of the three possible parametrizations
for SMRS, we choose the one in which the sea quarks carry 15% of the pion
momentum at $Q^{2}$= 4 GeV2. Their valence, sea and gluon momentum
distributions $xf(x)$ at the scale of $J/\psi$ mass are compared in Fig. 1.
Their ratios to SMRS are shown in the bottom panel. Within the range of
$x\sim$0.1–0.8, the valence-quark distributions of SMRS, JAM and xFitter are
close to each other, whereas GRV is lower by up to 20%–30%. Not surprisingly,
the sea distribution is essentially unknown, as illustrated by the large
variations between the four PDFs. The gluon distributions also show sizable
differences; e.g., in the region of $x>0.2$ the xFitter and JAM distributions
are smaller in comparison with SMRS and GRV, by up to a factor of 2-3.

##  III Heavy-Quark Pair Production and NRQCD Model

Within the NRQCD theoretical framework, the heavy quarkonium production is
factorized into production of a heavy-quark pair ($Q\bar{Q}$) at the parton
level, and its subsequent hadronization into quarkonium states. The $Q\bar{Q}$
production cross section can be calculated perturbatively Nason  _et al._
(1988, 1989); Mangano  _et al._ (1993), whereas the hadronization probability
of the $Q\bar{Q}$ pair is encoded in the nonperturbative LDME parameters
$\langle\mathcal{O}_{n}^{H}[^{2S+1}L_{J}]\rangle$, depending on the spin,
orbital , and total angular momentum quantum numbers, $S$, $L$ and $J$,
respectively, and on the color configuration ($n$). Parity, charge conjugation
and angular momentum conservation limit the allowed quantum numbers to only a
few. The LDMEs are assumed to be universal, i.e., independent of the beam and
target hadrons and of the energy scale. The color singlet (CS) LDMEs are
typically determined from decay rate measurements using a potential model
Eichten and Quigg (1995), while the color octet (CO) LDMEs are obtained from a
fit to the experimental data.

In NRQCD, the differential cross section $d\sigma/dx_{F}$ for the production
of a charmonium state $H$ ($H$ = $J/\psi$, $\psi(2S)$, or $\chi_{cJ}$) from
the $hN$ collisions, where $h$ is the beam hadron ($h$ = $p$, $\bar{p}$, or
$\pi$) and $N$ the target nucleon, is expressed as Vogt (2000)

$\displaystyle\frac{d\sigma^{H}}{dx_{F}}=$
$\displaystyle\sum\limits_{i,j=q,\bar{q},G}\int_{0}^{1}dx_{1}dx_{2}\delta(x_{F}-x_{1}+x_{2})$
$\displaystyle\times$ $\displaystyle
f^{h}_{i}(x_{1},\mu_{F})f^{N}_{j}(x_{2},\mu_{F})$ $\displaystyle\times$
$\displaystyle\hat{\sigma}[ij\rightarrow
H](x_{1}P_{h},x_{2}P_{N},\mu_{F},\mu_{R},m_{c}),$ (1)
$\displaystyle\hat{\sigma}[ij\rightarrow H]=$
$\displaystyle\sum\limits_{n}C^{ij}_{c\bar{c}[n]}(x_{1}P_{h},x_{2}P_{N},\mu_{F},\mu_{R},m_{c})$
$\displaystyle\times$
$\displaystyle\langle\mathcal{O}_{n}^{H}[^{2S+1}L_{J}]\rangle$ (2)

where the indexes $i$ and $j$ run over the type of interacting partons
(gluons, quarks and antiquarks), and $C^{ij}_{c\bar{c}[n]}$ denotes the hard-
QCD production cross section for $c\bar{c}$ pair. The parameter $m_{c}$ is the
charm quark mass; $f^{h}$ and $f^{N}$ are the incoming hadron and the target
nucleon parton distribution functions, evaluated at their respective
Bjorken-$x$ values, $x_{1}$ and $x_{2}$. The $\mu_{F}$ and $\mu_{R}$ are the
factorization and renormalization scales. The Feynman variable $x_{F}$ and the
beam and target parton momentum fractions $x_{1}$ and $x_{2}$ are:

$\displaystyle x_{F}=\frac{2p_{L}}{\sqrt{s}}\mbox{,
}x_{1,2}=\frac{(x_{F}^{2}+4{M_{c\bar{c}}}^{2}/s)^{1/2}\pm x_{F}}{2}.$ (3)

Here $M_{c\bar{c}}$ and $p_{L}$ are the mass and longitudinal momentum of the
$c\bar{c}$ pair in the center-of-mass frame. The total cross sections are
obtained by integrating over $x_{F}$.

In this study, we use the formula given in Ref. Beneke and Rothstein (1996)
for computation of $J/\psi$, $\psi(2S)$, and $\chi_{cJ}$ production via $GG$,
$q\bar{q}$ and $qG$ subprocesses. The scattering subprocesses $q\bar{q}\to
Q\bar{Q}$ and $GG\to Q\bar{Q}$ at $\mathcal{O}(\alpha_{s}^{2})$ produce
$Q\bar{Q}$ pairs in an $S$-wave CO state or $P$-wave CS state. Table 1
summarizes the relationships between the LDMEs and the scattering subprocesses
for $J/\psi$, $\psi(2S)$, $\chi_{c0}$, $\chi_{c1}$, and $\chi_{c2}$, up to
$\mathcal{O}(\alpha_{s}^{3})$. For the $q\bar{q}$ subprocess, the $c\bar{c}$
pairs are produced at $\mathcal{O}(\alpha_{s}^{2})$ in color octet states,
which then hadronize into various charmonium states with the LDMEs
$\langle\mathcal{O}_{8}^{H}[^{3}S_{1}]\rangle$. For the $GG$ subprocess, both
$J/\psi$, and $\psi(2S)$ can be produced from either the CO $c\bar{c}$ at
$\mathcal{O}(\alpha_{s}^{2})$ or the CS $c\bar{c}$ at
$\mathcal{O}(\alpha_{s}^{3})$. The CO ${}^{1}S_{0}$, ${}^{3}P_{0}$, and
${}^{3}P_{2}$ are combined into a single LDME, $\Delta_{8}^{H}$, via the
relation:
$\Delta_{8}^{H}=\langle\mathcal{O}_{8}^{H}[^{1}S_{0}]\rangle+\frac{3}{m_{c}^{2}}\langle\mathcal{O}_{8}^{H}[^{3}P_{0}]\rangle+\frac{4}{5m_{c}^{2}}\langle\mathcal{O}_{8}^{H}[^{3}P_{2}]\rangle$.

$H$ | $q\bar{q}$ | $GG$ | $qG$
---|---|---|---
$J/\psi$, $\psi(2S)$ |  $\langle\mathcal{O}_{8}^{H}[^{3}S_{1}]\rangle$ ($\alpha_{s}^{2}$) |  $\Delta_{8}^{H}$ ($\alpha_{s}^{2}$) |
|  |  $\langle\mathcal{O}_{1}^{H}[^{3}S_{1}]\rangle$ ($\alpha_{s}^{3}$) |
$\chi_{c0}$ |  $\langle\mathcal{O}_{8}^{H}[^{3}S_{1}]\rangle$ ($\alpha_{s}^{2}$) |  $\langle\mathcal{O}_{1}^{H}[^{3}P_{0}]\rangle$ ($\alpha_{s}^{2}$) |
$\chi_{c1}$ |  $\langle\mathcal{O}_{8}^{H}[^{3}S_{1}]\rangle$ ($\alpha_{s}^{2}$) |  $\langle\mathcal{O}_{1}^{H}[^{3}P_{1}]\rangle$ ($\alpha_{s}^{3}$) |  $\langle\mathcal{O}_{1}^{H}[^{3}P_{1}]\rangle$ ($\alpha_{s}^{3}$)
$\chi_{c2}$ |  $\langle\mathcal{O}_{8}^{H}[^{3}S_{1}]\rangle$ ($\alpha_{s}^{2}$) |  $\langle\mathcal{O}_{1}^{H}[^{3}P_{2}]\rangle$ ($\alpha_{s}^{2}$) |
Table 1: Relationship of LDMEs and the associated orders of $\alpha_{s}$ to
the scattering subprocesses for various charmonium states in the NRQCD
framework of Ref. Beneke and Rothstein (1996). Here
$\Delta_{8}^{H}=\langle\mathcal{O}_{8}^{H}[^{1}S_{0}]\rangle+\frac{3}{m_{c}^{2}}\langle\mathcal{O}_{8}^{H}[^{3}P_{0}]\rangle+\frac{4}{5m_{c}^{2}}\langle\mathcal{O}_{8}^{H}[^{3}P_{2}]\rangle$.

The number of independent LDMEs is further reduced by applying the spin
symmetry relations Beneke and Rothstein (1996); Maltoni  _et al._ (2006):

$\displaystyle\langle\mathcal{O}_{8}^{J/\psi,\psi(2S)}[^{3}P_{J}]\rangle$
$\displaystyle=(2J+1)\langle\mathcal{O}_{8}^{J/\psi,\psi(2S)}[^{3}P_{0}]\rangle\mbox{
for $J=2$}$
$\displaystyle\langle\mathcal{O}_{8}^{\chi_{cJ}}[^{3}S_{1}]\rangle$
$\displaystyle=(2J+1)\langle\mathcal{O}_{8}^{\chi_{c0}}[^{3}S_{1}]\rangle\mbox{
for $J=1,2$}$
$\displaystyle\langle\mathcal{O}_{1}^{\chi_{cJ}}[^{3}P_{J}]\rangle$
$\displaystyle=(2J+1)\langle\mathcal{O}_{1}^{\chi_{c0}}[^{3}P_{0}]\rangle\mbox{
for $J=1,2$}.$ (4)

The LDMEs used in the present work exhibit sensitivity to different elementary
scattering subprocesses contributing to the charmonium production. In the
cases of $J/\psi$ and $\psi(2S)$ production, the CO
$\langle\mathcal{O}_{8}^{H}[^{3}S_{1}]\rangle$ LDME is related to the
$q\bar{q}\to Q\bar{Q}$ subprocess, while the $GG\to Q\bar{Q}$ subprocess is
strongly dependent on the $\Delta_{8}^{H}$ term. More details on the NRQCD
framework used in this work can be found in Refs. Beneke and Rothstein (1996);
Hsieh  _et al._ (2021). In the following study, the CS
$\langle\mathcal{O}_{1}^{H}[^{3}S_{1}]\rangle$ LDMEs for $J/\psi$ and
$\psi(2S)$ and the CS $\langle\mathcal{O}_{1}^{H}[^{3}P_{0}]\rangle$ and CO
$\langle\mathcal{O}_{8}^{H}[^{3}S_{1}]\rangle$ LDMEs for $\chi_{c0}$ are fixed
to be 1.16, 0.76, 0.044 and 0.0032, respectively, which are the values used in
Refs. Beneke and Rothstein (1996); Hsieh  _et al._ (2021).

With the information of LDMEs, the direct production cross sections of
$J/\psi$, $\psi(2S)$ and three $\chi_{cJ}$ states as a function of $x_{F}$ can
be evaluated as shown in Eq.(III). The $J/\psi$ cross section is estimated
taking into account the direct production of $J/\psi$ and the feed-down from
hadronic decays of $\psi(2S)$ and radiative decays of three $\chi_{cJ}$ states
as follows,

$\displaystyle\sigma_{J/\psi}$ $\displaystyle=\sigma_{J/\psi}^{direct}$
$\displaystyle+Br(\psi(2S)\rightarrow J/\psi X)\sigma_{\psi(2S)}$
$\displaystyle+\sum\limits_{J=0}^{2}Br(\chi_{cJ}\rightarrow
J/\psi\gamma)\sigma_{\chi_{cJ}}$ (5)

The various branching ratios $Br$ are taken from the PDG 2020 Zyla  _et al._
(2020): $Br(\psi(2S)\rightarrow J/\psi X)=61.4\%$, $Br(\chi_{c0}\rightarrow
J/\psi\gamma)=1.4\%$, $Br(\chi_{c1}\rightarrow J/\psi\gamma)=34.3\%$, and
$Br(\chi_{c2}\rightarrow J/\psi\gamma)=19.0\%$.

In the present analysis we use the convention of charm quark mass,
factorization and renormalization scales in Ref. Beneke and Rothstein (1996)
for fixed-target hadroproduction of charmonium: $m_{c}=1.5$ GeV/$c^{2}$ and
$\mu_{F}=\mu_{R}=2m_{c}$. The uncertainties associated with this choice are
evaluated by changing the reference scale from $m_{c}$ to $3m_{c}$. The
nucleon PDFs are taken from CTEQ14nlo Dulat  _et al._ (2016). For the lithium
and beryllium targets, the nuclear EPPS16 PDFs Eskola  _et al._ (2017) is
used.

##  IV OVERVIEW OF DATA USED

The present analysis is based on pion and proton-induced total and
differential cross sections for $J/\psi$ and $\psi(2S)$ production, and on the
differential
$R_{\psi}(x_{F})=\sigma_{\psi(2S)}(x_{F})/\sigma_{(J/\psi)}(x_{F})$ ratios.
The total cross sections for the pion-induced data were taken from the
compilations made in Refs. Schuler (1994) and Vogt (1999). The proton-induced
total cross sections and ratios were taken from Ref. Maltoni  _et al._ (2006).
The proton-induced values for $R_{\psi}$ from HERA-B Abt  _et al._ (2007) and
NA38 Baglin  _et al._ (1995) and the pion-induced ones from WA92 Alexandrov
_et al._ (1999) and WA39 Corden  _et al._ (1980) were added to the selection.
The $x_{F}$-differential cross sections for pion-induced $J/\psi$ production
Gribushin  _et al._ (1996); Antoniazzi  _et al._ (1992); Badier  _et al._
(1983); McEwen  _et al._ (1983); Akerlof  _et al._ (1993); Corden  _et al._
(1981) and $\psi(2S)$ production Gribushin  _et al._ (1996) were selected
according to the targets used: hydrogen, lithium and beryllium. Data sets with
heavier targets were not included. The same criterion was applied to the
proton-induced $J/\psi$ production Antoniazzi  _et al._ (1992); Badier  _et
al._ (1983). The $R_{\psi}(x_{F})$ ratios were taken from Ref. Heinrich  _et
al._ (1991) for the pion-induced production and from Refs. Abt  _et al._
(2007); Alessandro  _et al._ (2004); Schub  _et al._ (1995); Alexopoulos  _et
al._ (1996) for the proton-induced one. Assuming that nuclear effects are
identical for both charmonium states, no restriction on the target employed
was applied.

The data sets with $x_{F}$ dependent measurements are listed in Table 2. In
terms of pion-induced (proton-induced) data sets, there are 8 (2) for $J/\psi$
production, 2 (0) for $\psi(2S)$ production and 1 (4) for $R_{\psi}(x_{F})$.
In total, there are 164 and 82 data points for the pion-induced and proton-
induced data, respectively. The beam momenta of the data sets cover the range
of 39.5–515 GeV/$c$, corresponding to $\sqrt{s}$ values ranging from 8.6 to
31.1 GeV.

Experiment | Beam |  $P_{beam}$ (GeV/$c$) | Target | Data | $x_{F}$ | ndf | Norma.a | Ref.
---|---|---|---|---|---|---|---|---
FNAL E672, E706 | $\pi$ | 515 | Be | $\sigma^{J/\psi}$ | [0.11, 0.79] | 35 | 12.0 | Gribushin  _et al._ (1996)
FNAL E705 | $\pi$ | 300 | Li | $\sigma^{J/\psi}$ | [-0.10, 0.45] | 12 | 9.5 | Antoniazzi  _et al._ (1992)
CERN NA3b | $\pi$ | 280 | p | $\sigma^{J/\psi}$ | [0.025, 0.825] | 17 | 13.0 | Badier  _et al._ (1983)
CERN NA3b | $\pi$ | 200 | p | $\sigma^{J/\psi}$ | [0.05, 0.75] | 8 | 13.0 | Badier  _et al._ (1983)
CERN WA11b | $\pi$ | 190 | Be | $\sigma^{J/\psi}$ | [-0.35, 0.75] | 12 |  c10.0 | McEwen  _et al._ (1983)
CERN NA3b | $\pi$ | 150 | p | $\sigma^{J/\psi}$ | [0.025, 0.925] | 19 | 13.0 | Badier  _et al._ (1983)
FNAL E537 | $\pi$ | 125 | Be | $\sigma^{J/\psi}$ | [0.05, 0.95] | 10 | 6.0 | Akerlof  _et al._ (1993)
CERN WA39b | $\pi$ | 39.5 | p | $\sigma^{J/\psi}$ | [0.05, 0.85] | 9 | 15.0 | Corden  _et al._ (1981)
FNAL E672, E706 | $\pi$ | 515 | Be | $\sigma^{\psi(2S)}$ | [0.17, 0.73] | 5 | 16.0 | Gribushin  _et al._ (1996)
FNAL E615 | $\pi$ | 253 | W | $\sigma^{\psi(2S)}/\sigma^{J/\psi}$ | [0.275, 0.975] | 15 |  | Heinrich  _et al._ (1991)
HERA-B | p | 920 | W | $\sigma^{\psi(2S)}/\sigma^{J/\psi}$ | [-0.3, 0.075] | 8 |  | Abt  _et al._ (2007)
CERN NA50 | p | 450 | W | $\sigma^{\psi(2S)}/\sigma^{J/\psi}$ | [-0.075, 0.075] | 4 |  | Alessandro  _et al._ (2004)
FNAL E789 | p | 800 | Au | $\sigma^{\psi(2S)}/\sigma^{J/\psi}$ | [0.00, 0.12] | 5 |  | Schub  _et al._ (1995)
FNAL E771 | p | 800 | Si | $\sigma^{\psi(2S)}/\sigma^{J/\psi}$ | [0.00, 0.20] | 6 |  | Alexopoulos  _et al._ (1996)
FNAL E705 | p | 300 | Li | $\sigma^{J/\psi}$ | [-0.10, 0.45] | 12 | 10.1 | Antoniazzi  _et al._ (1992)
CERN NA3b | p | 200 | p | $\sigma^{J/\psi}$ | [0.05, 0.75] | 8 | 13.0 | Badier  _et al._ (1983)
Table 2: Differential cross sections datasets for charmonium production
($J/\psi$, $\psi(2S)$ and $R_{\psi}(x_{F})$) used in the study, listed in
order of decreasing beam momentum.
aPercentage of uncertainty in the cross section normalization.
bThe numerical information was extracted from the published figures.
cInformation not available but an educated guess .

##  V Results of NRQCD calculations

###  V.1 Reference NRQCD calculations

| SMRS | GRV | JAM | xFitter
---|---|---|---|---
| REF | FIT | REF | FIT | REF | FIT | REF | FIT
$\chi^{2}_{total}/\text{ndf}$ | 5.7 | 1.9 | 7.0 | 2.4 | 17.7 | 5.6 | 14.3 | 4.2
$\chi^{2}/\text{ndp}|^{\pi^{-}}_{x_{F}}$ | 5.3 | 1.8 | 7.6 | 2.4 | 25.5 | 5.9 | 19.5 | 4.5
$\chi^{2}/\text{ndp}|^{p}_{x_{F}}$ | 10.7 | 1.6 | 10.5 | 1.7 | 11.2 | 2.7 | 11.5 | 1.9
$\chi^{2}/\text{ndp}|^{\pi^{-}}_{\sqrt{s}}$ | 2.1 | 8.7 | 2.9 | 5.6 | 5.3 | 11.4 | 4.8 | 4.4
$\chi^{2}/\text{ndp}|^{p}_{\sqrt{s}}$ | 3.8 | 8.1 | 3.4 | 8.1 | 3.5 | 5.1 | 3.6 | 6.9
$\langle\mathcal{O}_{8}^{J/\psi}[^{3}S_{1}]\rangle$ | 0.0690 | 0.0259$\pm$0.0023 | 0.0950 | 0.0432$\pm$0.0038 | 0.0830 | 0.1192$\pm$0.0021 | 0.0740 | 0.0849$\pm$0.0041
$\Delta_{8}^{J/\psi}$ | 0.0250 | 0.0560$\pm$0.0016 | 0.0180 | 0.0521$\pm$0.0017 | 0.0200 | 0.0244$\pm$0.0016 | 0.0220 | 0.0393$\pm$0.0034
$\langle\mathcal{O}_{8}^{\psi(2S)}[^{3}S_{1}]\rangle$ | 0.0210 | 0.0132$\pm$0.0009 | 0.0260 | 0.0210$\pm$0.0013 | 0.0260 | 0.0237$\pm$0.0009 | 0.0230 | 0.0186$\pm$0.0012
$\Delta_{8}^{\psi(2S)}$ | 0.0017 | 0.0057$\pm$0.0003 | 0.0004 | 0.0042$\pm$0.0003 | 0.0004 | 0.0021$\pm$0.0003 | 0.0009 | 0.0040$\pm$0.0006
Table 3: Results of the NRQCD calculation using the reference values of the
LDMEs (columns labeled “REF”) and of the fit of the LDMEs to the differential
cross sections (columns “FIT”). The upper part of the table gives the values
of the reduced $\chi^{2}/\text{ndf}$ of the entire data set and the $\chi^{2}$
divided by the number of data point (ndp) for the pion-induced and proton-
induced data sets separately. The subscript $x_{F}$ or $\sqrt{s}$ for
$\chi^{2}/\text{ndp}$ refers to $x_{F}$-dependent data or $\sqrt{s}$-dependent
$x_{F}$-integrated data. The lower part of the table displays the values of
the reference and fitted LDMEs for SMRS, GRV, JAM and xFitter pion PDFs. All
LDMEs are in units of $\rm{GeV}^{3}$.

Before performing a fit to the data listed in Table 2 to obtain the best-fit
LDMEs for the four pion PDFs, we first carry out NRQCD calculations using the
LDMEs found in a recent study Hsieh  _et al._ (2021), where only the pion and
proton total cross section data were fitted. The values of the LDMEs, obtained
separately for each pion PDF, are listed in Table 3. We then compare the
results of the NRQCD calculations for the $x_{F}$ dependent charmonium
production cross sections with the data listed in Table 2 and shown in Fig. 2.
We call these “Reference NRQCD calculation” (REF), which provides the
reference information to be compared with that obtained later from a fit to
the $x_{F}$-dependent cross section data. Note that Fig. 2 is for the SMRS
pion PDFs, and similar figures for the other three pion PDFs can be found in
the supplemental material Sup (2022).

The total $\chi^{2}$/ndf, as well as the $\chi^{2}$/ndp (ndp denotes “number
of data points”) for individual pion or proton data sets, are listed in Table
3 under the label “REF”. Table 3 shows that the reduced $\chi^{2}$/ndf for
“REF” are quite large, suggesting that the LDMEs deduced from the fit to total
cross section data are not optimal for describing the $x_{F}$-dependent data.
A further investigation shows that a significant contribution to the overall
$\chi^{2}$ comes from the absolute normalization of the measured cross
sections relative to the NRQCD calculations. Despite the poor agreement
between the data and the calculation, it is interesting to note that
calculations using the SMRS and GRV pion PDFs are in a better agreement with
the data than the JAM and xFitter PDFs.

###  V.2 NRQCD FITS

We now proceed to a refined determination of the color-octet
$\langle\mathcal{O}_{8}^{H}[^{3}S_{1}]\rangle$ and $\Delta_{8}^{H}$ LDMEs for
$J/\psi$ and $\psi(2S)$ production by fitting the $x_{F}$ differential cross
sections and $R_{\psi}(x_{F})$ ratios for proton and pion beams. To avoid
double counting, total cross sections data that result from an integration
over the associated differential cross sections are not included in the fit.
We note that the NRQCD calculations do not require a normalization factor, as
they predict absolute cross sections. However, the experimental
$x_{F}$-dependent $J/\psi$ and $\psi(2S)$ cross sections are associated with
experimental normalization uncertainties $\delta_{\sigma}$, as quoted in Table
2. An attempt to fit the data without taking into account the normalization
uncertainties only marginally reduces the total $\chi^{2}/\text{ndf}$. In
order to take into account these uncertainties, a normalization parameter $F$
is added for each of the $x_{F}$-differential data sets. Accordingly, a
penalty term of $((F-1)/\delta_{\sigma})^{2}$ is included in the calculation
of the overall $\chi^{2}$. To avoid unrealistic values of $F$, we limit the
deviation of $F$ from 1.0 to be less than 2 $\delta_{\sigma}$. The results of
this approach are labeled as “FIT” below.

Figure 2: The $x_{F}$-dependent cross sections for $J/\psi$ and $\psi(2S)$
production and $R_{\psi}(x_{F})$ ratios in $\pi^{-}N$ and $pN$ interactions,
following the order given in Table. 2. The symbol and value in parenthesis
denote the particle type and momentum of beam. The solid red and dotted black
curves represent the NRQCD results of SMRS pion PDFs from the fit described in
the text (“FIT”) and from the calculation using the LDMEs obtained in Ref.
Hsieh  _et al._ (2021) (“REF”), respectively. The values of $\chi^{2}$ divided
by the number of data point (ndp) for each data set are also shown.

Figure 2 shows the new fit to the data for $x_{F}$-differential data and
ratios using the SMRS pion PDFs. The newly-determined LDMEs parameters are
shown in Table 3. Except for the $J/\psi$ data of WA11, the new NRQCD fit
provides a reasonably good description of data for both pion and proton beams.
Table 3 shows that for all four pion PDFs and for nearly all data sets the
individual $\chi^{2}/\text{ndp}$ are significantly improved. Similar figures
for GRV, JAM and xFitter are available in the supplemental material Sup
(2022).

Table 3 also lists the $\chi^{2}$ values for both the “REF” and “FIT”
calculations. The $\chi^{2}/\text{ndp}$ and the fitted normalization factors
for each data set are summarized in Table 4. The improved description of the
differential cross sections is also confirmed by the overall
$\chi^{2}_{total}/\text{ndf}$ and the $\chi^{2}/\text{ndp}$ values for various
data sets. The $\chi^{2}$/ndp of the pion-induced $x_{F}$ data sets,
$\chi^{2}/\text{ndp}|^{\pi}_{x_{F}}$, are 1.8, 2.4, 5.9 and 4.5 for the SMRS,
GRV, JAM and xFitter PDFs, respectively, an improvement of about a factor of
three over that of “REF”. As expected, the $\chi^{2}$/ndp of the proton-
induced $x_{F}$ data sets, $\chi^{2}/\text{ndp}|^{p}_{x_{F}}$, are of similar
values, around 2.0 for all four pion PDFs. In contrast, the
$\chi^{2}/\text{ndp}$ of the integrated cross sections
($\chi^{2}/\text{ndp}|^{\pi^{-},p}_{\sqrt{s}}$ ) are now larger since these
data are not included in the global fit.

Table 3 also shows the newly fitted LDMEs. In comparison with the “REF”
calculation, the “FIT” results give smaller
$\langle\mathcal{O}_{8}^{H}[^{3}S_{1}]\rangle$ values for both SMRS and GRV
PDFs, while the corresponding $\Delta_{8}^{H}$ LDMEs are slightly larger. For
the JAM and xFitter PDFs the “REF” and “FIT” LDMEs remain consistent within
their uncertainties. The $\chi^{2}/\text{ndp}|^{p}_{x_{F},\sqrt{s}}$ have a
mild dependence on the pion PDFs, only through the correlation of LDMEs and
PDFs in the global fit.

In NRQCD, the relative weighting between $q\bar{q}$ and $GG$ subprocesses is
set by a convolution of the pQCD partonic cross sections, the associated
parton densities, and the LDMEs. The $F$ factor does not modify the shape of
$d\sigma/dx_{F}$. Therefore, adequate shapes of $d\sigma/dx_{F}$ distributions
of individual $GG$ and $q\bar{q}$ contributions from NRQCD calculations are
required to achieve a reasonable description of the data, particularly for
$x_{F}>0.5$. Since the partonic cross sections and the nucleon PDFs involved
in the calculations of the cross sections remain the same, the variation of
the results originates from the difference in the pion PDFs and the LDMEs.

Data | SMRS | GRV | JAM | xFitter
---|---|---|---|---
Exp ($P_{beam}$) |  $\chi^{2}$/ndp | $F$ |  $\chi^{2}$/ndp | $F$ |  $\chi^{2}$/ndp | $F$ |  $\chi^{2}$/ndp | $F$
E672, E706 ($\sigma^{J/\psi}$) | 1.3 | 0.80 $\pm$ 0.01 | 2.6 | 0.79 $\pm$ 0.01 | 8.8 | 1.21 $\pm$ 0.02 | 4.2 | 1.08 $\pm$ 0.02
E705 ($\sigma^{J/\psi}$) | 2.0 | 0.98 $\pm$ 0.02 | 1.7 | 0.96 $\pm$ 0.02 | 6.9 | 1.19 $\pm$ 0.00 | 2.6 | 1.18 $\pm$ 0.01
NA3 ($\sigma^{J/\psi}$) | 2.1 | 0.86 $\pm$ 0.02 | 2.3 | 0.87 $\pm$ 0.02 | 2.2 | 1.02 $\pm$ 0.02 | 2.9 | 1.01 $\pm$ 0.02
NA3 ($\sigma^{J/\psi}$) | 1.3 | 0.87 $\pm$ 0.02 | 0.9 | 0.89 $\pm$ 0.02 | 2.3 | 0.92 $\pm$ 0.02 | 1.5 | 0.95 $\pm$ 0.02
WA11 ($\sigma^{J/\psi}$) | 3.7 | 1.02 $\pm$ 0.02 | 8.5 | 1.02 $\pm$ 0.02 | 37.2 | 1.10 $\pm$ 0.02 | 22.0 | 1.12 $\pm$ 0.02
NA3 ($\sigma^{J/\psi}$) | 1.6 | 1.24 $\pm$ 0.03 | 1.3 | 1.23 $\pm$ 0.03 | 1.3 | 1.07 $\pm$ 0.03 | 1.6 | 1.18 $\pm$ 0.03
E537 ($\sigma^{J/\psi}$) | 3.3 | 0.88 $\pm$ 0.00 | 1.6 | 0.88 $\pm$ 0.01 | 2.8 | 0.88 $\pm$ 0.00 | 2.1 | 0.88 $\pm$ 0.01
WA39 ($\sigma^{J/\psi}$) | 1.4 | 1.30 $\pm$ 0.04 | 1.4 | 1.18 $\pm$ 0.07 | 5.1 | 0.70 $\pm$ 0.00 | 1.3 | 0.70 $\pm$ 0.05
E672, E706 ($\sigma^{\psi(2S)}$) | 0.2 | 0.80 $\pm$ 0.01 | 0.2 | 0.79 $\pm$ 0.01 | 0.4 | 1.21 $\pm$ 0.02 | 0.2 | 1.08 $\pm$ 0.02
E615 ($\sigma^{\psi(2S)}/\sigma^{J/\psi}$) | 1.6 | 1 $\pm$ 0 | 1.7 | 1 $\pm$ 0 | 5.6 | 1 $\pm$ 0 | 4.3 | 1 $\pm$ 0
HERA-B ($\sigma^{\psi(2S)}/\sigma^{J/\psi}$) | 1.4 | 1 $\pm$ 0 | 1.5 | 1 $\pm$ 0 | 1.2 | 1 $\pm$ 0 | 1.2 | 1 $\pm$ 0
NA50 ($\sigma^{\psi(2S)}/\sigma^{J/\psi}$) | 1.0 | 1 $\pm$ 0 | 1.6 | 1 $\pm$ 0 | 1.2 | 1 $\pm$ 0 | 1.1 | 1 $\pm$ 0
E789 ($\sigma^{\psi(2S)}/\sigma^{J/\psi}$) | 3.1 | 1 $\pm$ 0 | 3.3 | 1 $\pm$ 0 | 2.8 | 1 $\pm$ 0 | 2.9 | 1 $\pm$ 0
E771 ($\sigma^{\psi(2S)}/\sigma^{J/\psi}$) | 0.3 | 1 $\pm$ 0 | 0.3 | 1 $\pm$ 0 | 0.3 | 1 $\pm$ 0 | 0.3 | 1 $\pm$ 0
E705 ($\sigma^{J/\psi}$) | 2.3 | 1.20 $\pm$ 0.00 | 2.2 | 1.20 $\pm$ 0.00 | 4.0 | 1.20 $\pm$ 0.00 | 3.1 | 1.20 $\pm$ 0.00
NA3 ($\sigma^{J/\psi}$) | 1.0 | 1.00 $\pm$ 0.01 | 1.2 | 1.00 $\pm$ 0.01 | 2.9 | 0.99 $\pm$ 0.01 | 1.6 | 1.00 $\pm$ 0.01
Table 4: Results of the NRQCD fits. The columns display the $\chi^{2}$/ndp
values and fitted normalization factors $F$ for each of the selected data sets
and for SMRS, GRV, JAM and xFitter pion PDFs.

###  V.3 Differential cross sections for $J/\psi$

Figure 3: Differential cross sections for $J/\psi$ production with a
39.5-GeV/$c$ $\pi^{-}$ beam Corden  _et al._ (1981). The data are compared to
the NRQCD fit results for the SMRS, GRV, xFitter, and JAM PDFs. The total
cross sections and $q\bar{q}$, $GG$, and $qG$ contributions are denoted as
solid black, dashed blue, dotted red, and dot-dashed green lines,
respectively. The uncertainty bands associated with JAM and xFitter PDFs are
also shown. Figure 4: Same as Fig. 4 for $J/\psi$ production data with a a 515
GeV/$c$ $\pi^{-}$ beam Gribushin  _et al._ (1996).

A comparison of the $J/\psi$ production data and the NRQCD calculations in
terms of the subprocess contributions has been made for all of the data sets
included in the fit. Irrespective of the pion PDFs, the relative weighting of
$q\bar{q}$ and $GG$ shows a strong energy dependence. At the lowest energy,
the $q\bar{q}$ term provides the major contribution to the cross section,
similar to the DY production, while the $GG$ contribution is dominant at the
highest beam energies. A global analysis of charmonium data sets with a wide
range of beam energy could simultaneously constrain both pion’s valence quark
and gluon distributions. It is instructive to compare the results obtained
with each of the four pion PDFs. This comparison is illustrated in Figs. 3 and
4 for the data with pion beam momenta of 39.5 GeV/$c$ Corden  _et al._ (1981)
and 515 GeV/$c$ Gribushin  _et al._ (1996). The $\chi^{2}$/ndf values are
displayed in the plots.

At the lowest beam momentum of 39.5 GeV/$c$ (Fig. 3), the $q\bar{q}$
subprocess provides the largest contribution to the cross section over the
whole $x_{F}$ region. The $GG$ contribution is much reduced, so that the shape
of the $x_{F}$ distribution is essentially determined by the shape of the
$q\bar{q}$ contribution. Since the pion valence-quark distribution is well
determined from the DY data, good $\chi^{2}/\text{ndf}$ values are obtained
for the four PDFs. Nevertheless, the agreement with the data is less
satisfactory for JAM. Figure 3 also suggests that future $J/\Psi$ data at
negative $x_{F}$ with low beam energies could further constrain the pion
valence-quark distribution at lower $x$.

At the highest beam momentum of 515 GeV/$c$, where the $GG$ contribution
becomes dominant, Fig. 4 shows that SMRS and GRV are favored over JAM and
xFitter. The fraction of the $GG$ component is maximized around $x_{F}=0$,
corresponding to the gluon distribution $G_{\pi}(x)$ around $x\sim$0.1–0.2. As
a result of the rapid drop of the $G_{\pi}(x)$ toward $x=1$, the $GG$
contribution quickly decreases at large $x_{F}$. In contrast, the $q\bar{q}$
contribution has a slower fall-off toward high $x_{F}$ because of a relatively
strong pion valence antiquark density, in comparison with the gluon one, at
large $x$. Consequently, the $q\bar{q}$ contribution has a broader $x_{F}$
distribution than that of the $GG$ contribution and the relative importance of
$q\bar{q}$ rises at the large $x_{F}$ region. The ratio of $q\bar{q}$ to $GG$
shows a strong $x_{F}$ dependence, making the $x_{F}$-differential cross
sections at high energies particularly sensitive to the shape of pion
$G_{\pi}(x)$.

Similar conclusions can be drawn for the intermediate energies used in this
analysis. The corresponding figures are available in the supplemental material
Sup (2022). As a general observation, the $q\bar{q}$ and $GG$ contributions
have quite similar strengths for the fits with SMRS and GRV, whereas the
$q\bar{q}$ contribution is the dominant component for the fits with JAM and
xFitter. In terms of $\chi^{2}$/ndf, the data show a slight preference for GRV
and SMRS.

###  V.4 Differential cross sections for $\psi(2S)$ and the $R_{\psi}(x_{F})$
ratios

Figure 5: Same as Fig. 4 for $\psi(2S)$ production with a 515-GeV/$c$
$\pi^{-}$ beam Gribushin  _et al._ (1996).

Additional information on the charmonium production mechanism can be obtained
by comparing the production of the two charmonium states, $J/\psi$ and
$\psi(2S)$. The strengths of their $q\bar{q}$ and $GG$ subprocesses are
controlled by the associated LDMEs. In comparison with the $J/\psi$, the
smaller cross section for the $\psi(2S)$ production implies also smaller
LDMEs. The fitted LDMEs are indeed smaller, but interestingly, not in the same
proportion. As shown in Table 3, the values of the
$\langle\mathcal{O}_{8}^{\psi(2S)}[^{3}S_{1}]\rangle$ LDMEs for $\psi(2S)$ are
smaller than that for $J/\psi$ by about a factor of two. In contrast, the
$\Delta_{8}^{\psi(2S)}$ values for $\psi(2S)$ are an order of magnitude
smaller. This is illustrated in Fig. 5 for the E672/E706 $\psi(2S)$ data taken
at 515 GeV/$c$ Gribushin  _et al._ (1996). In comparison with the production
of $J/\psi$ at the same energy (Fig. 4), the $q\bar{q}$ contribution is
greatly enhanced in $\psi(2S)$ production. Figure 5 shows that this
observation is valid for all pion PDFs, and the $q\bar{q}$ term is even
dominant for JAM and xFitter. For the fit with the SMRS pion PDFs around
$x_{F}=0$, the $q\bar{q}$ component accounts for about 15% of the direct part
(feed-down excluded) of the $J/\psi$ cross section. Its fraction rises to
nearly 30% for $\psi(2S)$. Obviously, the increase of the $q\bar{q}$ term is
compensated by a decrease of the $GG$ term. This significant difference
between the two charmonium states can only be partially explained by the
larger $\psi(2S)$ mass. Its full understanding would require further
investigations.

Figure 6: The $\psi(2S)$ to $J/\psi$ cross section ratios $R_{\psi}(x_{F})$
for $J/\psi$ and $\psi(2S)$ production with a 252-GeV/$c$ $\pi^{-}$ beam
Heinrich  _et al._ (1991). The data are compared to the NRQCD fit results for
the SMRS, GRV, xFitter, and JAM PDFs. The ratios of total cross sections and
individual $R^{q\bar{q}}_{\psi}(x_{F})$ and $R^{GG}_{\psi}(x_{F})$
contributions are denoted as solid black, dashed blue, and dotted red lines,
respectively.

The observations above are consistent with the measurements of the $\psi(2S)$
to $J/\psi$ ratios, $R_{\psi}(x_{F})$. The largest statistics on
$R_{\psi}(x_{F})$ have been collected by the E615 experiment for an incident
pion momentum of 252 GeV/$c$ Heinrich  _et al._ (1991). The data are compared
to the NRQCD fits with each of the four pion PDFs in Fig. 6. The
$R_{\psi}(x_{F})$ shows a strong $x_{F}$ dependence and this suggests that the
relative weights of the individual subprocesses $q\bar{q}$ and $GG$ components
in $J/\psi$ and $\psi(2S)$ production are distinctly different. We note that
the CEM models predicts an $x_{F}$-independent $R_{\psi}(x_{F})$ Abt  _et al._
(2007), since the fractions of $q\bar{q}$ and $GG$ components are identical
for each charmonium state. In NRQCD, an $x_{F}$-dependent $R_{\psi}(x_{F})$ is
possible due to the different LDMEs associated with the $q\bar{q}$ and $GG$
channels in producing $J/\psi$ and $\psi(2S)$. The pronounced $x_{F}$
dependence of $R_{\psi}(x_{F})$ in Fig. 6 clearly disfavors the CEM model.

As shown in Figs. 4 and 5, the $q\bar{q}$ subprocess gives a significantly
broader $x_{F}$ distribution than the $GG$ subprocess. This is caused by the
slower fall-off of the valence-quark distribution than the gluon distribution
toward $x=1$. Therefore, the pronounced rise in the $R_{\psi}(x_{F})$ data at
forward $x_{F}$, shown in Fig. 6, clearly indicates that the $q\bar{q}$
subprocess is more important for the $\psi(2S)$ production than for the
$J/\psi$ production.

It is also instructive to examine the $x_{F}$ dependence of $R_{\psi}(x_{F})$
from the $q\bar{q}$ and $GG$ subprocesses separately. In Fig. 6, the dashed
blue and dotted red curves correspond, respectively, to

$R^{q\bar{q}}_{\psi}(x_{F})\equiv\frac{\sigma^{q\bar{q}}_{\psi(2S)}(x_{F})}{\sigma^{q\bar{q}}_{J/\psi}(x_{F})};\leavevmode\nobreak\
\leavevmode\nobreak\
R^{GG}_{\psi}(x_{F})\equiv\frac{\sigma^{GG}_{\psi(2S)}(x_{F})}{\sigma^{GG}_{J/\psi}(x_{F})},$
(6)

where the superscripts $q\bar{q}$ and $GG$ denote the two subprocesses.
Neglecting the tiny contribution from the $qG$ subprocess, one can then obtain

$\displaystyle R_{\psi}(x_{F})\equiv$
$\displaystyle\frac{\sigma_{\psi(2S)}(x_{F})}{\sigma_{J/\psi}(x_{F})}$ (7)
$\displaystyle=$
$\displaystyle[A(x_{F})R^{q\bar{q}}_{\psi}(x_{F})+B(x_{F})R^{GG}_{\psi}(x_{F})],$

where

$A(x_{F})=\frac{\sigma^{q\bar{q}}_{J/\psi}(x_{F})}{\sigma_{J/\psi}(x_{F})}\mbox{\leavevmode\nobreak\
\leavevmode\nobreak\ and\leavevmode\nobreak\ \leavevmode\nobreak\
}B(x_{F})=\frac{\sigma^{GG}_{J/\psi}(x_{F})}{\sigma_{J/\psi}(x_{F})}\\\ $

have the property $0\leq A(x_{F})\leq 1$ and $0\leq B(x_{F})\leq 1$. It
follows that $R_{\psi}(x_{F})$ must be bounded by $R^{GG}_{\psi}(x_{F})$ and
$R^{q\bar{q}}_{\psi}(x_{F})$ in Fig. 6. As shown in Fig. 6, the
$R_{\psi}(x_{F})$ data largely fall within these two bounds for calculations
with the SMRS and GRV PDFs, while a large fraction of the data are outside of
these bounds for the calculations using the JAM and xFitter PDFs. The striking
contrast between the SMRS/GRV and the JAM/xFitter PDFs in their ability to
describe the $R_{\psi}(x_{F})$ data in Fig. 6 illustrates the advantages of
the $R_{\psi}(x_{F})$ data in constraining the pion PDFs. We also note that
none of the pion PDFs can explain the sharp rise of the $R_{\psi}(x_{F})$ data
beyond $x_{F}=0.8$. This incompatibility at large $x_{F}$ could be due to
either higher-twist effects Heinrich  _et al._ (1991) or higher-order QCD
processes that are beyond the present leading-order NRQCD analysis.

Figure 7: Differential cross sections for $J/\psi$ production with a 300
GeV/$c$ $\pi^{-}$ beam Antoniazzi  _et al._ (1992). The data are compared to
the fit results with SMRS, GRV, xFitter, and JAM PDFs. The total cross section
and its decomposition into contributions from CS and CO subprocesses are
denoted as solid black, dashed blue and dotted red lines, respectively.

Our analysis also shows that fixed-target charmonium production data are
particularly sensitive to the color octet contribution to the cross section.
This is illustrated in Fig. 7 which displays the decomposition of the $J/\psi$
$x_{F}$-dependent cross sections from the E705 experiment Antoniazzi  _et al._
(1992) into color octet and color singlet contributions. The CO contribution
plays a dominant role in the $J/\psi$ production across the entire $x_{F}$
range, and this observation is valid for any of the four pion PDFs.

###  V.5 Integrated cross sections

Because of the presence of valence antiquarks in the pion, the $q\bar{q}$ and
$GG$ subprocesses to the $J/\psi$ production with proton and pion beams have
different contributions to the integrated cross sections. In the production
with a proton beam the $GG$ contribution is dominant across all center-of-mass
energies $\sqrt{s}$ except near threshold. With pion beams the $q\bar{q}$
contribution is significantly enhanced. It dominates at low energies, with
$GG$ contribution gradually becoming important as $\sqrt{s}$ increases.

Our analysis of the differential cross sections shows that the relative
contributions of the $q\bar{q}$ and $GG$ subprocesses in the production of
$J/\psi$ and $\psi(2S)$ differ considerably. The same conclusion can be drawn
from the integrated cross sections. Figure 8 shows the comparison of data and
NRQCD calculations for the $J/\psi$ and $\psi(2S)$ production cross sections
in $pN$ and $\pi^{-}N$ collisions with the SMRS pion PDFs and the “FIT” LDMEs
in Table 3. The fractions of $q\bar{q}$ and $GG$ contributions as a function
of $\sqrt{s}$ vary considerably, reflecting the differences of the
corresponding gluon and quark parton distributions between the pion PDFs. For
SMRS, whose gluon strength at large $x$ is relatively strong, the $GG$
contribution starts to dominate the cross section beyond $\sqrt{s}=$18 and 10
GeV for the production of $\psi(2S)$ and $J/\psi$, respectively, while the
transition happens at larger $\sqrt{s}$ for the results with JAM, in
consequence of a relatively weak gluon strength. The plots for GRV, JAM and
xFitter pion PDFs are provided in the supplemental material Sup (2022). All
these observations confirm our previous conclusion: the $q\bar{q}$
contribution plays a much more important role in the $\psi(2S)$ production,
compared to $J/\psi$.

Figure 8: Integrated charmonium cross sections in $pN$ and $\pi^{-}N$
collisions. The data for $J/\psi$ and $\psi(2S)$ production are compared to
the results of NRQCD calculations with the SMRS pion PDFs and the “FIT” LDMEs
in Table 3. The total cross section and its $q\bar{q}$ and $GG$ contributions
are denoted as solid black, dashed blue and dotted red lines, respectively.

###  V.6 Systematic studies

So far only the uncertainties associated with the parametrizations of JAM and
xFitter PDFs have been taken into account. Our results are also sensitive to
the NRQCD input parameters and to the choice of the nuclear PDFs. We have
checked that fits performed with the nCTEQ15 Kovarik  _et al._ (2016)
parametrization instead of EPPS16 result in negligible differences. Fits with
the factorization/renormalization scale parameter $\mu_{R}$ set to 1, 2, and
3$m_{c}$, with $m_{c}$ = 1.5 GeV/$c^{2}$, have also been made. The values of
the total $\chi^{2}$/ndf do not vary much: they remain nearly unchanged
between $\mu_{R}=m_{c}$ and $\mu_{R}=3m_{c}$. The effect on the values of the
LDMEs is more important. For both $J/\psi$ and $\psi(2S)$ the fitted LDMEs
increase by nearly a factor of two when $\mu_{R}$ increases from
$\mu_{R}=m_{c}$ to $\mu_{R}=3m_{c}$. Nevertheless, the shape and the magnitude
of the final cross section remain nearly unchanged, as illustrated in Fig. 9
for the fit with the SMRS PDFs. The relative contributions of the $q\bar{q}$
and $GG$ subprocesses for the three values of $\mu_{R}$ are only slightly
modified. Changing the charm quark mass from 1.5 to 1.2 GeV/$c^{2}$ produces a
much stronger effect. The total $\chi^{2}$/ndf increases by up to an order of
magnitude, depending on the pion PDF and the $\mu_{R}$ value used. The overall
quality of the fits is strongly degraded and the LDMEs reach unreasonable
levels. The systematic studies with GRV, JAM and xFitter pion PDFs lead to
results fully consistent with these conclusions. The corresponding figures are
available in the supplemental material Sup (2022).

Figure 9: The NRQCD results with variation of charm quark mass $m_{c}$ and
renormalization scale $\mu_{R}$, compared with the $d\sigma/dx_{F}$ data of
$J/\psi$ production off the beryllium target with a 515-GeV/$c$ $\pi^{-}$ beam
from the E672/E706 experiment Gribushin  _et al._ (1996). The pion PDFs used
for the calculation is SMRS. The total cross sections and $q\bar{q}$, $GG$,
and $qG$ contributions are denoted as solid black, dashed blue, dotted red and
dot-dashed green lines, respectively. The charm quark mass $m_{c}$,
factorization scale $\mu_{F}$, and renormalization scale $\mu_{R}$ used for
the NRQCD calculation as well as the fit $\chi^{2}$/ndf are displayed in each
plot.

##  VI Discussion

The values of the color-octet LDMEs, resulting from the fits to the data may
contain model uncertainties, although they provide a good description of the
data. The formalism used is limited to LO and is able to determine
individually the CO $\langle\mathcal{O}_{8}^{H}[^{3}S_{1}]\rangle$ LDMEs for
$J/\psi$ and $\psi(2S)$ only. The $\Delta_{8}^{J/\psi}$ and
$\Delta_{8}^{\psi(2S)}$ terms combine each three additional color octet LDMEs.
Furthermore, most of the data included in the analysis have transverse momenta
$p_{T}$ smaller than 3 GeV/$c$. This is in sharp contrast with most of the
available LDMEs that result from fits at much larger energies and for
transverse momenta $p_{T}$ usually larger than 5 GeV/$c$. Therefore, the
comparison with these LDMEs remains qualitative.

For the fits on the $J/\psi$ data the
$\langle\mathcal{O}_{8}^{J/\psi}[^{3}S_{1}]\rangle$ values obtained for SMRS
and GRV are nearly an order of magnitude larger than most of the published
LDMEs Butenschoen and Kniehl (2011); Chao  _et al._ (2012). Within the
experimental uncertainties, they are quite close to what were reported in
Refs. Zhang  _et al._ (2015); Bodwin  _et al._ (2014) derived from the data on
$\eta_{c}$ production using spin symmetry relations. Similarly for $\psi(2S)$,
the fitted $\langle\mathcal{O}_{8}^{\psi(2S)}[^{3}S_{1}]\rangle$ values are
nearly an order of magnitude larger than the values quoted in Refs. Ma  _et
al._ (2011); Gong  _et al._ (2013) and more recently in Ref. Butenschoen and
Kniehl (2022). Within the systematic uncertainties associated with the fits
and given the exploratory status of the present analysis, the results obtained
are considered satisfactory.

Our analysis is performed using a leading-order NRQCD framework only. The
results obtained may vary if a more advanced NRQCD formalism with higher order
terms is applied. In addition, for most of the fixed-target data considered
here, the mean transverse momenta are smaller than the $J/\psi$ mass.
Inclusion of higher-order corrections could therefore provide a better
description, but probably would not change the general conclusions. The
analysis has been also limited to data taken with only light targets. A large
amount of data of $x_{F}$-differential cross sections with heavier targets
have been collected in the past. These data could be included in a more
complete analysis if the energy loss effects Arleo  _et al._ (2019)
responsible for the suppression of the charmonium cross section in hadron-
nucleus collisions are reliably accounted for.

##  VII Conclusion

We have analyzed fixed-target experimental cross sections for $J/\psi$ and
$\psi(2S)$ production using the NRQCD framework. To minimize nuclear matter
effects, only data on hydrogen, lithium and beryllium targets were selected.
Heavier targets were only considered for the data on the $J/\psi$ to
$\psi(2S)$ ratios. Assuming the universality of the NRQCD approach, both pion
and proton-induced data sets were included in the analysis. Fits to the
individual $x_{F}$-differential cross sections and their ratios have been
made, using four different pion PDF parametrizations. The proton data,
although not directly sensitive to the pion PDFs, enrich the selection and
contribute to the stability of the final results.

A simultaneous fit to all pion and proton data sets has been achieved. The
results of these common fits show that the relative fractions of the
$q\bar{q}$ and $GG$ contributions to the cross sections strongly depend on the
beam particle, on its incoming energy and on the $x_{F}$ region considered. A
strong dependence on the pion PDF parametrization used is observed and
particularly on the magnitude of the pion gluon distribution. The results
indicate a clear preference for parametrizations with larger gluon
distributions at relatively large $x$. Good agreement with the data is
obtained with the SMRS and GRV PDFs. The fits with the recent JAM and xFitter
parametrizations turn out to show much larger deviations for most of the data
sets.

The comparison between the results for $J/\psi$ and $\psi(2S)$ production
leads to an important new observation: the strengths of the $q\bar{q}$ and
$GG$ contributions to these two charmonium states are - unexpectedly - quite
different. The $q\bar{q}$ component of the $\psi(2S)$ cross section is,
proportionally, few times larger than the $q\bar{q}$ component of the $J/\psi$
cross section. This interesting feature is confirmed for both differential and
integrated cross sections and for both pion and proton beams. The production
of $\psi(2S)$ appears to be more sensitive to the pion’s valence quark
distribution than that of $J/\psi$. This observation could be relevant for a
better understanding of the charmonium production mechanism.

Within the NRQCD framework our analysis shows that the fixed-target pion-
induced $J/\psi$ and $\psi(2S)$ data provide strong constraints of the pion
gluon density, particularly in the large-$x$ region. This conclusion is fully
supported by the results from our previous study done with the Color
Evaporation Model Chang  _et al._ (2020). While further theoretical efforts
are required to reduce the model dependence in describing the charmonium
production, we argue that, an extension of the global analysis of pion PDFs to
include charmonium production data in the future is warranted Bourrely  _et
al._ (2022).

New results of Drell-Yan as well as $J/\psi$ measurements in $\pi A$ reactions
will be available from the CERN COMPASS Aghasyan  _et al._ (2017) and AMBER
Adams  _et al._ (2018) experiments in the near future. These data will be
important in providing better knowledge of the pion PDFs. For the longer-term
electron-ion collider projects in U.S. and China, the pion as well kaon
structures are planned to be explored using the tagged DIS process Aguilar
_et al._ (2019); Anderle  _et al._ (2021); Chávez  _et al._ (2022).

## Acknowledgments

This work was supported in part by the U.S. National Science Foundation and
National Science and Technology Council of Taiwan (R.O.C.).

## References

  * Horn and Roberts (2016) T. Horn and C. D. Roberts, J. Phys. G 43, 073001 (2016), arXiv:1602.04016 [nucl-th] .
  * Nam (2012) S.-i. Nam, Phys. Rev. D 86, 074005 (2012), arXiv:1205.4156 [hep-ph] .
  * Watanabe  _et al._ (2016) A. Watanabe, C. W. Kao, and K. Suzuki, Phys. Rev. D 94, 114008 (2016), arXiv:1610.08817 [hep-ph] .
  * Watanabe  _et al._ (2018) A. Watanabe, T. Sawada, and C. W. Kao, Phys. Rev. D 97, 074015 (2018), arXiv:1710.09529 [hep-ph] .
  * Hutauruk  _et al._ (2016) P. T. P. Hutauruk, I. C. Cloet, and A. W. Thomas, Phys. Rev. C 94, 035201 (2016), arXiv:1604.02853 [nucl-th] .
  * Lan  _et al._ (2019) J. Lan, C. Mondal, S. Jia, X. Zhao, and J. P. Vary, Phys. Rev. Lett. 122, 172001 (2019), arXiv:1901.11430 [nucl-th] .
  * Lan  _et al._ (2020) J. Lan, C. Mondal, S. Jia, X. Zhao, and J. P. Vary, Phys. Rev. D 101, 034024 (2020), arXiv:1907.01509 [nucl-th] .
  * Lan  _et al._ (2022) J. Lan, K. Fu, C. Mondal, X. Zhao, and j. P. Vary (BLFQ), Phys. Lett. B 825, 136890 (2022), arXiv:2106.04954 [hep-ph] .
  * de Teramond  _et al._ (2018) G. F. de Teramond, T. Liu, R. S. Sufian, H. G. Dosch, S. J. Brodsky, and A. Deur (HLFHS), Phys. Rev. Lett. 120, 182001 (2018), arXiv:1801.09154 [hep-ph] .
  * Watanabe  _et al._ (2020) A. Watanabe, T. Sawada, and M. Huang, Phys. Lett. B 805, 135470 (2020), arXiv:1910.10008 [hep-ph] .
  * Lan and Mondal (2020) J. Lan and C. Mondal, Phys. Lett. B 807, 135613 (2020), arXiv:2007.05858 [hep-ph] .
  * Han  _et al._ (2020) C. Han, H. Xing, X. Wang, Q. Fu, R. Wang, and X. Chen, Phys. Lett. B 800, 135066 (2020), arXiv:1809.01549 [hep-ph] .
  * Han  _et al._ (2021) C. Han, G. Xie, R. Wang, and X. Chen, Eur. Phys. J. C 81, 302 (2021), arXiv:2010.14284 [hep-ph] .
  * Chang  _et al._ (2014) L. Chang, C. Mezrag, H. Moutarde, C. D. Roberts, J. Rodríguez-Quintero, and P. C. Tandy, Phys. Lett. B 737, 23 (2014), arXiv:1406.5450 [nucl-th] .
  * Chang and Thomas (2015) L. Chang and A. W. Thomas, Phys. Lett. B 749, 547 (2015), arXiv:1410.8250 [nucl-th] .
  * Chen  _et al._ (2016) C. Chen, L. Chang, C. D. Roberts, S. Wan, and H.-S. Zong, Phys. Rev. D 93, 074021 (2016), arXiv:1602.01502 [nucl-th] .
  * Shi  _et al._ (2018) C. Shi, C. Mezrag, and H.-s. Zong, Phys. Rev. D 98, 054029 (2018), arXiv:1806.10232 [nucl-th] .
  * Bednar  _et al._ (2020) K. D. Bednar, I. C. Cloët, and P. C. Tandy, Phys. Rev. Lett. 124, 042002 (2020), arXiv:1811.12310 [nucl-th] .
  * Ding  _et al._ (2020) M. Ding, K. Raya, D. Binosi, L. Chang, C. D. Roberts, and S. M. Schmidt, Phys. Rev. D 101, 054014 (2020), arXiv:1905.05208 [nucl-th] .
  * Cui  _et al._ (2020) Z.-F. Cui, M. Ding, F. Gao, K. Raya, D. Binosi, L. Chang, C. D. Roberts, J. Rodríguez-Quintero, and S. M. Schmidt, Eur. Phys. J. C 80, 1064 (2020).
  * Freese  _et al._ (2021) A. Freese, I. C. Cloët, and P. C. Tandy, Phys. Lett. B 823, 136719 (2021), arXiv:2103.05839 [hep-ph] .
  * Chang and Roberts (2021) L. Chang and C. D. Roberts, Chin. Phys. Lett. 38, 081101 (2021), arXiv:2106.08451 [hep-ph] .
  * Cui  _et al._ (2022a) Z. F. Cui, M. Ding, J. M. Morgado, K. Raya, D. Binosi, L. Chang, J. Papavassiliou, C. D. Roberts, J. Rodríguez-Quintero, and S. M. Schmidt, Eur. Phys. J. A 58, 10 (2022a), arXiv:2112.09210 [hep-ph] .
  * Cui  _et al._ (2022b) Z. F. Cui, M. Ding, J. M. Morgado, K. Raya, D. Binosi, L. Chang, F. De Soto, C. D. Roberts, J. Rodríguez-Quintero, and S. M. Schmidt, Phys. Rev. D 105, L091502 (2022b), arXiv:2201.00884 [hep-ph] .
  * Zhang  _et al._ (2019) J.-H. Zhang, J.-W. Chen, L. Jin, H.-W. Lin, A. Schäfer, and Y. Zhao, Phys. Rev. D 100, 034505 (2019), arXiv:1804.01483 [hep-lat] .
  * Sufian  _et al._ (2019) R. S. Sufian, J. Karpie, C. Egerer, K. Orginos, J.-W. Qiu, and D. G. Richards, Phys. Rev. D 99, 074507 (2019), arXiv:1901.03921 [hep-lat] .
  * Izubuchi  _et al._ (2019) T. Izubuchi, L. Jin, C. Kallidonis, N. Karthik, S. Mukherjee, P. Petreczky, C. Shugert, and S. Syritsyn, Phys. Rev. D 100, 034516 (2019), arXiv:1905.06349 [hep-lat] .
  * Joó  _et al._ (2019) B. Joó, J. Karpie, K. Orginos, A. V. Radyushkin, D. G. Richards, R. S. Sufian, and S. Zafeiropoulos, Phys. Rev. D 100, 114512 (2019), arXiv:1909.08517 [hep-lat] .
  * Sufian  _et al._ (2020) R. S. Sufian, C. Egerer, J. Karpie, R. G. Edwards, B. Joó, Y.-Q. Ma, K. Orginos, J.-W. Qiu, and D. G. Richards, Phys. Rev. D 102, 054508 (2020), arXiv:2001.04960 [hep-lat] .
  * Chen  _et al._ (2020) J.-W. Chen, H.-W. Lin, and J.-H. Zhang, Nucl. Phys. B 952, 114940 (2020), arXiv:1904.12376 [hep-lat] .
  * Gao  _et al._ (2020) X. Gao, L. Jin, C. Kallidonis, N. Karthik, S. Mukherjee, P. Petreczky, C. Shugert, S. Syritsyn, and Y. Zhao, Phys. Rev. D 102, 094513 (2020), arXiv:2007.06590 [hep-lat] .
  * Alexandrou  _et al._ (2021a) C. Alexandrou, S. Bacchio, I. Cloet, M. Constantinou, K. Hadjiyiannakou, G. Koutsou, and C. Lauer (ETM), Phys. Rev. D 103, 014508 (2021a), arXiv:2010.03495 [hep-lat] .
  * Alexandrou  _et al._ (2021b) C. Alexandrou, S. Bacchio, I. Cloët, M. Constantinou, K. Hadjiyiannakou, G. Koutsou, and C. Lauer (ETM), Phys. Rev. D 104, 054504 (2021b), arXiv:2104.02247 [hep-lat] .
  * Fan and Lin (2021) Z. Fan and H.-W. Lin, Phys. Lett. B 823, 136778 (2021), arXiv:2104.06372 [hep-lat] .
  * Detmold  _et al._ (2022) W. Detmold, A. V. Grebe, I. Kanamori, C. J. D. Lin, S. Mondal, R. J. Perry, and Y. Zhao (HOPE), Phys. Rev. D 105, 034506 (2022), arXiv:2109.15241 [hep-lat] .
  * Gao  _et al._ (2022) X. Gao, A. D. Hanlon, N. Karthik, S. Mukherjee, P. Petreczky, P. Scior, S. Shi, S. Syritsyn, Y. Zhao, and K. Zhou, (2022), arXiv:2208.02297 [hep-lat] .
  * Chang and Dutta (2013) W.-C. Chang and D. Dutta, Int. J. Mod. Phys. E 22, 1330020 (2013), arXiv:1306.3971 [nucl-th] .
  * Londergan  _et al._ (1995) J. T. Londergan, G. Q. Liu, E. N. Rodionov, and A. W. Thomas, Phys. Lett. B 361, 110 (1995).
  * Bonesini  _et al._ (1988) M. Bonesini _et al._ (WA70), Z. Phys. C 37, 535 (1988).
  * Khoze  _et al._ (2006) V. A. Khoze, A. D. Martin, and M. G. Ryskin, Eur. Phys. J. C 48, 797 (2006), arXiv:hep-ph/0606213 .
  * McKenney  _et al._ (2016) J. R. McKenney, N. Sato, W. Melnitchouk, and C.-R. Ji, Phys. Rev. D 93, 054011 (2016), arXiv:1512.04459 [hep-ph] .
  * Gluck  _et al._ (1978) M. Gluck, J. F. Owens, and E. Reya, Phys. Rev. D 17, 2324 (1978).
  * Barger  _et al._ (1980) V. D. Barger, W.-Y. Keung, and R. J. N. Phillips, Z. Phys. C 6, 169 (1980).
  * Owens (1984) J. F. Owens, Phys. Rev. D 30, 943 (1984).
  * Einhorn and Ellis (1975) M. B. Einhorn and S. D. Ellis, Phys. Rev. D 12, 2007 (1975).
  * Fritzsch (1977) H. Fritzsch, Phys. Lett. B 67, 217 (1977).
  * Halzen (1977) F. Halzen, Phys. Lett. B 69, 105 (1977).
  * Bodwin  _et al._ (2005) G. T. Bodwin, E. Braaten, and J. Lee, Phys. Rev. D 72, 014004 (2005), arXiv:hep-ph/0504014 .
  * Bodwin  _et al._ (1995) G. T. Bodwin, E. Braaten, and G. P. Lepage, Phys. Rev. D 51, 1125 (1995), [Erratum: Phys.Rev.D 55, 5853 (1997)], arXiv:hep-ph/9407339 .
  * Schuler (1994) G. A. Schuler,  _Quarkonium production and decays_ , Ph.D. thesis, Hamburg U. (1994), arXiv:hep-ph/9403387 .
  * Vogt (1999) R. Vogt, Phys. Rept. 310, 197 (1999).
  * Chang  _et al._ (2020) W.-C. Chang, J.-C. Peng, S. Platchkov, and T. Sawada, Phys. Rev. D 102, 054024 (2020), arXiv:2006.06947 [hep-ph] .
  * Hsieh  _et al._ (2021) C.-Y. Hsieh, Y.-S. Lian, W.-C. Chang, J.-C. Peng, S. Platchkov, and T. Sawada, Chin. J. Phys. 73, 13 (2021), arXiv:2103.11660 [hep-ph] .
  * Beneke and Rothstein (1996) M. Beneke and I. Z. Rothstein, Phys. Rev. D 54, 2005 (1996), [Erratum: Phys.Rev.D 54, 7082 (1996)], arXiv:hep-ph/9603400 .
  * Barry  _et al._ (2018) P. C. Barry, N. Sato, W. Melnitchouk, and C.-R. Ji, Phys. Rev. Lett. 121, 152001 (2018), arXiv:1804.01965 [hep-ph] .
  * Cao  _et al._ (2021) N. Y. Cao, P. C. Barry, N. Sato, and W. Melnitchouk (Jefferson Lab Angular Momentum), Phys. Rev. D 103, 114014 (2021), arXiv:2103.02159 [hep-ph] .
  * Barry  _et al._ (2021) P. C. Barry, C.-R. Ji, N. Sato, and W. Melnitchouk (Jefferson Lab Angular Momentum (JAM)), Phys. Rev. Lett. 127, 232001 (2021), arXiv:2108.05822 [hep-ph] .
  * Novikov  _et al._ (2020) I. Novikov _et al._ , Phys. Rev. D 102, 014040 (2020), arXiv:2002.02902 [hep-ph] .
  * Sutton  _et al._ (1992) P. J. Sutton, A. D. Martin, R. G. Roberts, and W. J. Stirling, Phys. Rev. D 45, 2349 (1992).
  * Gluck  _et al._ (1992) M. Gluck, E. Reya, and A. Vogt, Z. Phys. C 53, 651 (1992).
  * Whalley  _et al._ (2005) M. R. Whalley, D. Bourilkov, and R. C. Group, in  _HERA and the LHC: A Workshop on the Implications of HERA and LHC Physics (Startup Meeting, CERN, 26-27 March 2004; Midterm Meeting, CERN, 11-13 October 2004)_ (2005) pp. 575–581, arXiv:hep-ph/0508110 .
  * Buckley  _et al._ (2015) 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), arXiv:1412.7420 [hep-ph] .
  * Nason  _et al._ (1988) P. Nason, S. Dawson, and R. K. Ellis, Nucl. Phys. B 303, 607 (1988).
  * Nason  _et al._ (1989) P. Nason, S. Dawson, and R. K. Ellis, Nucl. Phys. B 327, 49 (1989), [Erratum: Nucl.Phys.B 335, 260–260 (1990)].
  * Mangano  _et al._ (1993) M. L. Mangano, P. Nason, and G. Ridolfi, Nucl. Phys. B 405, 507 (1993).
  * Eichten and Quigg (1995) E. J. Eichten and C. Quigg, Phys. Rev. D 52, 1726 (1995), arXiv:hep-ph/9503356 .
  * Vogt (2000) R. Vogt, Phys. Rev. C 61, 035203 (2000), arXiv:hep-ph/9907317 .
  * Maltoni  _et al._ (2006) F. Maltoni _et al._ , Phys. Lett. B 638, 202 (2006), arXiv:hep-ph/0601203 .
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Dulat  _et al._ (2016) S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin, C. Schmidt, D. Stump, and C. P. Yuan, Phys. Rev. D 93, 033006 (2016), arXiv:1506.07443 [hep-ph] .
  * Eskola  _et al._ (2017) K. J. Eskola, P. Paakkinen, H. Paukkunen, and C. A. Salgado, Eur. Phys. J. C 77, 163 (2017), arXiv:1612.05741 [hep-ph] .
  * Abt  _et al._ (2007) I. Abt _et al._ (HERA-B), Eur. Phys. J. C 49, 545 (2007), arXiv:hep-ex/0607046 .
  * Baglin  _et al._ (1995) C. Baglin _et al._ (NA38), Phys. Lett. B 345, 617 (1995).
  * Alexandrov  _et al._ (1999) Y. Alexandrov _et al._ (BEATRICE), Nucl. Phys. B 557, 3 (1999).
  * Corden  _et al._ (1980) M. J. Corden _et al._ (WA39), Phys. Lett. B 96, 411 (1980).
  * Gribushin  _et al._ (1996) A. Gribushin _et al._ (E672, E706), Phys. Rev. D 53, 4723 (1996).
  * Antoniazzi  _et al._ (1992) L. Antoniazzi _et al._ (E705), Phys. Rev. D 46, 4828 (1992).
  * Badier  _et al._ (1983) J. Badier _et al._ (NA3), Z. Phys. C 20, 101 (1983).
  * McEwen  _et al._ (1983) J. G. McEwen _et al._ (WA11), Phys. Lett. B 121, 198 (1983).
  * Akerlof  _et al._ (1993) C. Akerlof _et al._ (E537), Phys. Rev. D 48, 5067 (1993).
  * Corden  _et al._ (1981) M. J. Corden _et al._ (WA39), Phys. Lett. B 98, 220 (1981).
  * Heinrich  _et al._ (1991) J. G. Heinrich _et al._ (E615), Phys. Rev. D 44, 1909 (1991).
  * Alessandro  _et al._ (2004) B. Alessandro _et al._ (NA50), Eur. Phys. J. C 33, 31 (2004).
  * Schub  _et al._ (1995) M. H. Schub _et al._ (E789), Phys. Rev. D 52, 1307 (1995), [Erratum: Phys.Rev.D 53, 570 (1996)].
  * Alexopoulos  _et al._ (1996) T. Alexopoulos _et al._ (E771), Phys. Lett. B 374, 271 (1996).
  * Sup (2022) Supplemental materials (2022).
  * Kovarik  _et al._ (2016) K. Kovarik _et al._ , Phys. Rev. D 93, 085037 (2016), arXiv:1509.00792 [hep-ph] .
  * Butenschoen and Kniehl (2011) M. Butenschoen and B. A. Kniehl, Phys. Rev. Lett. 106, 022003 (2011), arXiv:1009.5662 [hep-ph] .
  * Chao  _et al._ (2012) K.-T. Chao, Y.-Q. Ma, H.-S. Shao, K. Wang, and Y.-J. Zhang, Phys. Rev. Lett. 108, 242004 (2012), arXiv:1201.2675 [hep-ph] .
  * Zhang  _et al._ (2015) H.-F. Zhang, Z. Sun, W.-L. Sang, and R. Li, Phys. Rev. Lett. 114, 092006 (2015), arXiv:1412.0508 [hep-ph] .
  * Bodwin  _et al._ (2014) G. T. Bodwin, H. S. Chung, U.-R. Kim, and J. Lee, Phys. Rev. Lett. 113, 022001 (2014), arXiv:1403.3612 [hep-ph] .
  * Ma  _et al._ (2011) Y.-Q. Ma, K. Wang, and K.-T. Chao, Phys. Rev. D 84, 114001 (2011), arXiv:1012.1030 [hep-ph] .
  * Gong  _et al._ (2013) B. Gong, L.-P. Wan, J.-X. Wang, and H.-F. Zhang, Phys. Rev. Lett. 110, 042002 (2013), arXiv:1205.6682 [hep-ph] .
  * Butenschoen and Kniehl (2022) M. Butenschoen and B. A. Kniehl, (2022), arXiv:2207.09346 [hep-ph] .
  * Arleo  _et al._ (2019) F. Arleo, C.-J. Naïm, and S. Platchkov, JHEP 01, 129 (2019), arXiv:1810.05120 [hep-ph] .
  * Bourrely  _et al._ (2022) C. Bourrely, W.-C. Chang, and J.-C. Peng, Phys. Rev. D 105, 076018 (2022), arXiv:2202.12547 [hep-ph] .
  * Aghasyan  _et al._ (2017) M. Aghasyan _et al._ (COMPASS), Phys. Rev. Lett. 119, 112002 (2017), arXiv:1704.00488 [hep-ex] .
  * Adams  _et al._ (2018) B. Adams _et al._ , (2018), arXiv:1808.00848 [hep-ex] .
  * Aguilar  _et al._ (2019) A. C. Aguilar _et al._ , Eur. Phys. J. A 55, 190 (2019), arXiv:1907.08218 [nucl-ex] .
  * Anderle  _et al._ (2021) D. P. Anderle _et al._ , Front. Phys. (Beijing) 16, 64701 (2021), arXiv:2102.09222 [nucl-ex] .
  * Chávez  _et al._ (2022) J. M. M. Chávez, V. Bertone, F. De Soto Borrero, M. Defurne, C. Mezrag, H. Moutarde, J. Rodríguez-Quintero, and J. Segovia, Phys. Rev. Lett. 128, 202501 (2022), arXiv:2110.09462 [hep-ph] .

## Supplemental materials

Figure 1: The $x_{F}$-dependent cross sections for $J/\psi$ and $\psi(2S)$ production and $R_{\psi}(x_{F})$ ratios in $\pi^{-}N$ and $pN$ interactions, following the order given in Table. 2. The symbol and value in parenthesis denote the particle type and momentum of beam. The solid red and dotted black curves represent the NRQCD results of GRV pion PDFs from the fit described in the text (“FIT”) and from the calculation using the LDMEs obtained in Ref. Hsieh  _et al._ (2021) (“REF”), respectively. The values of $\chi^{2}$ divided by the number of data point (ndp) for each data set are also shown. Figure 2: The $x_{F}$-dependent cross sections for $J/\psi$ and $\psi(2S)$ production and $R_{\psi}(x_{F})$ ratios in $\pi^{-}N$ and $pN$ interactions, following the order given in Table. II. The symbol and value in parenthesis denote the particle type and momentum of beam. The solid red and dotted black curves represent the NRQCD results of JAM pion PDFs from the fit described in the text (“FIT”) and from the calculation using the LDMEs obtained in Ref. Hsieh  _et al._ (2021) (“REF”), respectively. The values of $\chi^{2}$ divided by the number of data point (ndp) for each data set are also shown. Figure 3: The $x_{F}$-dependent cross sections for $J/\psi$ and $\psi(2S)$ production and $R_{\psi}(x_{F})$ ratios in $\pi^{-}N$ and $pN$ interactions, following the order given in Table. II. The symbol and value in parenthesis denote the particle type and momentum of beam. The solid red and dotted black curves represent the NRQCD results of xFitter pion PDFs from the fit described in the text (“FIT”) and from the calculation using the LDMEs obtained in Ref. Hsieh  _et al._ (2021) (“REF”), respectively. The values of $\chi^{2}$ divided by the number of data point (ndp) for each data set are also shown. Figure 4: Differential cross sections for $J/\psi$ production with a 300-GeV/$c$ $\pi^{-}$ beam from the E705 experiment Antoniazzi  _et al._ (1992). The data are compared to the NRQCD fit results for the SMRS, GRV, xFitter, and JAM PDFs. The total cross sections and $q\bar{q}$, $GG$, and $qG$ contributions are denoted as solid black, dashed blue, dotted red, and dot-dashed green lines, respectively. The uncertainty bands associated with JAM and xFitter PDFs are also shown. Figure 5: Differential cross sections for $J/\psi$ production with a 280-GeV/$c$ $\pi^{-}$ beam from the NA3 experiment Badier  _et al._ (1983). The data are compared to the NRQCD fit results for the SMRS, GRV, xFitter, and JAM PDFs. The total cross sections and $q\bar{q}$, $GG$, and $qG$ contributions are denoted as solid black, dashed blue, dotted red, and dot-dashed green lines, respectively. The uncertainty bands associated with JAM and xFitter PDFs are also shown. Figure 6: Differential cross sections for $J/\psi$ production with a 200-GeV/$c$ $\pi^{-}$ beam from the NA3 experiment Badier  _et al._ (1983). The data are compared to the NRQCD fit results for the SMRS, GRV, xFitter, and JAM PDFs. The total cross sections and $q\bar{q}$, $GG$, and $qG$ contributions are denoted as solid black, dashed blue, dotted red, and dot-dashed green lines, respectively. The uncertainty bands associated with JAM and xFitter PDFs are also shown. Figure 7: Differential cross sections for $J/\psi$ production with a 190-GeV/$c$ $\pi^{-}$ beam from the WA11 experiment McEwen  _et al._ (1983). The data are compared to the NRQCD fit results for the SMRS, GRV, xFitter, and JAM PDFs. The total cross sections and $q\bar{q}$, $GG$, and $qG$ contributions are denoted as solid black, dashed blue, dotted red, and dot-dashed green lines, respectively. The uncertainty bands associated with JAM and xFitter PDFs are also shown. Figure 8: Differential cross sections for $J/\psi$ production with a 150-GeV/$c$ $\pi^{-}$ beam from the NA3 experiment Badier  _et al._ (1983). The data are compared to the NRQCD fit results for the SMRS, GRV, xFitter, and JAM PDFs. The total cross sections and $q\bar{q}$, $GG$, and $qG$ contributions are denoted as solid black, dashed blue, dotted red, and dot-dashed green lines, respectively. The uncertainty bands associated with JAM and xFitter PDFs are also shown. Figure 9: Differential cross sections for $J/\psi$ production with a 125-GeV/$c$ $\pi^{-}$ beam from the E537 experiment Akerlof  _et al._ (1993). The data are compared to the NRQCD fit results for the SMRS, GRV, xFitter, and JAM PDFs. The total cross sections and $q\bar{q}$, $GG$, and $qG$ contributions are denoted as solid black, dashed blue, dotted red, and dot-dashed green lines, respectively. The uncertainty bands associated with JAM and xFitter PDFs are also shown. Figure 10: The $\psi(2S)$ to $J/\psi$ cross section ratios $R_{\psi}(x_{F})$ for $J/\psi$ and $\psi(2S)$ production with a 920-GeV/$c$ proton beam from the HERA-B experiment Abt  _et al._ (2007). The data are compared to the NRQCD fit results for the SMRS, GRV, xFitter, and JAM PDFs. The ratios of total cross sections and individual $R^{q\bar{q}}_{\psi}(x_{F})$ and $R^{GG}_{\psi}(x_{F})$ contributions are denoted as solid black, dashed blue, and dotted red lines, respectively. Figure 11: The $\psi(2S)$ to $J/\psi$ cross section ratios $R_{\psi}(x_{F})$ for $J/\psi$ and $\psi(2S)$ production with a 450-GeV/$c$ proton beam from the NA50 experiment Alessandro  _et al._ (2004). The data are compared to the NRQCD fit results for the SMRS, GRV, xFitter, and JAM PDFs. The ratios of total cross sections and individual $R^{q\bar{q}}_{\psi}(x_{F})$ and $R^{GG}_{\psi}(x_{F})$ contributions are denoted as solid black, dashed blue, and dotted red lines, respectively. Figure 12: The $\psi(2S)$ to $J/\psi$ cross section ratios $R_{\psi}(x_{F})$ for $J/\psi$ and $\psi(2S)$ production with a 450-GeV/$c$ proton beam from the NA50 experiment Alessandro  _et al._ (2004). The data are compared to the NRQCD fit results for the SMRS, GRV, xFitter, and JAM PDFs. The ratios of total cross sections and individual $R^{q\bar{q}}_{\psi}(x_{F})$ and $R^{GG}_{\psi}(x_{F})$ contributions are denoted as solid black, dashed blue, and dotted red lines, respectively. Figure 13: The $\psi(2S)$ to $J/\psi$ cross section ratios $R_{\psi}(x_{F})$ for $J/\psi$ and $\psi(2S)$ production with a 800-GeV/$c$ proton beam from the E771 experiment Alexopoulos  _et al._ (1996). The ratios of total cross sections and individual $R^{q\bar{q}}_{\psi}(x_{F})$ and $R^{GG}_{\psi}(x_{F})$ contributions are denoted as solid black, dashed blue, and dotted red lines, respectively. Figure 14: Differential cross sections for $J/\psi$ production with a 300-GeV/$c$ proton beam from the E705 experiment Antoniazzi  _et al._ (1992). The data are compared to the NRQCD fit results of LDMEs for the SMRS, GRV, xFitter, and JAM PDFs. The total cross sections and $q\bar{q}$, $GG$, and $qG$ contributions are denoted as solid black, dashed blue, dotted red, and dot-dashed green lines, respectively. Figure 15: Differential cross sections for $J/\psi$ production with a 200-GeV/$c$ proton beam from the NA3 experiment Badier  _et al._ (1983). The data are compared to the NRQCD fit results of LDMEs for the SMRS, GRV, xFitter, and JAM PDFs. The total cross sections and $q\bar{q}$, $GG$, and $qG$ contributions are denoted as solid black, dashed blue, dotted red, and dot-dashed green lines, respectively. Figure 16: Integrated charmonium cross sections in $pN$ and $\pi^{-}N$ collisions. The data for $J/\psi$ and $\psi(2S)$ production are compared to the fit made using the GRV pion PDFs. The total cross section and its $q\bar{q}$ and $GG$ contributions are denoted as solid black, dashed blue and dotted red lines, respectively. Figure 17: Integrated charmonium cross sections in $pN$ and $\pi^{-}N$ collisions. The data for $J/\psi$ and $\psi(2S)$ production are compared to the fit made using the JAM pion PDFs. The total cross section and its $q\bar{q}$ and $GG$ contributions are denoted as solid black, dashed blue and dotted red lines, respectively. The uncertainty bands associated with JAM PDFs are also shown. Figure 18: Integrated charmonium cross sections in $pN$ and $\pi^{-}N$ collisions. The data for $J/\psi$ and $\psi(2S)$ production are compared to the fit made using the xFitter pion PDFs. The total cross section and its $q\bar{q}$ and $GG$ contributions are denoted as solid black, dashed blue and dotted red lines, respectively. The uncertainty bands associated with xFitter PDFs are also shown. Figure 19: The NRQCD results with variation of charm quark mass $m_{c}$ and renormalization scale $\mu_{R}$, compared with the $d\sigma/dx_{F}$ data of $J/\psi$ production off the beryllium target with a 515-GeV/$c$ $\pi^{-}$ beam from the E672/E706 experiment Gribushin  _et al._ (1996). The pion PDFs used for the calculation is GRV. The total cross sections and $q\bar{q}$, $GG$, and $qG$ contributions are denoted as solid black, dashed blue, dotted red and dot-dashed green lines, respectively. The charm quark mass $m_{c}$, factorization scale $\mu_{F}$, and renormalization scale $\mu_{R}$ used for the NRQCD calculation as well as the fit $\chi^{2}$/ndf are displayed in each plot. Figure 20: Same as Fig. 19 but with the input of JAM pion PDFs. Figure 21: Same as Fig. 19 but with the input of xFitter pion PDFs. | SMRS
---|---
|  $m_{c}=1.5$ GeV/$c^{2}$ |  $m_{c}=1.2$ GeV/$c^{2}$
$\mu/m_{c}$ | 1 | 2 | 3 | 1 | 2 | 3
$\chi^{2}_{total}/\text{ndf}$ | 2.5 | 1.9 | 2.3 | 393.9 | 44.1 | 6.1
$\chi^{2}/\text{ndp}|^{\pi^{-}}_{x_{F}}$ | 2.3 | 1.8 | 2.2 | 413.0 | 48.6 | 6.4
$\chi^{2}/\text{ndp}|^{p}_{x_{F}}$ | 2.3 | 1.6 | 2.0 | 202.6 | 15.0 | 2.9
$\chi^{2}/\text{ndp}|^{\pi^{-}}_{\sqrt{s}}$ | 4.6 | 8.7 | 7.3 | 240.6 | 16.6 | 5.2
$\chi^{2}/\text{ndp}|^{p}_{\sqrt{s}}$ | 5.0 | 8.1 | 9.2 | 1305.1 | 93.4 | 25.4
$\langle\mathcal{O}_{8}^{J/\psi}[^{3}S_{1}]\rangle$ | 1.6E-02 | 2.6E-02 | 5.5E-02 | 2.5E-10 | 5.5E-13 | 2.1E-09
|  $\pm$1.6E-03 |  $\pm$2.3E-03 |  $\pm$4.1E-03 |  $\pm$1.5E-07 |  $\pm$2.2E-06 |  $\pm$1.6E-05
$\Delta_{8}^{J/\psi}$ | 1.3E-02 | 5.6E-02 | 8.1E-02 | 1.1E-14 | 5.0E-13 | 1.9E-11
|  $\pm$7.9E-04 |  $\pm$1.6E-03 |  $\pm$2.0E-03 |  $\pm$5.3E-08 |  $\pm$7.1E-07 |  $\pm$4.7E-06
$\langle\mathcal{O}_{8}^{\psi(2S)}[^{3}S_{1}]\rangle$ | 7.7E-03 | 1.3E-02 | 2.0E-02 | 1.9E-10 | 1.5E-11 | 2.8E-03
|  $\pm$4.0E-04 |  $\pm$8.6E-04 |  $\pm$1.3E-03 |  $\pm$2.0E-07 |  $\pm$6.4E-06 |  $\pm$2.0E-04
$\Delta_{8}^{\psi(2S)}$ | 2.5E-04 | 5.7E-03 | 8.3E-03 | 5.7E-12 | 8.2E-12 | 6.3E-10
|  $\pm$1.9E-04 |  $\pm$2.9E-04 |  $\pm$4.9E-04 |  $\pm$2.1E-07 |  $\pm$1.4E-06 |  $\pm$9.6E-06
Table 5: The reduced $\chi^{2}/\text{ndf}$ of values for the whole data sets and the $\chi^{2}$ divided by the number of data point (ndf) for the pion-induced and proton-induced data sets with the systematic variation of charm quark mass $m_{c}$ between 1.5 and 1.2 GeV/$c^{2}$, and $\mu=\mu_{R}=\mu_{F}$ at 0.5, 1.0, and 2.0 $m_{c}$ in NRQCD calculations and the corresponding input or best-fit LDMEs for SMRS pion PDFs. All LDMEs are in units of $\rm{GeV}^{3}$. | GRV
---|---
|  $m_{c}=1.5$ GeV/$c^{2}$ |  $m_{c}=1.2$ GeV/$c^{2}$
$\mu/m_{c}$ | 1 | 2 | 3 | 1 | 2 | 3
$\chi^{2}_{total}/\text{ndf}$ | 2.4 | 2.4 | 2.6 | 616.4 | 54.2 | 7.7
$\chi^{2}/\text{ndp}|^{\pi^{-}}_{x_{F}}$ | 2.3 | 2.4 | 2.6 | 681.0 | 61.1 | 8.4
$\chi^{2}/\text{ndp}|^{p}_{x_{F}}$ | 2.0 | 1.7 | 1.8 | 202.6 | 14.1 | 2.7
$\chi^{2}/\text{ndp}|^{\pi^{-}}_{\sqrt{s}}$ | 8.4 | 5.6 | 3.6 | 479.5 | 16.8 | 6.7
$\chi^{2}/\text{ndp}|^{p}_{\sqrt{s}}$ | 5.5 | 8.1 | 9.6 | 1305.1 | 93.4 | 25.8
$\langle\mathcal{O}_{8}^{J/\psi}[^{3}S_{1}]\rangle$ | 1.5E-04 | 4.3E-02 | 9.4E-02 | 2.7E-10 | 3.2E-13 | 7.7E-09
|  $\pm$1.3E-04 |  $\pm$3.8E-03 |  $\pm$5.5E-03 |  $\pm$1.0E-07 |  $\pm$2.5E-06 |  $\pm$2.5E-05
$\Delta_{8}^{J/\psi}$ | 1.9E-02 | 5.2E-02 | 7.4E-02 | 3.6E-12 | 3.2E-14 | 6.8E-09
|  $\pm$1.4E-04 |  $\pm$1.7E-03 |  $\pm$2.0E-03 |  $\pm$9.0E-08 |  $\pm$5.9E-07 |  $\pm$4.8E-06
$\langle\mathcal{O}_{8}^{\psi(2S)}[^{3}S_{1}]\rangle$ | 8.4E-03 | 2.1E-02 | 3.2E-02 | 1.1E-10 | 5.5E-11 | 3.1E-03
|  $\pm$5.8E-04 |  $\pm$1.3E-03 |  $\pm$1.4E-03 |  $\pm$5.5E-07 |  $\pm$5.6E-06 |  $\pm$2.4E-04
$\Delta_{8}^{\psi(2S)}$ | 5.2E-04 | 4.2E-03 | 6.0E-03 | 1.4E-10 | 8.4E-14 | 2.1E-08
|  $\pm$2.5E-04 |  $\pm$2.9E-04 |  $\pm$3.2E-04 |  $\pm$3.8E-08 |  $\pm$1.0E-06 |  $\pm$6.7E-06
Table 6: The reduced $\chi^{2}/\text{ndf}$ of values for the whole data sets and the $\chi^{2}$ divided by the number of data point (ndf) for the pion-induced and proton-induced data sets with the systematic variation of charm quark mass $m_{c}$ between 1.5 and 1.2 GeV/$c^{2}$, and $\mu=\mu_{R}=\mu_{F}$ at 0.5, 1.0, and 2.0 $m_{c}$ in NRQCD calculations and the corresponding input or best-fit LDMEs for GRV pion PDFs. All LDMEs are in units of $\rm{GeV}^{3}$. | JAM
---|---
|  $m_{c}=1.5$ GeV/$c^{2}$ |  $m_{c}=1.2$ GeV/$c^{2}$
$\mu/m_{c}$ | 1 | 2 | 3 | 1 | 2 | 3
$\chi^{2}_{total}/\text{ndf}$ | 8.2 | 5.6 | 5.0 | 229.6 | 30.1 | 16.6
$\chi^{2}/\text{ndp}|^{\pi^{-}}_{x_{F}}$ | 8.7 | 5.9 | 5.2 | 215.1 | 30.7 | 17.6
$\chi^{2}/\text{ndp}|^{p}_{x_{F}}$ | 4.0 | 2.7 | 2.6 | 202.6 | 18.4 | 7.7
$\chi^{2}/\text{ndp}|^{\pi^{-}}_{\sqrt{s}}$ | 31.3 | 11.4 | 6.7 | 39.9 | 3.3 | 15.6
$\chi^{2}/\text{ndp}|^{p}_{\sqrt{s}}$ | 6.8 | 5.1 | 6.4 | 1305.1 | 114.0 | 62.3
$\langle\mathcal{O}_{8}^{J/\psi}[^{3}S_{1}]\rangle$ | 6.0E-02 | 1.2E-01 | 1.7E-01 | 2.0E-11 | 3.0E-03 | 2.2E-02
|  $\pm$1.2E-03 |  $\pm$2.1E-03 |  $\pm$3.7E-03 |  $\pm$8.0E-07 |  $\pm$3.6E-04 |  $\pm$6.0E-04
$\Delta_{8}^{J/\psi}$ | 3.5E-03 | 2.4E-02 | 4.4E-02 | 9.9E-12 | 4.5E-08 | 5.2E-08
|  $\pm$6.9E-04 |  $\pm$1.6E-03 |  $\pm$2.6E-03 |  $\pm$2.6E-07 |  $\pm$3.0E-06 |  $\pm$1.6E-05
$\langle\mathcal{O}_{8}^{\psi(2S)}[^{3}S_{1}]\rangle$ | 1.1E-02 | 2.4E-02 | 3.3E-02 | 9.9E-10 | 2.8E-03 | 6.3E-03
|  $\pm$3.6E-04 |  $\pm$8.5E-04 |  $\pm$1.3E-03 |  $\pm$2.8E-06 |  $\pm$1.3E-04 |  $\pm$2.2E-04
$\Delta_{8}^{\psi(2S)}$ | 1.9E-09 | 2.1E-03 | 4.2E-03 | 2.5E-13 | 3.3E-08 | 1.1E-09
|  $\pm$1.1E-05 |  $\pm$3.2E-04 |  $\pm$5.0E-04 |  $\pm$3.7E-07 |  $\pm$2.9E-06 |  $\pm$7.8E-06
Table 7: The reduced $\chi^{2}/\text{ndf}$ of values for the whole data sets and the $\chi^{2}$ divided by the number of data point (ndf) for the pion-induced and proton-induced data sets with the systematic variation of charm quark mass $m_{c}$ between 1.5 and 1.2 GeV/$c^{2}$, and $\mu=\mu_{R}=\mu_{F}$ at 0.5, 1.0, and 2.0 $m_{c}$ in NRQCD calculations and the corresponding input or best-fit LDMEs for JAM pion PDFs. All LDMEs are in units of $\rm{GeV}^{3}$. | xFitter
---|---
|  $m_{c}=1.5$ GeV/$c^{2}$ |  $m_{c}=1.2$ GeV/$c^{2}$
$\mu/m_{c}$ | 1 | 2 | 3 | 1 | 2 | 3
$\chi^{2}_{total}/\text{ndf}$ | 4.8 | 4.2 | 4.2 | 99.1 | 15.1 | 10.2
$\chi^{2}/\text{ndp}|^{\pi^{-}}_{x_{F}}$ | 4.6 | 4.5 | 4.5 | 58.0 | 13.1 | 10.8
$\chi^{2}/\text{ndp}|^{p}_{x_{F}}$ | 3.9 | 1.9 | 1.9 | 202.6 | 16.9 | 5.0
$\chi^{2}/\text{ndp}|^{\pi^{-}}_{\sqrt{s}}$ | 9.7 | 4.4 | 3.1 | 49.5 | 3.4 | 3.8
$\chi^{2}/\text{ndp}|^{p}_{\sqrt{s}}$ | 11.0 | 6.9 | 8.4 | 1305.1 | 101.3 | 42.2
$\langle\mathcal{O}_{8}^{J/\psi}[^{3}S_{1}]\rangle$ | 3.6E-02 | 8.5E-02 | 1.3E-01 | 5.2E-11 | 3.4E-08 | 1.2E-02
|  $\pm$9.9E-04 |  $\pm$4.1E-03 |  $\pm$1.0E+00 |  $\pm$9.1E-07 |  $\pm$6.8E-05 |  $\pm$5.2E-04
$\Delta_{8}^{J/\psi}$ | 1.7E-02 | 3.9E-02 | 6.3E-02 | 1.3E-10 | 8.7E-09 | 1.3E-08
|  $\pm$9.5E-04 |  $\pm$3.4E-03 |  $\pm$9.4E+00 |  $\pm$5.5E-07 |  $\pm$6.1E-06 |  $\pm$1.2E-05
$\langle\mathcal{O}_{8}^{\psi(2S)}[^{3}S_{1}]\rangle$ | 9.3E-03 | 1.9E-02 | 2.7E-02 | 2.1E-09 | 2.3E-03 | 4.6E-03
|  $\pm$3.6E-04 |  $\pm$1.2E-03 |  $\pm$1.1E+00 |  $\pm$9.4E-07 |  $\pm$1.2E-04 |  $\pm$1.8E-04
$\Delta_{8}^{\psi(2S)}$ | 9.1E-04 | 4.0E-03 | 6.6E-03 | 6.8E-11 | 3.2E-13 | 2.8E-09
|  $\pm$1.8E-04 |  $\pm$6.4E-04 |  $\pm$9.7E+00 |  $\pm$7.0E-07 |  $\pm$4.1E-06 |  $\pm$1.1E-05
Table 8: The reduced $\chi^{2}/\text{ndf}$ of values for the whole data sets
and the $\chi^{2}$ divided by the number of data point (ndf) for the pion-
induced and proton-induced data sets with the systematic variation of charm
quark mass $m_{c}$ between 1.5 and 1.2 GeV/$c^{2}$, and $\mu=\mu_{R}=\mu_{F}$
at 0.5, 1.0, and 2.0 $m_{c}$ in NRQCD calculations and the corresponding input
or best-fit LDMEs for xFitter pion PDFs. All LDMEs are in units of
$\rm{GeV}^{3}$.