Contribution from Intrinsic Charm Production to Fixed-Target Interactions at the LHC

Background: Intrinsic charm, nonperturbative charm in the hadron wavefunction, has long been speculated but has never been satisfactorily proven. Open charm and $J/\psi$ measurements in a fixed-taget configuration at the LHCb searched for this contribution but reported no evidence. Purpose: $\overline D$ meson and $J/\psi$ production is calculated for the SMOG fixed-target configuration in the LHCb experiment using a combination of perturbative QCD and intrinsic charm to see whether intrinsic chaarm would indeed be observable in the SMOG kinematics. Methods: Open charm and $J/\psi$ production is calculated to next-to-leading order in perturbative QCD. Because a gas jet nuclear target is used, cold nuclear matter effects are included in the perturbative calculations. The intrinsic charm is calculated assuming production from a $|uud c \overline c \rangle$ Fock state. Results: The differential rapidity and transverse momentum distributions in $p+{\rm Ne}$, $p+{\rm He}$ and $p+{\rm Ar}$ fixed-target interactions are calculated in the SMOG acceptance and compared to data. The predicted asymmetries between $\overline D$ (leading charm) and $D$ (nonleading charm) are also shown. Conclusions: The contribution from intrinisic charm is small and decreases with center of mass energy. The calculations agree well with the current SMOG data, with or without intrinsic charm.

and Ar so far, which can interact with either proton or nuclear beams circulating in the LHC. To this point, p + He at √ s N N = 86.6 GeV and p + Ar collisions at √ s N N = 110.4 GeV [4] and p + Ne collisions at √ s N N = 68.5 GeV have been reported [5,6].
Because of the boost, the LHCb acceptance covers a range from backward rapidity to near central rapidity in the fixed-target mode, y cm from ∼ −2.5 to ∼ 0. Data were taken for proton beams of 2.55, 4 and 6.5 TeV on neon, helium and argon gas targets respectively. Data on J/ψ and D 0 production were taken in all cases. The difference in center of mass energies for the three systems resulted in a slight difference in rapidity coverage with −2.29 ≤ y cm ≤ 0 for p + Ne, −2.53 ≤ y cm ≤ 0.07 for p + He, and −2.77 ≤ y cm ≤ −0.17 for p + Ar. In all cases the transverse momentum range probed was p T ≤ 8 GeV.
The spectrometer permits the collaboration to study J/ψ through their decay to muon pairs while the siliconstrip detector makes it possible to reconstruct D 0 mesons via their two-body decays to K − π + (as well as their D 0 charge conjugates). The data reported in Ref. [4] were collected under specific beam conditions with no proton bunch crossings at the p + p interaction point to reduce background.
The J/ψ and D 0 production cross sections were obtained for p + He collisions because a luminosity determination was available for that system. However, no such determination was made for the p + Ar system. Therefore, the yields, normalized to unity, were reported instead. The p + Ne data are also reported in terms of production cross sections [5,6].
The rapidity range covered in this fixed-target setup should allow coverage up to x ∼ 0.37 for D 0 mesons, permitting a test of intrinsic charm production. In their first paper [4], the LHCb collaboration found their data to be consistent with perturbative calculations that did not include any intrinsic charm production. The calculations were performed with HELAC-ONIA [7][8][9]. The only cold nuclear matter effect included was the modification of the parton distributions in nuclei, implemented using the nCTEQ15 nuclear parton densities (nPDFs) [10].
The rapidity and p T distributions have been reported for J/ψ and D 0 /D 0 production for all systems studied [4][5][6].
In addition, the asymmetry between D 0 and D 0 mesons has been reported for p + Ne collisions [6].

III. J/ψ AND D 0 PRODUCTION IN PERTURBATIVE QCD
The open charm and charmonium production cross sections in perturbative QCD are treated similarly in this work. Because charmonium is calculated in the color evaporation model, the main difference is in the partonic center of mass energy range,ŝ, of the integration. The open charm cross section is integrated over the full available energy range ofŝ while the upper limit of theŝ integration range for J/ψ is the square of the DD mass threshold. The next-to-leading order (NLO) heavy flavor cross section is obtained using the HVQMNR code [11], both for open charm and charmonium.
More recently, the improved color evaporation model has been developed [12] which uses the quarkonium mass itself as the limit of integration and the quarkonium momentum is modified according to the quarkonium mass, originally still within the context of the HVQMNR code [12] and later to LO in k T factorization [13] and NLO in collinear factorization to study the polarization as well [14]. These changes slightly modify the slope of the J/ψ p T distribution and also result in different kinematic distributions for the J/ψ and ψ(2S) [12] although these changes are not large [12,14]. Open charm is discussed first, followed by quarkonium. The cold nuclear matter effects, which will be applied to both open charm and charmonium production, will be discussed in Sec. IV.
The perturbative open heavy flavor (OHF) cross section can be schematically represented as where ij = g + g, q + q or q(q) + g andσ ij (ŝ, µ 2 F , µ 2 R ) is the partonic cross section for initial state ij with the q(q) + g process oappearing at next-to-leading order in α s . The cross section and parton distribution functions are calculated at factorization scale µ F and renormalization scale µ R . The next-to-leading order heavy flavor cross section is obtained using the HVQMNR code [11]. Both calculations employ the same set of values for the charm quark mass, m, and scales µ F and µ R , determined from a fit to the total cc cross section at NLO in Ref. [15]: 55 −0.85 , 1.6 +0.11 −0.12 ). The scales are defined relative to the transverse mass of the pair, µ F,R ∝ m T = m 2 + p 2 T where the p T is the cc pair p T , p 2 The charm quarks become D mesons by applying a fragmentation function, D(z). The default fragmentation function in HVQMNR, applied to open heavy flavor production only, is the Peterson function [16], where z represents the fraction of the parent heavy flavor quark momentum carried by the resulting heavy flavor hadron. The default parameter ǫ P , ǫ P = 0.06, employed in HVQMNR produces softer charm meson p T distributions than supported by data, even when intrinsic parton transverse momentum is included [17], In Ref. [18] the value of ǫ P was decreased to match the FONLL D meson p T distributions with intrinsic k T included. The same procedure is followed here with the same value of ǫ P , 0.008 [18]. Note that there is some correlation in the choices of ǫ P and k 2 T . A larger ǫ P , closer to the default value of 0.06 in the HVQMNR code, results in a more steeply falling p T distribution, the effect of which could be reversed by a sufficiently large k 2 T , as discussed in Ref. [17]. However, a strong correlation only exists at low √ s N N when the average p T of the heavy quark is low and p 2 T ≈ k 2 T . Above √ s N N of a few tens of GeV, p 2 T > k 2 T and the effect of k T broadening on the p T distribution is decreased, even though k 2 T growns slowly with √ s N N . On the other hand, a fixed ǫ P reduces the charm quark momentum by the same fraction at all energies. Thus increasing √ s N N weakens the correlation between ǫ P and k 2 T . See Ref. [18] for comparisons of D and B meson distributions with different choices of ǫ P and k 2 T compared to FONLL at √ s = 7 TeV. The parton densities in Eq. (1) include intrinsic k T , required to keep the pair cross section finite as p T QQ → 0. They are assumed to factorize into the normal collinear factorization parton densities and a k T -dependent function, The CT10 proton parton distribution functions (PDFs) [19] are employed in the calculations of f p (x, µ 2 F ). Results on open heavy flavors at fixed-target energies indicated that some level of transverse momentum broadening was needed to obtain agreement with the fixed-target data once fragmentation was included [20]. Broadening has typically been modeled by intrinsic transverse momentum, k T , added to the parton densities and playing the role of low transverse momentum QCD resummation [21].
In the HVQMNR code, an intrinsic k T is added each final state charm quark, rather than to the initial state, as in the case of Drell-Yan production [21]. In the initial-state, the intrinsic k T function multiplies the parton distribution functions for both hadrons, assuming the x and k T dependencies factorize, as in Eq. (3). At leading order, there is no difference between an initial (on the partons) or final-state (on the produced charm quarks) k T kick. However, at NLO, when there is a light parton in the final state, the correspondence can be inexact. The difference between the two implementations is small if k 2 T ≤ 2 − 3 GeV 2 [20]. The rapidity distributions are independent of the intrinsic k T .
If the k T kick is not too large, it does not matter whether it is added in the initial or final state. A Gaussian distribution is employed for G p (k T ) in Eq. (3) [20], The effect of the k T kick alone hardens the single charm meson p T distribution. In Ref. [20], k 2 T p = 1 GeV 2 , in combination with fragmentation using the default Peterson parameter ǫ P , was chosen to describe low energy fixedtarget charm production.
The broadening is applied by boosting the transverse momentum of the cc pair (plus light parton at NLO) to its rest frame from the longitudinal center-of-mass frame. The transverse momenta of the incident partons, k T1 and k T2 , or, in this case, the final-state c and c, are redistributed isotropically with unit modulus, according to Eq. (4), preserving momentum conservation. Once boosted back to the initial frame, transverse momentum of the cc pair changes from p T QQ to p T QQ + k T 1 + k T 2 [17].
The broadening effect decreases as √ s increases because the perturbatively-calculated average p T of the cc pair also increases with √ s. The value of k 2 T p is assumed to increase with √ s so that effect is non-negligible for low p T heavy flavor production at higher energies. The energy dependence of k 2 T in Ref. [15] is with n = 12 for J/ψ production [15] so that k 2 T increases slowly with energy. At SMOG energies, k 2 T = 1.10, 1.12, and 1.14 GeV 2 for √ s = 68.5, 86.6 and 110.4 GeV respectively. The values of k 2 T are thus well below the limit of applicability proposed in Ref. [20].
The model calculation of J/ψ production is now described. The J/ψ production mechanism remains an unsettled question, with a number of approaches having been introduced [12,22,23]. In the calculations presented here, the Color Evaporation Model [22] is employed. This model, together with the Improved Color Evaporation Model [12], can describe the J/ψ rapidity and transverse momentum distributions, including at low p T where other approaches have some difficulties and may require a p T cut [24].
The CEM assumes that some fraction, F C , of the cc pairs produced with a pair mass below the DD pair mass threshold will go on mass shell as a J/ψ, The same mass and scale parameters are employed as in Eq. (1). However, now an upper limit of 4m 2 D is applied and the normalization factor F C is obtained by fitting the energy dependence of the J/ψ forward cross section [15].
At LO in the CEM, the J/ψ p T , p T QQ above, is zero. Thus k T broadening is required to keep the p T distribution finite as p T → 0. The intrinsic k T broadening in J/ψ production is handled the same way as for open heavy flavor production, outlined above. However, in this case, D(z) is not applied, hadronization is implied by the factor F C . (Note that, for simplicity, in the rest of this paper, when open charm meson distributions are presented, p T refers to the single charm hadron transverse momentum distribution while, when J/ψ distributions are discussed, p T refers to the transverse momentum distribution of the J/ψ.) The mass and scale uncertainties on the p + p cross sections are calculated from the one standard deviation uncertainties on (m, µ F /m T , µ R /m T ). If the central, upper and lower limits of µ R,F /m T are denoted as C, H, and L respectively, the seven sets of scale values used to calculate the uncertainty are }. At each point, the set giving the highest (lowest) value of the cross section is used to calculated the maximum (minimum) contribution to the scale uncertainty. The mass uncertainty is based on the one standard deviation uncertainty on the charm quark mass, ±0.09 GeV. The higher mass contributes to the lower limit of the cross section while the lower mass contributes to the upper limit. The uncertainty band can be obtained for the best fit sets [15] by adding the uncertainties from the mass and scale variations in quadrature. The result at each point, defines the uncertainty on the cross section. The kinematic observables, denoted by X, are y and p T in this case. The calculation labeled "cent" employs the central values of m, µ F /m T and µ R /m T . The calculations with subscript µ keep the mass fixed to the central value while the scales are varied. On the other hand, in the calculations with subscript m, the scales are fixed to their central values while the mass is varied between its upper and lower limits. Figure 1 shows the rapidity and p T distributions for D 0 production in p + p collisions at √ s N N = 68.5, 86.6, and 110.4 GeV. The central value at each energy is given by the solid curve while the limits on the uncertainty bands are given by the dashed curves. The rapidity distributions are shown in the SMOG range. The p T distributions are integrated over these rapidity ranges. The uncertainty on the rapidity distribution, shown in Fig. 1(a), at √ s = 68.5 GeV is ∼ 33%, rising to ∼ 40% at 110.4 GeV. The cross section clearly rises with energy, such that the upper limit of one uncertainty band coincides with the central cross section at the next energy. The increase in the cross section at low p T in Fig. 1 The main difference between the p T distributions at different energies is in the high p T tail which hardens with energy. Figure 2 shows the rapidity and p T distributions for J/ψ production in p + p collisions in the CEM at the same energies. Similar trends are seen as for open charm production except the uncertainty at the most central rapidity, shown in Fig. 2(a), is on the order of 60%. The J/ψ p T distributions also harden with increased energy. It is notable, however, that they are harder overall than the D meson distributions which fall off more steeply over the same p T range.

IV. COLD NUCLEAR MATTER EFFECTS
Here nuclear modifications of the parton densities, transverse momentum broadening, and J/ψ absorption in the nucleus, all cold nuclear matter effects in p + A collisions are briefly discussed. More details can be found in Ref. [3]. Only a brief discussion of each effect, specific to the SMOG setup, is included in this section.
The combined cold nuclear matter effects on perturbative QCD production of open heavy flavor and J/ψ are modified from Eqs. (1) and (6) respectively, described in Sec. III, are where The nuclear modifications of the parton distributions, R j (x 2 , µ 2 F , A) are discussed in Sec. IV A. The k T broadening in the nuclear target, G A (k T ), is discussed in Sec. IV B. Finally, J/ψ absorption by nucleons, represented by the survival probability S abs A , is described in Sec. IV C.

A. Nuclear Effects on the Parton Densities
A number of global analyses have been made to describe the modification as a function of x and factorization scale µ F , assuming collinear factorization and starting from a minimum scale, µ 0 F . Nuclear PDF (nPDF) effects generated in this scheme are generally implemented by a parameterization as a function of x, µ F and A. The k T -independent proton parton distribution functions in Eqs. (1) and (6) are replaced by the nPDFs, The NLO EPPS16 [25] nPDF parameterization is employed in these calculations for A = 4 (helium), 20 (neon) and 40 (argon). EPPS16 has 20 fit parameters, for 41 total sets: one central set and 40 error sets. The error sets are determined by varying each parameter individually within one standard deviation of its best fit value. The uncertainties on R j (x 2 , µ 2 F , A) are calculated by summing the excursions of each of the error sets from the central value in quadrature.
The nPDF uncertainties on the J/ψ and D 0 distributions are obtained by calculating the perturbative cross sections at the central values assumed for the charm mass and the factorization and renormalization scales employing the central EPPS16 set as well as the 40 error sets and summing the differences in quadrature. The resulting uncertainty bands deviate from the central cross section on the order of 20%, significantly less than the mass and scale uncertainties shown in Figs. 1 and 2.
The EPPS16 ratios for gluons are shown at the J/ψ mass scale in Fig. 3. The central sets, along with the uncertainty bands, are shown for He (A = 4), Ne (A = 20), and Ar (A = 40) targets. The EPPS16 gluon sets are shown as a function of x for all three values of A. The nPDF effects increase with A, as expected, but the effect is most significant for lower values of x than covered by the SMOG device. The x range covered by the SMOG device for near midrapidity to backward rapidity, y ∼ −2.5, is indicated by the vertical black bars. The near midrapidity value of x is near the peak of the antishadowing region, x ∼ 0.075, while the most backward rapidity reaches into the EMC region, x ∼ 0.44. The modification due to EPPS16 is smallest for A = 4 and largest for A = 40, as expected.
The ratios are all shown for the J/ψ mass, µ F = m J/ψ . At this value of the factorization scale, within a factor of three of the minimum scale employed in the EPPS16 fits, µ F = 1.3 GeV, the uncertainty band is relatively narrow in the x region spanned by SMOG. The nPDF effect for the central set is a 10-15% enhancement at y ∼ 0 and a 5-10% depletion at y ∼ −2.5. Including the range indicated by the error sets, increases the spread, particularly at negative rapidity, but the effect is still generally smaller than the mass and scale uncertainties on the p + p cross sections themselves.

B. kT Broadening
The effect and magnitude of intrinsic k T broadening on the p T distribution in p + p collisions was discussed in Sec. III. Here further broadening due to the presence of a nuclear target, single to multiple parton scatterings in the nucleus, known as the Cronin effect [26] is described. The effect is implemented by replacing The total broadening in a nucleus relative to a nucleon can be expressed as where δk 2 T is [27,28] The amount of broadening, ∆ 2 (µ), depends on the interaction scale [28] and the number of scatterings the incident proton undergoes while passing through the nucleus [3,29], For helium, neon, and argon targets, δk 2 T = 0.05, 0.15, and 0.22 GeV 2 respectively, giving an average broadening of k 2 T A = 1.17, 1.25, and 1.36 GeV 2 for the p + He, p + Ne and p + Ar systems respectively. The effect of k T broadening in nuclei, relative to the proton, is to reduce the ratio p + A/p + p at low p T and enhance it at high p T . The rapidity distributions are unaffected because they do not depend on p T .

C. Nuclear Absorption of J/ψ in p + A Interactions
In p + A collisions, the proto-J/ψ may interact with nucleons and be dissociated before it can escape the target, referred to as nuclear absorption. The effect of nuclear absorption on the J/ψ production cross section in p + A collisions may be expressed as [30] σ pA = σ pN S abs where b is the impact parameter, z is the cc production point, S abs (b) is the nuclear absorption survival probability, and σ abs (z ′ − z) is the nucleon absorption cross section. The absorption cross section here is assumed to be constant at a given energy. It is written as a function of the path length through the nucleus in Eq. (18) because other functional forms may be chosen, see e.g. Ref. [31]. The energy dependence of σ abs was studied in Ref. [32]. Using those results and following Ref. [3], values of σ abs = 4, 3.5, and 3 mb are used at √ s N N = 68.5, 86.6 and 110.4 GeV. Absorption by comoving particles, whether characterized as hadrons or partons [33] has not been included. It has the same nuclear dependence in minimum bias collisions as absorption by nucleons [34].

V. INTRINSIC CHARM
As first proposed in 1980, the proton wave function in QCD can be represented as a superposition of Fock state fluctuations, e.g. |uudg , |uudqq , |uudQQ , . . . of the |uud state. When charm quarks, Q = c, are part of the state, it is referred to as intrinsic charm, or IC. If a proton in such a state scatters in a target, the coherence of the Fock components is broken and the fluctuations can hadronize [35][36][37]. These intrinsic QQ Fock states are dominated by configurations with equal rapidity constituents, so that the heavy quarks carry a large fraction of the proton momentum [35,36]. (While proton projectiles are emphasized here, any hadron wave function can be so described.) The frame-independent probability distribution of a 5-particle IC Fock state in the proton is dP ic 5 = P 0 ic 5 N 5 dx 1 · · · dx 5 dk x 1 · · · dk x 5 dk y 1 · · · dk y 5 where i = 1, 2, 3 are the light quarks (u, u, d) and i = 4 and 5 are the c and c quarks respectively. The factor N 5 normalizes the |uudcc probability to unity and P 0 ic 5 scales the unit-normalized probability to the assumed probability of IC in the proton. The delta functions in Eq. Additional delta functions can be employed to describe hadronization by simple coalescence when the Fock state is disrupted. For example, the J/ψ x F distribution can be calculated by the addition of the delta functions, δ(x F − x c − x c ), in the longitudinal, z, direction. The summed x c and x c momentum fractions are equivalent to the x F of the J/ψ assuming that it is brought on-shell by a soft scattering with the target. Similarly, the J/ψ p T distribution is described by δ(p T − k x c − k x c )δ(k y c + k y c ) where the J/ψ p T is chosen to be along the x direction for simplicity and without loss of generality.
Likewise D mesons (D − (cd) and D 0 (cu)) mesons can be directly produced from the disrupted Fock state employing for the p T where the light parton i can be either a u or d quark. The remaining partons in the state could coalesce into a Λ c (udc) with a D 0 or a Σ ++ c (uuc) with a D − . Thus the 5-particle proton IC state could produce D 9 Λ c or D − Σ ++ c through coalescence. The J/ψ and D x F and p T distributions, integrated over all phase space, are independent of the proton energy. The p T distribution from IC only varies when phase space cuts are considered, as shown in Refs. [3,29]. The rapidity distribution, however, depends strongly on √ s because x F = (2m T / √ s) sinh y. Thus even though the x F distribution is invariant with √ s, m i and k T -integration range, the rapidity distribution is not [3]. Indeed, the x F distribution depends only weakly on the heavy quark mass [38,39].
The J/ψ and D p T distributions have long tails at large p T rather than decreasing strongly with increasing p T as the perturbative QCD calculations do [3]. This is due to the nature of Eq. (19). The requirement that the heavy quarks have higher velocities not only gives results in more forward production at low p T but also allows them to carry most of the transverse momentum at low x. This also explains why the p T distributions are suppressed at low p T when only the midrapidity (generally low x) region is considered [3]. There is thus more than one way to satisfy minimization of the denominator of Eq. (19) while still satisfying momentum conservation.
In the case of open charm, the preference for D production in a 5-particle IC state makes the D a "leading" particle relative to D + (cd) and D 0 (cu). These "non-leading" D mesons could only be produced from a 5-particle IC state by standard fragmentation and would thus be produced at lower x F or rapidity than the D mesons. In order to produce a D meson by coalescence from a IC state, as is the case for D in the |uudcc state, a higher particle-number Fock state is required, such as the seven-particle state |uudccdd , resulting in D + production [40].
Previous studies of leading D meson production, including asymmetries between D + and D − production in fixedtarget π − A interactions have shown significant differences between leading and non-leading production [41][42][43]. These asymmetries have been reproduced by IC [44] as well as by string-breaking mechanisms such as in PYTHIA [45]. Later work has shown that there can be an asymmetry in c and c distributions themselves [46]. This asymmetry arises from QCD diagrams where two gluons from two different valence quarks in the nucleon couple to a heavy quark pair, gg → QQ, with charge conjugation C = +1 [47]. This amplitude interferes with QCD diagrams where an odd number of gluons attach to the heavy quark pair, e.g. g → QQ and ggg → QQ, with C = −1. The interference of amplitudes with the same final state but different C for the QQ pair produces asymmetric distribution functions. An analogous interference term is seen in e − and e + distributions in e + e − pair production [48]. In this work, it is assumed that the asymmetry between leading and non-leading D and D meson arises from their manifestation as final-state charm hadrons from the 5-particle IC state considered here. The asymmetry due to charge conjugation effects is not taken into account.
Note that if higher Fock states are considered, only equal rapidity D and D mesons would be produced from these states and at a lower average momentum fraction for both the D and D. The probability to produces these higher Fock states would also be reduced, see e.g. Ref. [40]. Here only the 5-particle proton Fock state is considered since it gives the most forward J/ψ and D production from IC. This assumption also maximizes the possible asymmetry between D and D production from IC. The IC production cross section from the |uudcc state can be written as The factor of µ 2 /4 m 2 c is from the soft interaction which breaks the coherence of the Fock state. Here µ 2 = 0.1 GeV 2 is assumed, see Ref. [49], and an inelastic pN cross section, σ in pN = 30 mb, is employed. Although σ in pN can change slowly with √ s, it is held constant here. Equation (20) is used for open charm production, σ D ic (pp) = σ ic (pp). The J/ψ cross section from the same IC state is calculated by scaling Eq. (20) by the factor F C used in the CEM calculation in Eq. (6), The nuclear, A, dependence of the IC cross section is assumed to be that extracted for the nuclear surface-like component of J/ψ production by the NA3 Collaboration [50]. The A dependence is the same for both open charm and J/ψ, with β = 0.71 [50]. Several values of P 0 ic 5 have been employed previously. Here, a value of 1% is assumed to maximize the potential effect. The form of IC postulated by Brodsky and collaborators in Refs. [35,36], used in Eq. (19), has been adopted. Other variants of the IC distribution in the proton have been proposed, including meson-cloud models where the proton fluctuates into a D(uc)Λ c (udc) [51][52][53][54] and a sea-like distribution [55,56]. The D distribution here is similar to that in the meson-cloud model while a sea-like distribution would not produce leading D mesons and, indeed, would not result in forward charm production.
IC has been included in global analyses of the proton parton densities [55][56][57][58][59]. The range of P 0 ic 5 explored here is consistent with the upper limits of the results of these analyses. For more details of these other works, see the review of Ref. [60]. Since the discussion of Ref. [60], new work has appeared claiming evidence of IC at the 1% level [1,2], as mentioned in the introduction. This finding remains contentious [61]. (See Ref. [62] for a discussion of a possible kinematic constraint on intrinsic charm in deep-inelastic scattering.) New evidence for a finite charm quark asymmetry in the nucleon wavefunction from lattice gauge theory, consistent with IC, was presented in Ref. [46]. See also the recent review in Ref. [47] for more applications of intrinsic heavy quark states.

VI. RESULTS
In this section, the rapidity and p T distributions for J/ψ and open charm production in p+A collisions are compared to those in p + p collisions. The cold nuclear matter effects introduced in Sec. IV are included, as well as the IC contribution described in Sec. V. The calculations are compared to the SMOG data. The asymmetry between D 0 and D 0 mesons is also computed.
The cross sections for D and J/ψ production in p+p collisions including the perturbative QCD and IC contributions are: Here σ OHF (pp) and σ CEM (pp) are defined in Eqs. (1) and (6) respectively. The IC contributions can be found in Eq. (20) for D while σ J/ψ ic (pp) is given in Eq. (21). Likewise, the cross sections for D and J/ψ production in p + A collisions are: Now σ OHF (pA) and σ CEM (pA) are defined in Eqs. (9) and (10) respectively while the D and J/ψ IC cross sections are in Eqs. (22) and (23). The p+A cross sections are given for several different scenarios. While the nPDF uncertainties have been calculated for all of the p + A scenarios, only the central values are shown to avoid cluttering the figures. The uncertainties due to the nPDFs lie well within those of the p + p cross sections. The ratios of the p + A to p + p c ross sections are not presented because LHCb did not so far measure the p + p cross sections at the same energies and so do not present the data in terms of a nuclear modification factor.
The J/ψ y and p T distributions are presented in Fig. 4. The rapidity distributions are discussed first. The solid black curves show the p + A result with EPPS16 nPDF effects only. These results are generally compatible with the p + p cross section at the most negative rapidity where the modifications are small but are enhanced relative to the p + p cross section at y ∼ 0 where the momentum fraction, x, probed is in the antishadowing region. The smallest enhancement is for p + He collisions, with the lowest A, while the largest is for p + Ar collisions with the greatest value of A. The black dashed curve includes both the EPPS16 modifications and nucleon absorption. When absorption is included, suppression of the cross section relative to the central p + p cross section is observed at all energies. The dot-dashed and dotted curves show the addition of IC to the cross section. IC does not change the distribution at midrapidity but introduces a small enhancement for y < −1. The IC contribution tends to decrease with energy because the distribution is pushed to more negative rapidity with increasing energy [3]. Thus the effect should be smaller at √ s N N = 110.4 GeV than at 68.5 GeV. On the other hand, IC is more suppressed for larger mass targets so that the reduction is smallest for the He target. Nonetheless, the overall contribution from IC is small at SMOG energies and to distinguish the presence of IC, high statistics data in smaller rapidity bins at the most negative rapidity are required.
Four p + A calculations are shown for the p T dependence, also accompanied by the uncertainty band for p + p collisions including perturbative production alone. The solid curve is again the nPDF effects alone with the central value of EPPS16. A small enhancement due to antishadowing is observed for p + A relative to p + p, particularly at low p T . The smallest nPDF effect is for p + He while the largest is for p + Ar. The dashed curve includes enhanced k T broadening due to the passage of the proton through the nucleus. The effect reduces the maximum of the distribution at low p T but makes the cross section softer at higher p T . As is the case for nPDF effects, the smallest additional broadening is for the lowest mass target, thus it is smallest for p + He and largest for p + Ar. In addition, the baseline p + p intrinsic k T broadening increases slightly with √ s N N . When absorption is included, the p + A calculation shifts downward by a constant factor representative of the survival probability for that value of A.
Finally, IC is included. As previously mentioned, the IC contribution at midrapidity decreases with increasing energy so that the rapidity-integrated p T distribution is suppressed at low p T , see Ref. [3]. This suppression increases with increasing √ s N N . However, at higher p T , as the cross section calculated perturbatively begins to fall off more steeply, the long tail of the IC p T distribution appears, typically for p T > 4 GeV. This long tail is because the IC contribution to J/ψ production can take almost all of the proton momentum, especially near midrapidity. If the sum of the momenta carried by the c and c approaches the total proton momentum, the slope of the distribution will change and eventually reach an energy-dependent endpoint [3]. However, in the p T range covered by SMOG, especially with the rapidity acceptance including y ∼ 0, this energy constraint is not reached. Detection of an enhanced J/ψ distribution at high p T would be a clear signature of IC. However, high statistics data at large p T are needed, with widths no larger than 1 GeV for p T > 4 GeV. The overall agreement with the SMOG data is generally quite good. The calculated central p + A values match the data well. The data also lie within the mass and scale uncertainties of the p + p cross section. If all the uncertainties (mass, scale, and nPDF) are added in quadrature, the resulting uncertainty band would be wider still. No uncertainties have been estimated or included on ∆k 2 T or σ abs . Results for D production are shown in Fig. 5. Only two p + A curves are shown for the rapidity distribution, nPDF effects alone (solid curve) and including a 1% IC contribution. Antishadowing effects at midrapidity are apparent when compared to the p + p results. The effect of IC is again small, as expected. While the calculated rapidity distributions generally agree with the data, the calculations at the most central rapidities underestimate the SMOG data. The results are, however, still consistent with the mass and scale uncertainties on the p + pT cross section, an underestimate of the total uncertainty on the p + A cross section, as previously discussed.
Four curves are again shown for the D p T distribution. The solid and dashed black curves include only perturbative cold nuclear matter effects, showing an enhancement due to antishowing at low p T for the nPDF effects alone and show the CEM p + p calculations (no IC) at the same energy for the central value and the limits of the uncertainty band. The p + A rapidity distributions are shown for EPPS16 only (solid); EPPS16 with absorption (dashed); EPPS16 and P 0 ic 5 = 1% (dot-dashed); and EPPS16, absorption, and P 0 ic 5 = 1% (dotted). The pT distributions show EPPS16 only (solid); EPPS16 with kT kick (dashed); EPPS16, absorption, and kT kick (dot-dashed); and EPPS16, absorption, kT kick and P 0 ic 5 = 1% (dotted). The SMOG p + Ne data are from Ref. [5] while the p + He and p + Ar data are from Ref. [4]. show the NLO p + p calculations (no IC) at the same energy for the central value and the limits of the uncertainty band. The p + A rapidity distributions are shown for EPPS16 only (solid) and EPPS16 with P 0 ic 5 = 1% (dashed). The pT distributions show EPPS16 only (solid); EPPS16 with kT kick (dashed); EPPS16 and P 0 ic 5 = 1% (dot-dashed); and EPPS16, kT kick and P 0 ic 5 = 1% (dotted). The SMOG p + Ne data are from Ref. [6] while the p + He and p + Ar data are from Ref. [4]. a reduction at low p T with a harder cross section at high p T . The dot-dashed and dotted curves include a 1% IC contribution to the nPDF only and nPDF plus enhanced k T broadening respectively. The calculations without and with IC are both softer than the same distributions for J/ψ shown in the previous figure. This is because the D meson carries only a single charm quark in both cases. When IC is included, the high p T enhancement due to this effect is also reduced. The agreement with the SMOG data again, is generally good. No conclusion can be drawn about the highest p T data point because the bin width is so large. One might expect, however, that a small bin width would lie closer to the curves because the bulk of the cross section is at the lowest p T due to the steeply falling distribution.
Because the effects of IC are small on the scale of individual distributions, especially when the effect is larger at the edge of the kinematic coverage, it is worth looking at the asymmetries between distributions. Forming the asymmetry between two distributions shows the effect on a linear scale, enhancing it. There is no asymmetry for J/ψ since both the c and c are involved but there is one between D 0 (cu) and D 0 (cu). There would be a similar asymmetry between D − (cd) and D + (cd). As previously discussed, the D 0 and D − are leading charm mesons in a proton projectile because their valence quarks are shared with the valence quarks of the proton. In the context of this work, these mesons can be produced directly from the lowest 5-particle IC state of the proton. The D 0 and D + are non-leading D mesons because they would need to arise from a higher particle proton state, such as a 7-particle state and, in this case, the "leading" and "non-leading" D distributions, as defined from the 5-particle IC state would be identical. Rather than increasing a potential asymmetry, including these states would tend to dilute it, as already mentioned in Sec. V. Such asymmetries have been observed in collisions involving a π − (ud) projectile [41][42][43] where now D 0 (cu) is leading relative to D 0 (cu). (Note that D − (cd) is also now leading relative to D + (cd).) These asymmetries have not previously been clearly observed in collisions with a proton projectile, potentially because the Λ c (udc) is the natural partner to D 0 production from the 5-particle proton IC state.
The asymmetry between D 0 and D 0 production is written as where X = y or p T . The calculated asymmetry is shown in Fig. 6, both as a function of y (a) and p T (b). No uncertainties are shown because the perturbative cold nuclear matter effects, as well as the mass and scale uncertainties, cancel in the ratios as a function of rapidity. There is a slight difference as a function of p T , depending on whether or not enhanced k T broadening is included. (Note that here the asymmetry is defined as the diffreence between the leading and non-leading D 0 mesons while in the recent SMOG paper [6], the asymmetry is defined as between the c and c, giving it the opposite sign in the publication relative to this work.) The calculated asymmetry is negligible at midrapidity and increases toward backward rapidity. It is largest for p + Ne collisions at the lowest energy. The increase in the asymmetry is pushed back to more negative rapidity as the energy increases because the higher energy boosts the IC distribution to more backward rapidity, see Ref. [3]. The asymmetry in the measured range is not large, A(y) < 0.2, in all cases because the peak of the IC distribution in rapidity is at still more negative rapidity in the SMOG energy domain and the difference between the rapidity distributions with and without IC is small.
The measured asymmetry is generally larger and is finite near midrapidity. While the measurement follows the same trend as the calculation, the calculation, assuming that the asymmetry arises from IC alone, underestimates the measurement. No asymmetry is assumed to arise from the perturbative calculation itself.
There is a production asymmetry that can arise at next-to-leading order in perturbative QCD from the interference between q + g and q + g contributions to heavy quark production relative to q + q production but this is small [63]. As shown in Ref. [63], c production was enhanced over c production in π − N interactions at √ s = 23 GeV, at high x F , up to 15% at x F = 0.8, but negliglbe at x f ∼ 0. A rapidity of −2 for √ s = 68.5 GeV is equivalent to x F = −0.32. While taking these differences into account in p + A collisions here would increase the asymmetry rather than wash it out, it is very small compared to that of IC alone and would not help explain the discrepancy. Reference [63] noted that the c − c asymmetry arising from the intereference effects is separate from the previously measured leading vs. non-leading charm asymmetry [41][42][43].
On the other hand, the calculated asymmetry as a function of p T increases quickly above a few GeV. At high p T , where the IC distribution becomes harder than the perturbative distribution, A(p T ) → 1, especially for p + Ne collisions at √ s N N = 68.5 GeV. In this case, the asymmetry effectively saturates. At the higher energies, the reduction due to the further boost of the rapidity distribution to more negative rapidity delays the rise in A(p T ), particularly in p + Ar collisions. Including broadening, indicated by the dashed line denoted 'band limits' in the figure legend, slightly delays the overall increase in A(p T ) in all cases because the perturbative p T distribution is hardened. This effect is smallest for p + He collisions where A = 4. Note that while it might be possible to separate the rise in A(p T ) at different energies and target masses, especially at the largest A and highest energy, it is likely not possible to detect effects of enhanced k T broadening at high p T in the asymmetry. Likewise, the p T asymmetry is underestimated at low p T . Even if there was an enhancement of c over c at high x F or rapidity in p + A collisions, the p T distribution, integrated over the longitudinal acceptance, would not be affected. The highest p T data are relatively consistent with the calculations but this should be checked with higher statistics data and smaller bins at high p T . Measurements at lower center of mass energies by a high intensity fixed-target experiment, such as NA60+ [64], would produce an enhancement from IC at lower p T [3] and would provide a larger lever arm in energy for IC studies.
In Ref. [6], the asymmetry data at √ s N N = 68.5 GeV was also compared to calculations combining a 1% IC component with 10% recombination or coalescence which would favor D 0 production [65,66]. While this calculation gives a larger asymmetry at more negative rapidity than that shown here, they find a p T asymmetry close to zero for p T > 4 GeV. It would be interesting to test a potential coalescence effect by looking for an asymmetry between D 0 and D − production. In this case, if coalescence in the initial state were present, then for a proton, one would expect more D 0 to be produced because there are two u quarks in the proton and only one d quark. The SMOG targets are all noble gases with Z ≈ N and thus with approximately the same number of u and d quarks in the targets with a slight preference for d quarks in the Ar target. It could also be interesting to check other targets where N > Z to see if there are different asymmetries for non-isoscalar targets.

VII. CONCLUSIONS
The J/ψ and D meson distributions have been compared to fixed-target SMOG data for A = 4, 20 and 40 and √ s N N = 68.5, 86.6, and 110.4 GeV. Cold nuclear matter effects on the perturbative calculation have been included, as well as intrinsic charm. Agreement with the available SMOG data is very good overall, for both charmonium and open charm.
The results do not depend strongly on the nPDF set. While only the the p + A results for the central EPPS16 set are shown in Figs. 4 and 5, the variation in the cross sections due to the nPDF sets is on the order of 20%, smaller than the uncertainties on the p + p cross sections. The effects of IC are small and do not depend on the EPPS16 uncertainties. Choosing another nPDF set would not change this conclusion, especially in the region where IC is becoming important, at the most negative rapidities and highest p T , because the nPDFs typically differ more at low x than at high x where IC can be important. (See Ref. [67] for a discussion of nPDF variations.) The current SMOG data cannot clearly distinguish between the presence or absence of IC due to the boost of the effect to higher rapidity at increased energy. Distinguishing the effects of IC on J/ψ in the SMOG apparatus would be difficult as a function of rapidity because the effect is small in the SMOG rapidity acceptance but could be measured if SMOG could gather sufficient statistics at high p T to increase the number of p T bins above 3 GeV. This could be feasible with upgrades to SMOG and increased luminosity.
Similar statements can be made with respect to D meson measureemnts. However, studying the asymmetry between leading and non-leading D 0 mesons could help amplify the effect of IC with higher statistics data although intrinsic charm alone underestimates the measured SMOG asymmetry.