Constraining Gluon Density of Pions at Large $x$ by Pion-induced $J/\psi$ Production

The gluon distributions of the pion obtained from various global fits exhibit large variations among them. We show that the existing pion-induced $J/\psi$ production data, usually not included in the global fits, can impose useful additional constraints on the pion PDFs. In particular, these data can probe pion's gluon densities at large $x$. Existing pion-induced $J/\psi$ data covering a broad range of beam momenta are compared with calculations using various sets of pion PDFs. The $J/\psi$ production cross sections are evaluated in the framework of next-to-leading order color evaporation model. It is found that $J/\psi$ data measured at forward rapidity and at sufficiently high beam momentum are sensitive to the large-$x$ gluon distribution of pions. The current $J/\psi$ data strongly favor the SMRS and GRV pion PDFs, containing significant gluon contents at large $x$.


I. INTRODUCTION
The pion, as the Goldstone boson of dynamical chiral symmetry breaking of the strong interaction, is the lightest QCD bound state. Because of its light mass, pion plays a dominant role in the long-range nucleon-nucleon interaction [1]. Understanding pion's internal structure is important to investigate the low-energy, nonperturbative aspects of QCD [2]. Even though the pion is theoretically simpler than the proton, its partonic structure is much less explored. As scattering off a pion target is not feasible, current knowledge on pion parton distribution functions (PDFs) mostly relies on the pion-induced Drell-Yan data. Since these fixed-target data are mostly sensitive to the valence quark distributions at x > 0.2, the sea and gluon densities are essentially unconstrained.
In principle, the prompt photon production process πN → γX can constrain the gluon content of pions through the Gq → γq subprocess, but the experimental uncertainties are large. Production of heavy quarkonia, like J/ψ and Υ(1S), with the pion beam has distinctive advantages: the cross sections are large and they can be readily detected via the dimuon decay channel. These data sets have been shown to be sensitive to both the quark and gluon distributions of the incident pion with model-dependent assumptions of quarkonia fragmentation [3,4]. The interesting possibility of accessing the pion PDFs from leading neutron deep inelastic scattering (DIS) data, has been considered with promising results [5,6]. However, this method is subject to large systematic uncertainties and further studies on the uncertainties of pion splitting function and the off-shellness of virtual pion are required [7,8]. To precisely determine the sea quark content of pions, there was a suggestion of performing the Drell-Yan measurement with π + and π − beams on the isoscalar deuterium target [9] and such measurement is planned in a future experiment [10].
Until a couple of years ago, knowledge of the pion PDFs was limited to global analyses carried out more than two decades ago: OW [11], ABFKW [12], SMRS [13], GRV [14] and GRS [15]. These analyses were based mostly on pioninduced Drell-Yan, often on prompt-photon and in some cases on J/ψ production data. New analyses were performed only recently, using the same Drell-Yan data in BS [16] as well as both the Drell-Yan and direct-photon data in xFitter [17]. The analysis of JAM [18] included both the Drell-Yan data and, for the first time, the leading neutron tagged electroproduction data. The experimental situation is also evolving. After more than two decades, a new measurement of pion-induced Drell-Yan production cross sections was performed by the CERN COMPASS Collaboration [19]. The data are expected to be available in the near future. A proposal dedicated to investigating the pion and kaon structure at the future electron-ion collider in the U.S. [20] was recently described in Ref. [21].
In this work we investigate the sensitivity of the J/ψ production data to the pion PDFs. The theoretical challenge of this reaction comes from the treatment of the hadronization of cc pairs into a charmonium bound state. This non-perturbative process has been modeled in several theoretical approaches including color evaporation model (CEM) [42], color-singlet model (CSM) [43] and non-relativistic QCD (NRQCD) [44]. The CEM assumes a constant probability for cc pairs to hadronize into a given charmonium. In CSM, the production of J/ψ is assumed to be through the color-singlet cc channel of the same quantum numbers as J/ψ. The NRQCD expands the calculations by the powers of the average velocity of cc pairs in the rest frame of J/ψ. The hadronization probability of each cc pair depends on its color and spin state. More details about these theoretical frameworks can be found in Ref. [45]. In general, CSM and NRQCD provide a good description of data taken at collider energies, but fail to explain measurements at fixed-target energies [46]. In contrast, the more phenomenological CEM gives good account of many features of fixed-target J/ψ cross section data, including their longitudinal momentum (x F ) distributions [47,48]. Therefore, we adopt the CEM approach in this work, even though this model has some limitations [49].
In the fixed-target energy domain, where the transverse momentum of the charmonium is less than its mass, the charmonium production is dominated by the quark-antiquark (qq) and gluon-gluon fusion (GG) partonic processes. The shape of the longitudinal momentum x F cross section is therefore sensitive to the quark and gluon parton distributions of colliding hadrons. Since the nucleon PDFs are known with good accuracy, the measurement of the x F distribution of J/ψ production with the pion beam provides, within the theoretical model uncertainties, valuable information about the pion quark and gluon partonic distributions. Our study is performed using next-to-leading order (NLO) CEM calculation, including the recent nucleon PDFs. The available pion-induced J/ψ data on hydrogen and several light-mass nuclear targets are compared to calculations using the available pion PDFs. Over the broad energy range considered, all pion PDF sets provide reasonable agreement with the x F -integrated cross sections. In contrast, for the x F distributions, we find that the agreement between data and calculations strongly depends on the magnitude and shape of the pion gluon distribution. This paper is organized as follows. In Sec. II, we briefly describe the CEM framework for the calculations of J/ψ production cross sections in the collisions of pions and nucleons. Some distinctive features of parton densities in various pion PDFs used for the calculations are presented in Sec. III. The NLO CEM calculations using various pion PDFs are compared with the existing J/ψ production data in Sec. IV. Section V shows the results of systematic study for the CEM calculation. We discuss our findings from the comparison of CEM calculations with data in Sec. VI, followed by a summary in Sec. VII.

II. COLOR EVAPORATION MODEL AND HEAVY QUARK PAIR PRODUCTION
The theoretical treatment of heavy quarkonium production consists of the QCD description of the production of heavy quark pairs (QQ) at the parton level, and their subsequent hadronization into the quarkonium states. One of the theoretical approaches is the non-relativistic QCD (NRQCD) [44] where the cross section of quarkonium production is expanded in terms of the strong coupling constant α S and the QQ velocity. The cross section is factorized into the hard and soft parts for each color and spin state of the QQ pairs. The shortdistance hard part is calculated perturbatively as a series of α S in pQCD. The soft part consists of long-distance matrix elements (LDMEs) characterizing the probability of hadronization process for each color and spin state. The LDMEs are determined by a fit to the experimental data. In the color sin-glet model [43], the production channel is assumed to be the color-singlet QQ state with quantum numbers exactly matching those of the heavy quarkonium.
Based on quark-hadron duality, the color evaporation model (CEM) assumes a constant probability for QQ pairs to hadronize into a quarkonium state. Taking J/ψ as an example, one first produces a cc pair via various QCD hard processes. For cc with an invariant mass M cc less than the DD threshold, a constant probability F , specific for each quarkonium, accounts for the hadronization of cc pairs into the colorless J/ψ state.
In CEM, the differential cross section, dσ/dx F , for J/ψ from the πN collision are expressed as where i, j denote the interacting partons (gluons, quarks and antiquarks), and m c , m D and M cc are the masses of the charm quark, D meson, and cc pair, respectively. The f π and f N are the corresponding pion and nucleon parton distribution functions evaluated at the corresponding Bjorken-x, x 1 and x 2 , at the factorization scale µ F . The short-distance differential cross section of heavy-quark pair productionσ[ij → ccX] is calculable as a perturbation series in the strong coupling α s (µ R ) evaluated at the renormalization scale µ R . The variable s is the square of the centerof-mass energy of the colliding π-N system, and p L is the longitudinal momentum of detected dimuon pair in the centerof-mass frame of π-N . It is assumed that the momenta of J/ψ and cc are approximately the same.
As mentioned above, the hadronization factor F is assumed to be universal, independent of the kinematics and the spin state of cc and the production subprocess. Therefore a unique feature of the CEM calculation is that the relative weight of each subprocess in dσ/dx F , is solely fixed by the convolution of partonic level cross sectionsσ and associated parton density distributions f π and f N , and, in particular, is independent of the F factor. The F factor is to be determined as the normalization parameter in the fit to the experimental measurements. The assumption of a common F factor for different subprocesses greatly reduces the number of free parameters in CEM.
The leading-order (O(α 2 S )) calculations of hard QCD kernelσ[ij → ccX] include the quark-antiquark (qq) and gluongluon fusion (GG) diagrams. Additional quark-gluon Compton scattering (Gq, Gq) and virtual gluon corrections enter the NLO (O(α 3 S )) calculations. The contributing partonic subprocesses in the fixed-order LO and NLO calculations are listed explicitly below [50]: Inclusion of both real and virtual gluon emission diagrams is necessary for calculating the full O(α 3 S ) cross sections. In this work, we utilize the theoretical framework of NLO calculation of the total cross sections for production of heavy quark pair, developed by Nason et al. [50][51][52]. This framework has been widely used in the calculation of heavy quark production. For example, it has been adopted in the NLO calculation of CEM for J/ψ production in hadronic collisions [47,48,53]. With a few parameters including the heavy quark mass m c and hadronization factor F , the CEM calculations adequately reproduced the fixed-target data with proton, antiproton and pion beams [47,48], as well as the collider data [53].

III. PION PDFS
Pion-induced Drell-Yan data are included in all global analyses for the determination of the pion PDFs. However, Drell-Yan data [54][55][56][57] constrain mainly the valence quark distribution. Without additional observables, the sea and gluon distributions remain practically unknown. Their magnitude can only be inferred through the momentum sum rule and valence quark sum rule. Different approaches have been taken to access the gluon and sea quark distributions: i) utilizing J/ψ production data in OW [11]; ii) utilizing the direct photon production data [58] in ABFKW [12], SMRS [13], GRV [14] and xFitter [17]; iii) utilizing leading neutron DIS [59,60] in JAM [18]. In addition, some pion PDFs are based on theoretical modeling. For example, the GRS [15] utilized a constituent quark model to related the gluon and antiquark density, and the BS [16] assumed quantum statistical distributions for all parton species with an universal temperature. We note that the OW analysis was performed at LO, whereas a NLO fit was carried out for all other analyses. Uncertainty bands for the resulting parton density distributions are available for two most recent global fits, JAM and xFitter. It was recently shown that the soft-gluon threshold resummation correction modifies the extraction of valence quark distribution, and particularly its falloff towards x = 1 [61]. This correction has not been implemented in any of the pion global analyses yet and it should only affect the calculated shape at the highest x F region.   lower by 20-30%. The sizable error bands of the sea distributions provided by JAM and xFitter clearly indicate that the pion sea remains poorly known. As for the gluon distributions, the early PDF sets of OW, ABFKW, SMRS and GRV have relatively large densities for x > 0.1, at variance with the recent xFitter and JAM PDFs that lie significantly lower. The spread of the gluon distributions around x = 0.5 among these six PDFs is even larger than the uncertainties of xFitter and JAM PDFs. Table I lists the momentum fractions of valence quarks (ū val ), sea quarks (ū sea ) and gluons (G) of negative pions estimated by various pion PDFs at Q 2 = 9.6 GeV 2 , follow- The values for the valence quarks show differences of up to 15%-20%, but are nearly equal for the two more recent PDFs, JAM and xFitter. The gluon first moments vary from 0.29 for xFitter to 0.51 for GRV. The low gluon value in xFitter is compensated by a much larger sea contribution.

IV. RESULTS OF NLO CEM CALCULATIONS
In this section we explore the sensitivity of the NLO CEM calculations to the various global fit parametrizations of the pion PDFs. We select four of them, namely SMRS and GRV, as the most widely used for a long time, and the two most recent fits, xFitter and JAM. 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 GeV 2 . As illustrated in Fig. 1, SMRS, JAM and xFitter have quite similar valence quark distributions while the magnitude of the GRV distribution is lower, by up to 20-30%. As for the gluon distributions, SMRS and GRV have similar shapes and magnitudes, while the magnitudes of xFitter and JAM are significantly smaller, by a factor of 2-4.
As a first step, we compare the NLO CEM cross sections integrated over x F > 0 for the process π − N → J/ψX for each of the four pion PDFs with the available measurements as a function of the center of mass energy, √ s, of the reaction. The calculations are performed using the nucleon CT14nlo PDFs [62] under the LHAPDF framework [63,64]. The cross sections are evaluated with a charm quark mass m c = 1.5 GeV/c 2 and renormalization and factorization scales of µ R = m c and µ F = 2m c , respectively [52]. The experimental cross sections are taken from the compilation of Ref. [45]. For the sake of completeness, the subsequent measurement from the WA92 experiment [65] is also included, after correcting it for the nuclear dependence. The hadronization factors F are assumed to be energy independent and are determined by the best fit to the data for each pion PDF.
The results and the comparison with data are displayed in Fig. 2. The total cross sections for the four PDFs exhibit quite similar √ s dependencies and all agree reasonably with the data. The differences between them are visible through the F factors, which vary from 0.05 to 0.09. As a general feature, the qq contribution dominates at low energies, whereas the GG contribution becomes increasingly important with increasing √ s. However, the relative fractions of qq and GG contributions as a function of √ s vary considerably, reflecting the differences between the corresponding parton distributions. For and J/ψ production cross sections at xF > 0 for the π − N reaction, calculated with four pion PDFs: SMRS, GRV, xFitter and JAM, is compared with data (solid circles) [45,65]. The black, blue and red curves represent the calculated total cross section, the qq and GG contributions, respectively. The shaded bands on the xFitter and JAM calculations come from the uncertainties of the corresponding PDF sets. The SMRS and GRV PDFs contain no information on uncertainties.
SMRS and GRV, the GG contribution starts to dominate the cross section beyond √ s = 13 and √ s = 11 GeV, respectively, while for xFitter and JAM the corresponding values are much higher: √ s = 19 and 21 GeV, respectively. In order to investigate further the effect led by different pion PDFs, we compare the longitudinal x F distribution of the calculated pion-induced J/ψ production cross section with a selection of fixed-target data from Fermilab and CERN experiments. Among the data sets available for pion-induced J/ψ production [66][67][68][69][70][71][72][73][74][75], we choose the ones that have large-x F coverage for either hydrogen or light nuclear targets (lithium and beryllium) in order to minimize the effects of the nuclear environment. The selected eight data sets are listed in Table II. The beam momenta of the data sets cover the range of 39.5 to 515 GeV/c, corresponding to √ s values ranging from 8.6 to 31.1 GeV. Some of the data listed in Table II involve nuclear targets. The target PDFs parametrizations are CT14nlo for the hydrogen target and EPPS16 [76] for the lithium and beryllium targets. Contrary to the integrated cross sections, we now allow the hadronization factor F to be fine-tuned for each data set individually.
Within the CEM and heavy quark pair production framework introduced in Sec. II, we performed the LO and NLO calculations of the differential cross sections as a function of x F with the charm quark mass m c = 1.5 GeV/c 2 , renormalization and factorization scales µ R = m c and µ F = 2m c [52]. TABLE II. The J/ψ production data sets with π − beam used in the analysis, listed in the order of decreasing beam momentum. The hadronization factor F , as an overall normalization parameter, is determined by the best fit to the x F distributions of cross sections, shown as the black lines in Figs. 3-10. The χ 2 /ndf value of the best fit is also displayed in the plot. The estimated individual qq and GG contributions are denoted as the blue and red lines, respectively. There is a negligible additional contribution from the qG subprocess, shown as green lines, to the total cross sections in the NLO calculation. The calculated value of the qG contribution is negative [50]. The variations coming from the uncertainty of xFitter and JAM PDF sets are displayed as shaded bands. In the following subsections (Secs. IV A-IV H), we briefly comment on the features of each experimental measurement and discuss the comparison of the data with the CEM calculations. Our observations are summarized in Sec. IV I.

A. Fermilab E672/E706 experiment
The Fermilab E672/E706 experiment [66] used a 515 GeV/c π − beam scattered off 3.71 and 1.12 cm long 9 Be targets. About 9600 J/ψ events integrated in the mass region between 2.8 and 3.4 GeV/c 2 were collected. The final cross sections cover the range 0.1 ≤ x F ≤ 0.8 in bins of 0.02 and have a normalization uncertainty of 12%.
The comparison of our calculations at both LO and NLO to the E672/E706 data is shown in Fig. 3. Judging from the reduced χ 2 /ndf values, the NLO calculations with SMRS, GRV and xFitter are in reasonable agreement with the data. The NLO calculation generally improves the agreement with the E672/E706 data, except for JAM. In comparison with the LO, the NLO calculation has a large effect on the cross sections, increasing its magnitude by more than a factor of two. An interesting observation is that this increase in magnitude is nearly entirely compensated by the F factor, pointing to a nearly uniform increase along x F . We also note that the the GG contribution dominates the cross section up to values of x F as large as 0.5-0.7, depending on the particular pion PDF set. The additional qG term in the NLO calculation has a minor (and negative) contribution, although largely dependent on the particular PDF set.

B. Fermilab E705 experiment
The Fermilab E705 experiment [67] used a 300 GeV/c negative hadron beam (with 98% pions) scattered off a 33 cm long lithium target. Data were also collected with a positive hadron beam consisting of protons and positive pions. Thanks to open geometry spectrometer, an excellent mass resolution was achieved, allowing a measurement of the J/ψ peak in the mass range between 2.98 and 3.18 GeV/c 2 . Since the final number of J/ψ events was not explicitly given, we estimate it from the published statistical errors to about 6000 events for the negative pion beam. The final cross sections have a normalization uncertainty of 11.1% and cover the range −0.1 ≤ x F ≤ 0.45 in bins of 0.05.
The comparison of our calculations with the experimental cross sections is shown in Fig. 4. The best χ 2 /ndf value is obtained with the SMRS PDFs. In contrast, the use of the JAM PDFs results in a significantly degraded χ 2 /ndf. The GG contribution for the JAM PDFs has a fall-off in x F too fast to describe the data. We observe a trend similar to the one seen already in Fig. 3: the crossover between the qq and GG terms for SMRS and GRV occurs at values of x F much larger than the ones for xFitter and JAM.

C. CERN NA3 experiment, 280 GeV/c
The CERN NA3 experiment [68], performed nearly four decades ago, still has the largest pion-induced J/ψ production statistics available today. Data were taken at three different incident momenta, 280, 200 and 150 GeV/c with both positive and negative hadron beams. The beam components were identified using Cherenkov counters. Moreover, in addition to a heavy platinum target, a liquid hydrogen target was also used, thus eliminating all possible nuclear effects. For all three energies, the cross sections have a normalization uncertainty of 13%. In the present study we only consider the NA3 hydrogen data. Unfortunately, these invaluable numerical cross sections were never published and could only be retrieved from the figures in the published paper [68] and unpublished thesis [77].
For the 280 GeV/c data taking, the authors used a 50 cm long hydrogen target, resulting in 23350 J/ψ π − events in the dimuon mass region between 2.7 and 3.5 GeV/c 2 . The retrieved data are available in 17 x F bins of 0.05, between 0.025 and 0.825. The comparison with the NLO CEM calculation is shown in Fig. 5. The resulting χ 2 /ndf values repeat The data at 200 GeV/c incident momentum were taken with a 30 cm long hydrogen target. With the negative hadron beam 3157 pion-induced J/ψ events were collected. The retrieved data extend from x F = 0.05 to x F = 0.85.
The comparison of the NLO calculation with the data is shown in Fig. 6. The agreement with the data is fair for all PDF sets, although the general trend persists: the most recent xFitter and JAM global fits have slightly worse χ 2 /ndf values. We also note that as the incident momentum decreases, the importance of the qq term increases, particularly for the larger values of x F . The GG contribution dominates the cross section for the calculation with the GRV PDFs up to x F = 0.6. In contrast, for the JAM PDFs the corresponding value is much lower, x F = 0.2.

E. CERN WA11 experiment
The WA11 Collaboration at CERN measured J/ψ production cross sections [70] using a 190 GeV/c negative pion 18 GeV/c 2 , including 7% background. The same experiment had previously measured the feed-down contribution from the χ c decays. In the cross sections shown this contribution was subtracted. For consistency, the reported feed-down contributions were added to the prompt cross section values shown, using the described procedure in reverse order. The comparison of the NLO CEM calculations with the WA11 data is shown in Fig. 7. The resulting χ 2 /ndf values are larger than for the NA3 data, pointing to additional systematic errors either in the original data or in the procedure of retrieving them. Not surprisingly however, the overall conclusions are similar to the ones made previously for the 200 GeV/c data. The calculations with SMRS and GRV are in better agreement with the data than xFitter and JAM. The comparison with the NLO CEM calculation is shown in Fig. 8. The calculated χ 2 /ndf values are rather small, pointing to somewhat overestimated experimental error bars. Nevertheless, they remain larger for the two most recent PDF sets. The overall trend previously observed is confirmed.

G. Fermilab E537 experiment
The E537 experiment at Fermilab has measured J/ψ production cross sections induced by a hadron beam of 125 GeV/c containing 82% negative pions and 18% antiprotons. Three different targets have been used: beryllium, copper and tungsten. An experimental mass resolution of σ = 200 MeV/c 2 for the Be target is reported. The 2881 collected events with the Be target in the region of the J/ψ peak cover the x F region between 0.05 and 0.95, in bins of 0.10. The normalization uncertainty on the cross sections is 6%.
The NLO CEM calculation and the E537 data are shown in Fig. 9. The χ 2 /ndf values are good for all calculations, although again slightly better for SMRS and GRV. For values of x F 0, the magnitude of the qq term is similar to that of the GG term. We also observe the relatively quick decrease of the GG term for the calculation with the JAM gluon PDF.

H. CERN WA39 experiment
The CERN WA39 Collaboration measured the J/ψ production cross section with a 39.5 GeV/c hadron beam mo- mentum. Data for the 67 cm long liquid hydrogen target were taken with negative and positive hadron beams. Measurements are reported with incident π + , π − , K + , K − , p and p. Most of the 402 events reported for the negative hadron beam are pion-induced J/ψ's. The x F -differential cross sections, available as a figure in the published paper, cover the region 0.05 ≤ x F ≤ 0.85 in bins of 0.10. The normalization uncertainty on the cross sections is 15%. The comparison between data and calculations is shown in Fig. 10. The immediate observation is that for this low incident momentum the qq contribution is much larger than the GG term, by a factor of 5-8 around x F = 0. The χ 2 /ndf values for the four PDFs are all close to 1, and slightly larger for the calculation with SMRS.

I. Observations
As a general observation, both LO and NLO CEM calculations provide a reasonable description of x F distributions of J/ψ production in the energy range considered (Figs. 3-10). We note that the large difference in the magnitude between LO and NLO is compensated by the F factor. The F factors for the xFitter and JAM PDFs are relatively stable across the range of collision energies, while the factors for SMRS and GRV PDFs show a mild rise toward low energies. From the comparison between data and calculations, interesting observations are summarized below: 1) The importance of the GG contribution relative to that of qq is greatly enhanced in the NLO calculation. As for the description of the large-x F data points for the pion beam larger than 125 GeV/c, the χ 2 /ndf values with the NLO calculations generally improve for the results with SMRS and GRV, whereas those with xFitter and JAM become worse, compared to the LO ones, for example, seeing Fig. 3.
2) At low energies, the GG contribution is relatively small but it increases rapidly with the increase of energy. The fraction of GG component is maximized around x F = 0, corresponding to the gluon density of pions at x ∼ 0.1 − 0.2. As a result of the rapid drop of pion's gluon density toward x = 1 shown in Fig. 1(c), the GG contribution decreases dramatically toward large x F . In contrast, the qq contribution falls off slower at high x F because of the relatively strong pion valence antiquark density at large x. Consequently, the qq contribution has a broader x F distribution than that of the GG contribution and the relative importance of qq rises at the large-x F region. As mentioned before, CEM dictates the relative weighting between qq and GG subprocesses by the convolution of pQCD calculation and parton densities, and the F factor cannot modify the shape of dσ/dx F . Therefore, adequate shapes of dσ/dx F distributions of individual GG and qq contributions from CEM calculations are required to achieve a reasonable description of data points at x F > 0.5. Since the partonic cross sections and nucleon PDFs are basically common in the calculations (Eq. (1)), the variation of results shall originate from the difference in the folded pion partonic densities. The calculations with SMRS and GRV pion PDFs agree with the data overall, while significantly large χ 2 /ndf values are found in the description of data with beam momentum greater than 125 GeV/c for both xFitter and JAM pion PDFs. 3) At low beam energies such as 39.5 GeV/c in Fig. 10, the qq process is the dominant mechanism of J/ψ production over the whole x F region. The data are much less, if at all, sensitive to the variation of the GG contribution. Good χ 2 /ndf values are obtained for all four pion PDFs.

V. SYSTEMATIC STUDY
Through the comparison of data with calculations over a broad energy range, we have two major findings: i) the largex F distribution of J/ψ production is sensitive to the pion gluon density; ii) the gluon densities of the recently available JAM and xFitter falls off too rapidly at large x to describe the x F distributions of J/ψ data. Judging from the consistency of observation for the data sets with proton and nuclear targets, the unaccounted nuclear medium effects such as energy loss effect is unlikely to change the conclusions.
To check the sensitivity of the CEM calculation to various QCD parameters and the choice of nuclear PDFs, we have performed a systematic study. Taking the convention of the charm quark mass in Refs. [47,48,53], we test the variations of results by setting m c to be 1.2 GeV/c 2 . The dependence on the renormalization scale µ R is checked by varying at 0.5, 1.0 an 2.0 m c [52]. We also make a different choice of nCTEQ15 [78] as the nuclear PDF in the calculations with GRV, JAM and xFitter. Overall the above observations remain qualitatively valid with respect to all these systematic variations. Figures 11 and 12 show the systematic study of compar-  11. The NLO CEM results with variation of charm quark mass mc and renormalization scale µR, compared with the dσ/dxF data of J/ψ production off the beryllium target with 515-GeV/c π − beam from the E672/E706 experiment [66]. The pion PDFs used for the calculation is GRV. The total cross sections, qq, GG and qG × (−1) contributions are denoted as black, blue, red and green lines, respectively. The charm quark mass mc, factorization scale µF and renormalization scale µR used for the CEM calculation as well as the fit χ 2 /ndf and F factors are displayed in each plot.
ing the E672/E706 data and CEM NLO calculation with GRV and JAM pion PDFs with the variation of m c and µ R . In total there are 6 settings of parameters under investigation. Overall the charm quark mass m c plays a more visible role than the renormalization scale µ R in the systematic effect. With a smaller charm quark mass m c , the fractions of qq decrease while the fractions of GG increase. The hadronization factor F drops with the decrease of m c , in accordance with a large phase space of cc production in Eq. (1). The variation of the renormalization scale µ R shows a similar but much less significant trend. For this data set at the largest beam momentum of 515 GeV/c, the GG contribution is dominating in the CEM NLO calculation. A reduction of m c from 1.5 to 1.2 GeV/c 2 reduces the relative contribution of qq and leads to a deterioration of χ 2 /ndf for both GRV and JAM. Nevertheless, this effect is particularly significant in the case of JAM. With a reduction of qq contribution, the large-x gluon density of JAM PDFs is not strong enough to sustain enough GG contribution in accounting for the cross sections at large x F . The information of F factor and χ 2 /ndf for the systematic study of this data set with the SMRS, GRV, xFitter and JAM pion PDFs, is shown in Table III.
From the systematic study of all data sets, the NLO CEM results clearly favor SMRS and GRV, especially at high energies. The χ 2 /ndf, representing the performance of data description, strongly correlates with how large the magnitude of gluon density is at the valence region. As shown in Fig. 1(c), SMRS and GRV have significantly larger gluon density at large x than xFitter and JAM. Overall our studies indicate that high-energy J/ψ data have an increased sensitivity to the pion large-x gluon density in the NLO calculations, resulting from the enhanced importance of GG contribution. On the other hand, the relatively small difference in the valence quark distributions for various PDFs plays a minor role in J/ψ production if away from the threshold region, as seen in the comparison of results of SMRS and GRV.  1   TABLE III. Results of F factor and χ 2 /ndf value of the best fit of the CEM calculations for SMRS, GRV, xFitter and JAM pion PDFs, to the data of J/ψ production off the beryllium target with 515-GeV/c π − beam [67], with the systematic variation of charm quark mass mc between 1.2 and 1.5 GeV/c 2 , and renormalization scale µR at 0.5, 1.0 an 2.0 mc.

VI. DISCUSSION
From the early CEM LO studies [3,4], it was known that the fixed-target J/ψ production is sensitive to pion valence quark distribution at low energies via qq mechanism and to the GG contribution at high energies [3]. In this study we confirm the sensitivity of the fixed-target J/ψ data to the pion's gluon density in the valence-quark region. Moreover, we show that this sensitivity is further enhanced when the NLO calculations are performed.
The NLO CEM calculations show that the x F distributions of fixed-target J/ψ production can serve as a tool for accessing the pion partonic densities. At low energies the data are predominantly sensitive to pion's valence quark distributions, while at high energies the data become increasingly sensitive to the gluon distributions in pion. Thus a global fit taking into account the J/ψ data across a broad energy range shall help to pin down the large-x gluon density of pions. At low energies where the qq mechanism dominates, the pion-induced J/ψ production, having much larger cross sections than the Drell-Yan process, could be a powerful alternative to Drell-Yan process in probing the quark distributions of pions.
We note that the recent effort to include leading neutron DIS data in the JAM global analysis has provided new constraints on pion's sea and gluon distributions at x ∼ 0.001 − 0.1 [18]. Unfortunately, the existing leading neutron DIS data are not sensitive to the PDF at x > 0.1. It is also important to include the direct photon production as well as J/ψ production data in the future global fits to place stringent constraints on the gluon distributions at large x. As shown in this study, the JAM gluon density at large x is too low to reproduce the J/ψ data. The upcoming tagged DIS (TDIS) experiment at the Jefferson Lab will be able to extend the sensitive region up to x = 0.2 [79].

VII. SUMMARY
We have examined the available pion PDFs extracted from the global fit to Drell-Yan, prompt-photon production or lead-ing neutron DIS data. These PDFs present pronounced differences, particularly in the gluon distributions. We have calculated their total and x F differential cross sections for pioninduced J/ψ production using the CEM framework at NLO. The calculations are compared to the data using the hydrogen and light nuclear targets.
We confirm the importance of the gluon fusion process in J/ψ production, especially at high (fixed-target) energies, observed in the earlier LO CEM result. We find that this dominance is even more pronounced in the NLO calculation. Since the calculated shapes of x F distributions of GG and qq contributions are directly related to the parton x distributions of corresponding PDFs, a proper description of J/ψ production data, especially for x F > 0.5 imposes a strong constrain on the relevant pion's parton PDFs. Among four pion PDFs which we have examined, the CEM NLO calculations strongly favor SMRS and GRV PDFs whose gluon densities at x > 0.1 are higher, compared with xFitter and JAM PDFs. The GG contribution from the latter two pion PDFs drops too fast toward x F = 1 to describe the data.
Our study clearly indicates that the fixed-target pioninduced J/ψ data are extremely helpful in constraining the pion gluon density, particularly at the large-x region. In the near future, new measurements of Drell-Yan as well as J/ψ data in πA reactions will be available from the CERN COMPASS experiment. While further theoretical efforts are required to reduce the model dependence in describing the J/ψ production, we believe that it is important to include the existing large amount of pion-induced J/ψ data as well as the new ones in future pion global analysis.