Limits on Intrinsic Charm Production from the SeaQuest Experiment

Background: A nonperturbative charm production contribution, known as intrinsic charm, has long been speculated but has never been satisfactorily proven. The SeaQuest experiment at FNAL is in an ideal kinematic region to provide evidence of $J/\psi$ production by intrinsic charm. Purpose: $J/\psi$ production in the SeaQuest kinematics is calculated with a combination of perturbative QCD and intrinsic charm to see whether the SeaQuest data can put limits on an intrinsic charm contribution. Methods: $J/\psi$ production in perturbative QCD is calculated to next-to-leading order in the cross section. Cold nuclear matter effects include nuclear modification of the parton densities, absorption by nucleons, and $p_T$ broadening by multiple scattering. The $J/\psi$ contribution from intrinsic charm is calculated assuming production from a $|uud c \overline c \rangle$ Fock state. Results: The nuclear modification factor, $R_{pA}$, is calculated as a function of $x_F$ and $p_T$ for $p+$C, $p+$Fe, and $p+$W interactions relative to $p+$d. Conclusions: The SeaQuest kinematic acceptance is ideal for testing the limits on intrinsic charm in the proton.


I. INTRODUCTION
Production of J/ψ has long been studied in a variety of environments, from more elementary collisions such as p + p and p + p, to protons on nuclei, p + A, and in nucleus-nucleus collisions, A + A. While many of the earlier studies of J/ψ production in p + A collisions were made in fixed-target configurations [1][2][3][4][5][6][7][8], more recently production has predominantly been studied at hadron and nuclear colliders, in particular at the Relativistic Heavy-Ion Collider at Brookhaven [9,10] and the Large Hadron Collider at CERN [11][12][13][14]. However, in the last several years, J/ψ production has been studied by the SeaQuest experiment at Fermilab [15] in a fixed-target setup and at a lower incident proton energy than previous experiments, allowing unprecedented forward coverage probing partons carrying a large fraction, x 1 , of the incident proton momentum.
This large x 1 coverage is of interest because it is ideal for setting limits on any cc component of the proton wavefunction, proposed by Brodsky and collaborators [16][17][18]. There have been some tantalizing hints [1,19,20] but no concrete confirmation so far. Fixed-target data usually do not cover the far forward x F region or have done so with insufficient statistics to make any definitive statement. Production at collider energies is ineffective for such studies because the collider kinematics put large x partons from the projectiles outside the detector coverage [21]. Several additional new fixed-target experiments have been proposed [22,23] or have taken data [24], including using the LHC for fixed-target studies [23,24], but all are higher energy than SeaQuest.
Earlier studies of J/ψ production at fixed-target experiments parameterized the nuclear dependence of J/ψ production as a power law [1][2][3][4][5][6][7][8], where A is the target mass and σ pA and σ pN are the production cross sections in proton-nucleus and proton-nucleon interactions respectively. If J/ψ production is unaffected by the presence of the nucleus, σ pA would grow linearly with A and α = 1. However, a less than linear A dependence has been observed [1][2][3][4][5][6][7][8]. Typical values of α were between 0.9 and 1. This A dependence was previously used to determine an effective nuclear absorption cross section [5,25,26] under the assumption that the deviation of α from unity was solely due to J/ψ dissociation by nucleons. However, more detailed studies of α as a function of the longitudinal momentum fraction, x F , and transverse momentum, p T , revealed a more complex picture. Indeed, it has long been known that α decreases as a function of x F [1][2][3]8] and increases with p T [3,6,8], neither of which can be accounted for by nuclear absorption alone. A number of cold nuclear matter effects can contribute to this dependence. These include nuclear modifications of the parton distributions in nucleons in a nucleus relative to those of a free proton [27], and energy loss or transverse momentum broadening due to multiple scattering as the projectile transits the nucleus [28][29][30], in addition to absorption. It has however, proven difficult for most model calculations to describe J/ψ production on nuclear targets as a function of x F and p T with these cold nuclear matter effects alone. The NA3 Collaboration [1] proposed dividing J/ψ production into two components that they referred to as hard and diffractive for the nuclear volume-like and surface-like dependencies respectively, The hard component is the perturbative QCD contribution discussed so far while the so-called diffractive component has been attributed to intrinsic charm, cc pairs in the proton wavefunction [16][17][18]. Since the charm quark mass is large, these intrinsic heavy quark pairs carry a significant fraction of the longitudinal momentum and contribute at large x F whereas perturbative charm production decreases strongly with x F . The light spectator quarks in the intrinsic cc state interact on the nuclear surface, leading to an approximate A 2/3 dependence [18]. The NA3 Collaboration found α ′ = 0.97 and β = 0.71 for proton projectiles [1]. This value of β is approximately consistent with a nuclear dependence due to surface interactions relative to interactions with the nuclear volume, as in hard scattering [18]. Such an approach was employed to compare to the NA3 measurements [31] for incident protons at 200 GeV and other fixed-target experiments at higher energies [32].
The SeaQuest experiment at FNAL, with a 120 GeV incident proton beam, is in an ideal kinematic regime to test intrinsic charm production. Its J/ψ acceptance is in the region 0.4 < x F < 0.95 and p T < 2.3 GeV. They have measured the x F and p T dependence of J/ψ production in p + p, p + d, p + C, p + Fe and p + W interactions [33]. A similar two-component model is employed to make predictions for these measurements, at energies and with kinematic acceptance favorable to intrinsic charm, in this work. The results are expressed through the nuclear modification, R pA , factor, the ratio of the per nucleon cross section in p + A collisions relative to p + d interactions at the same energy, instead of α, as previously employed.
Here J/ψ production in perturbative QCD is presented in Sec. II and the cold nuclear matter effects included in the calculation are introduced. The x F and p T distributions from intrinsic charm is presented in Sec. III, along with a discussion of its nuclear dependence. Section IV presents the results for the modification of J/ψ production in nuclear targets at SeaQuest, both for perturbative QCD effects alone and including the intrinsic charm contribution. The conclusions are presented in Sec. V.

II. J/ψ PRODUCTION AND COLD NUCLEAR MATTER EFFECTS IN PERTURBATIVE QCD
This section describes J/ψ production in perturbative QCD and describes how cold nuclear matter effects are implemented in this work.
The J/ψ production mechanism remains an unsettled question, with a number of approaches having been introduced [34][35][36]. In the calculations presented here, the Color Evaporation Model [34] is employed. This model, together with the Improved Color Evaporation Model [36], can describe the x F and p T distributions of J/ψ production, including at low p T where other approaches have some difficulties and may require a p T cut [37]. The CEM calculations are described in Sec. II A. The remainder of this section discusses nuclear modifications of the parton densities, nPDF effects, Sec. II B; absorption by nucleons, Sec. II C; and transverse momentum broadening, Sec. II D; are implemented in this work.

A. Charmonium Production in the Color Evaporation Model
The Color Evaporation Model [34] is employed to calculate J/ψ production. This approach assumes that some fraction, F C , of the cc pairs produced in perturbative QCD with a pair mass below that of the DD pair mass threshold will go on mass shell as a J/ψ, where ij = gg, qq or q(q)g andσ ij (ŝ, µ 2 F , µ 2 R ) is the partonic cross section for initial state ij evaluated at factorization scale µ F and renormalization scale µ R . (Note that the q(q)g process only enters at next-to-leading order in α s .) The parton densities are written here to include intrinsic k T , required to keep the pair cross section finite as p T → 0. They are assumed to factorize into the normal parton densities in collinear factorization and a k T -dependent function, The CT10 proton parton densities [38] are employed in the calculations of f p (x, µ 2 F ). The charm quark mass, factorization scale and renormalization scale determined from a fit to the total cc cross section at NLO in Ref. [39] are employed here, (m, µ F /m T , µ R /m T ) = (1.27 ± 0.09 GeV, 2.1 +2. 55 −0.85 , 1.6 +0.11 −0.12 ) where µ F is the factorization scale and µ R is the renormalization scale. The scale factors, µ F and µ R , 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 LO in the CEM, the J/ψ p T is zero. Thus, the calculated p T distribution is LO while the x F distribution is at NLO.
The Improved Color Evaporation Model has recently been developed [36] and extended to studies of quarkonium polarization [40][41][42][43] in hadroproduction. The change in the mass integration range and in the definition of the p T of the quarkonium state changes the ψ(2S) p T distribution relative to the J/ψ at high p T and, due to the narrower integration range, increases F C for J/ψ in the ICEM by ∼ 40% [36]. The J/ψ p T distributions are compatible with each other in the two approaches, see Ref. [42] for a comparison.
Calculations in the CEM were carried out at next-to-leading order (NLO) in the heavy flavor cross section using the exclusive NLO QQ production code HVQMNR [44] after imposing a cut on the pair invariant mass. This exclusive calculation is applied because it can keep track of both the c and c quarks that constitute the J/ψ. However, to keep the p T distribution finite as p T → 0, the calculations require augmentation by k T broadening, as described in more detail below.
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 after fragmentation was applied [45]. Broadening has typically been modeled by including intrinsic transverse momentum, k T , smearing to the parton densities and plays the role of low p T QCD resummation [46].
In the HVQMNR code, an intrinsic k T is added in the final state, rather than the initial state, as was previously done for Drell-Yan production [46]. When applied in the intial-state it multiplies the parton distribution functions for both hadrons, assuming the x and k T dependencies factorize in the parton distributions, as in Eq. (4). If the k T kick is not too large, it does not matter whether the k T dependence is added in the initial or final state.
Because the kick is employed in the final state, effectively the factors G p (k T ) in Eq. (4) for each parton is replaced by A Gaussian is employed to define g p (k T ) [45], In Ref. [45], k 2 T p = 1 GeV 2 was chosen to describe the p T dependence of fixed-target charm production. The broadening is applied by first boosting the produced cc pair (plus light parton at NLO) to its rest frame from its longitudinal center-of-mass frame. The intrinsic transverse momenta, k T 1 and k T 2 , of the incident partons are chosen at random with k 2 T 1 and k 2 T 2 distributed according to Eq. (6), preserving the total kick. Boosting back to the original frame changes the initial transverse momentum of the cc pair from p T to p T + k T 1 + k T 2 . If the k T was applied to the incident parton distributions, instead of the final-state cc pair, there would be no difference at LO. At NLO there can also be a light parton in the final state, making the correspondence inexact. The difference between the two implementations should be small if k 2 T ≤ 2 GeV 2 , as is the case here [45]. While the rapidity distributions are independent of the intrinsic k T , there will be some modification of the x F distribution because x F = (2m T / √ s N N ) sinh y and m T = p 2 T + m 2 . The effect of the k T kick on the p T distribution can be expected to decrease as √ s increases because the average p T also increases with energy. However, the value of k 2 T p is assumed to increase with √ s so that effect remains particularly important for low p T production at higher energies. The energy dependence of k 2 T in Ref. [39] is with n = 12 for J/ψ production [39]. At the energy of the SeaQuest experiment, √ s N N = 15.4 GeV, k 2 T p = 0.97 GeV 2 . Figure 1 shows the x F and p T distributions for J/ψ production in p + p collisions in the CEM at 120 GeV. The mass and scale uncertainties are calculated using the one standard deviation uncertainties on the quark mass and scale parameters. If the central, upper and lower limits of µ R,F /m T are denoted as C, H, and L respectively, the seven sets used to determine the scale uncertainty are {(µ F /m T , µ F /m T )} = {(C, C), (H, H), (L, L), (C, L), (L, C), (C, H), (H, C)}. The uncertainty band can be obtained for the best fit sets [39] by adding the uncertainties from the mass and scale variations in quadrature. The envelope contained by the resulting curves, At the center of mass energy of the SeaQuest experiment, the J/ψ cross section is still growing rather quickly with √ s and the mass uncertainty is dominant. Changing the charm quark mass within its 7% uncertainty can change the forward J/ψ cross section (in the range 0 < x F < 1) by almost a factor of two. On the other hand, changing the factorization scale at the relatively large momentum fraction carried by the incident partons, gives a rather small effect on the cross sections, less than 30%.
Note that in the acceptance of the SeaQuest experiment, 0.4 < x F < 0.95 and p T < 2.3 GeV, ∼ 22% of the total cross section is captured. In particular, the far forward acceptance in x F is in a region where one might expect qq or qg/qg contributions to heavy flavor production to be non-negligible compared to gg production of J/ψ. Figure 2 shows the ratio of the partial cross sections for each partonic initial state, the squared amplitude weighted by the parton densities, relative to the total cross section. Indeed at √ s N N = 15.4 GeV, if J/ψ production is broken down into its gg, qq and qg/qg contributions, the gg contribution is ≈ 80% over the x F range, integrated over all p T , while the qq and qg contributions are ∼ 25% and ∼ −10% respectively, see Fig. 2(a). The contributions from the production channels change more as a function of p T , integrated over forward x F , as shown in Fig. 2(b), As p T increases the gg and qg contributions become comparable at p T ∼ 5 GeV while the qq contribution decreases. See Refs. [44,47] for discussions of the relative contributions of the production channels. Also, note that even though the qg contribution is negative over much of the kinematic range, the total cross section is positive. These results show that in the specific kinematics and the low center of mass energy of the SeaQuest experiment, the gg contribution, while it is larger than the other contributions, is not fully dominant. At higher energies, this component is larger than 95% due to both the larger gg partonic cross section at high energies and the large gg luminosity at high energies and small x.
These distributions are not significantly modified in p + d relative to p + p collisions. The p + d cross section is about 1-2% larger per nucleon than the p + p cross section. This difference can be attributed to the modified valence quark contributions in proton relative to deuterium targets and is non-negligible at the SeaQuest energy because the qq contribution is around 25% of the cross section in the SeaQuest acceptance.

B. Nuclear Effects on the Parton Densities
The parton distribution functions in nuclei are modified from those of a free proton. A characteristic dependence of this modification has been determined from deep inelastic scattering experiments of leptons (electrons, muons, or neutrinos) on nuclei. These leptonic measurements probe only charged partons, leaving the gluon modifications to be inferred through the evolution of the parton densities with scale and momentum sum rules. The characteristic dependence is usually presented as the ratio of the structure function in a heavier nucleus, F 2A relative to the structure function in the deuteron, F 2d . The deuteron is chosen because it has an equal number of up and down valence quarks between the constitutent neutron and proton. At high momentum fractions, x > 0.3, carried by the interacting parton, there is a depletion in the nucleus relative to a nucleon, known as the EMC Effect. At smaller values of x, x < 0.03, there is also a depletion, known as shadowing. In the intermediate range of x, bridging the two regions where the ratio R = F 2A /F 2d < 1, there is an enhancement referred to as antishadowing.
A number of global analyses have been made over time by several groups to describe the modification as a function of x and factorization scale µ F , assuming collinear factorization and starting from a minimum scale, µ 0 F . These analyses have evolved with time, similar to global analyses of the proton parton distribution functions, as more data become available.
Nuclear PDF effects are generally implemented by a parameterization of the modification as a function of x, µ F and A so that the k T -independent proton parton distribution function in Eq. (3) is replaced by the nuclear parton distribution function Note that the separate modifications for valence and sea quarks in the case of up and down quarks requires the modification to be applied separately to the valence and sea quark proton distributions. The isospin also needs to be taken into account for the light quark distributions as well to distinguish between proton and neutron contributions.
Here the EPPS16 [48] nPDF parameterization at next-to-leading order in the strong coupling constant α s is employed in the calculations. It is available for A = 12 (carbon), 56 (iron) and 184 (tungsten) but assumes no modification for the deuteron target.
The EPPS16 analysis used the LHC data from the 2012 p + Pb run at √ s N N = 5.02 TeV. Data on dijet production, as well as for W ± and Z gauge boson data, are still rather sparse but probe much higher scales than previously available from fixed-target nuclear deep-inelastic scattering experiments. This analysis also included the CHORUS data [49] with neutrino and antineutrino beams. All these additional data sets gave access to a wide range of scales, 4 ≤ µ 2 F ≤ 2 × 10 4 GeV 2 , at relatively large x. They were also sensitive to differences in the valence and sea quark distributions, allowing separation between u V and d V and u and d respectively.
EPPS16 has 20 fit parameters, giving 41 total sets with 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 the individual quark ratios are calculated by summing the excursions of each of the error sets from the central value in quadrature. The sets cover the range 1.3 < Q < 10000 GeV and 10 −7 < x < 1.
The uncertainties on the J/ψ distributions due to EPPS16 are obtained by calculating the J/ψ cross sections with the central set and the 40 error sets and summing the differences in quadrature. The resulting J/ψ uncertainty band deviates from the central cross section on the order of 20%, significantly less than the mass and scale uncertainties shown in Fig. 1.
The EPPS16 ratios for gluons, valence up quarks and u quarks are shown at the J/ψ mass scale in Fig. 3. The central sets, along with the uncertainty bands, are shown for all three targets for gluons (a), anti-up quarks (b) and valuence up quarks (c). The range x > 0.007 is emphasized since SeaQuest does not probe the small x region. Indeed, the average x range probed by SeaQuest is rather narrow and in the range where antishadowing is dominant. To set the limits in x shown on Fig. 3, the momentum fraction x 2 in the nucleus is calculated for x F = 0.95, p T = 0 and x F = 0.4, p T = 2.3 GeV in 2 → 2 kinematics where x 2 = 0.5(−x F + (x 2 F + 4m 2 T /s) 1/2 ) and m T = (m 2 J/ψ + p 2 T ) 1/2 , assuming that any light final-state parton associated with J/ψ production in the 2 → 3 next-to-leading order kinematics has a small effect on the resulting x 2 . These combinations of x F and p T give the minimum and maximum values of x 2 , 0.068 and 0.136 respectively. 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 quite large.
The fixed target neutrino deep-inelastic scattering data and the gauge boson data from the LHC reveal differences between the valence and sea quark ratios. The nuclear modification is weaker for d V than for u V while there is no antishadowing for d. On the other hand, the u ratio shows antishadowing on a level similar to the valence quarks, albeit over a narrower range in x. Only the u and u V modifications are shown in Fig. 3. Examples of the sea quark and valence quark ratios are shown because SeaQuest is in a kinematic region where J/ψ production is not dominated by production via the gg channel, see Fig. 2.
Finally, it is worth noting that these modifications assume minimum bias collisions, averaging over the volume of the nucleus. Measurements of J/ψ production in p + Pb collisions at much lower x such as of the ALICE Collaboration, have been shown as a function of collision centrality, as determined by event multiplicity [50]. These results exhibit a clear dependence on where the probe impacts the nucleus, for example, in peripheral collisions, closer to the 'edge' of the nucleus where there might be only one nucleon in the path of the proton or more central collisions where the proton may encounter multiple nucleons [51,52]. It would be interesting to see whether similar modifications could be observed in high-statistics fixed-target measurements such as at SeaQuest.

C. Nuclear Absorption in pA Interactions
In p + A collisions, the proto-J/ψ produced in the initial parton interactions may further interact with nucleons and be dissociated or absorbed before it can escape the target, referred to as nuclear absorption. The effect of nuclear absorption alone on the J/ψ production cross section in p + A collisions may be expressed as [25] 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. Evem though the absorption cross section used in this work is assumed to be constant, it is written as a function of the path length through the nucleus in Eq. (12) because other functional forms have previously been assumed, see e.g. Ref. [53]. Expanding the exponent in Eq. (12), integrating, and reexponentiating the results assuming A is large leads to Eq. (1) with S abs A = A α and α = 1 − 9σ abs /(16πr 2 0 ) [25]. There have been many studies of the absorption mechanism of the J/ψ in nuclei and the cross section obtained in various model studies depends on whether the J/ψ is assumed to be produced inside or outside the target and whether it is produced as a color singlet or color octet. For comparison of results assuming different absorption scenarios, see for example Refs. [32,[53][54][55] and references therein.
The extracted absorption cross section also depends on which cold nuclear matter effects were taken into account. If one assumes only absorption, the cross section could be underestimated at low energies because antishadowing effects have been neglected. At energies where antishadowing dominates, a larger effective absorption cross section is needed to reach the level of measured suppression in p + A relative to p + p collisions, as discussed in Ref. [54]. In that reference, the effective absorption cross section was extracted for each measured x F or rapidity point for several nPDF parameterizations available at the time. While the effective absorption showed some dependence on the chosen nPDF parameterization, the authors were able to discern some trends. At midrapidity or x F ∼ 0, absorption was generally largest and was seen to decrease forward of this central region, at least up to a point. Experiments such as Fermilab's E866 experiment [6] at √ s N N = 38.8 GeV or the NA3 collaboration [1] at CERN with √ s N N = 19.4 GeV showed a strong increase in effective absorption for more forward x F , x F > 0.3 [54]. It was also noted, including by the E866 Collaboration [6], that although high x F also corresponds to low x 2 and the shadowing region, results at different energies did not scale with x 2 , demonstrating that shadowing could not be the dominant effect in that region. Other effects such as energy loss in cold matter were speculated to lead to the observed behavior.
Based on the midrapidity results, an energy-dependent effective absorption cross section was obtained employing the available data [54]. Extrapolating the findings back to 120 GeV gives 7 ≤ σ abs ≤ 9 mb. A value of 9 mb is used here and is compared to results with no nucleon absorption, σ abs = 0 and α = 1.
Even though J/ψ absorption may be inextricably entangled with other cold nuclear matter effects in hadroproduction, it may be possible to make an independent measurement of the J/ψ-nucleon total cross section, effectively the absorption cross section, from measurement of the A dependence in photoproduction. Such an experiment was performed at SLAC at 20 GeV [56], giving a value of 3.5 ± 0.8 ± 0.5 mb, lower than estimates from hadroproduction at higher energies [54]. A new experiment of the A dependence at similar energies would be useful [57], especially if the coherent and incoherent contributions could be separated. Measurements of J/ψ production in ultraperipheral collisions, while not shedding new light on absorption by nucleons, can provide additional constraints on the nuclear gluon parton densities. Measurements reported by ALICE [58,59] and CMS [60] have already put some constraints on the gluon density in the nucleus at low x.
Absorption by comoving particles, whether characterized as hadrons or partons, has been included in some calculations of cold nuclear matter effects [61]. This type of suppression has not been included separately here because it has been shown to have the same nuclear dependence in minimum bias collisions [62].

D. kT Broadening
The effect and magnitude of intrinsic k T broadening on the J/ψ p T distribution in p + p collisions was discussed in Sec. II A. Here further broadening due to the presence of a nuclear target is considered. One might expect that a higher intrinsic k T is required in a nuclear medium relative to that in p + p due to multiple scattering in the nucleus, known as the Cronin effect [63]. The effect is implemented by replacing g p (k T ) in Eq.
The total broadening in a nucleus relative to a nucleon can be expressed as The same expression for δk 2 T used in Ref. [64], based on Ref. [65], is employed here, The strength of the broadening, ∆ 2 (µ), depends on the interaction scale [65], The scale at which the broadening is evaluated is taken to be µ = 2m c [64]. The scale dependence suggests that one might expect a larger k T kick in the nucleus for bottom quarks than charm quarks. The size of the effect depends on the number of scatterings the incident proton could undergo while passing through a nucleus. This is given by ν − 1, the number of collisions, less the first collision. The number of scatterings strongly depends on the impact parameter of the proton with respect to the nucleus. However, in minimum-bias collisions, the impact parameter is averaged over all possible paths. The average number of scatterings is where and ρ A is the nuclear density distribution. An average nuclear density, ρ 0 = 0.16/fm 3 , is assumed for a spherical nucleus of radius R A = 1.2A 1/3 ; and σ in pp is the inelastic p + p cross section, ∼ 32 mb at fixed-target energies.
With the given values of ρ 0 and σ in pp in Eq. (16) and ∆ 2 (µ = 2m c ) = 0.101 GeV 2 from Eq. (15),, For carbon, iron and tungsten targets, δk 2 T = 0.1, 0.25, and 0.39 GeV 2 respectively. This gives an average broadening in the nucleus of k 2 T A = 1.07, 1.22, and 1.36 GeV 2 for the SeaQuest targets. The change in k 2 T between carbon and tungsten targets is similar to that observed for J/ψ production in other fixed-target experiments [8]. That work also observed that this relative broadening was independent of center-of-mass energy.
The effect of k T broadening in nuclei, relative to the proton, is to reduce the ratio p + A/p + p or p + A/p + d at low p T and enhance it at high p T . Because x F depends on m T , one can expect a small effect as a function of x F while a more significant one should be seen as a function of p T .

III. INTRINSIC CHARM
The wave function of a proton in QCD can be represented as a superposition of Fock state fluctuations, e.g. |uudg , |uudqq , |uudQQ , . . . of the |uud state. When the projectile scatters in the target, the coherence of the Fock components is broken and the fluctuations can hadronize [16,17,66]. These intrinsic QQ Fock states are dominated by configurations with equal rapidity constituents. Thus the intrinsic heavy quarks carry a large fraction of the projectile momentum [16,17].
The frame-independent probability distribution of a 5-particle cc Fock state in the proton is where i = 1, 2, 3 are the interchangeable light quarks (u, u, d) and i = 4 and 5 are the c and c quarks respectively. Here N 5 normalizes the |uudcc probability to unity and P 0 ic 5 scales the unit-normalized probability to the assumed intrinsic charm content of the proton. The delta functions conserve longitudinal (z) and transverse (x and y) momentum. The denominator of Eq. 18 is minimized when the heaviest constituents carry the largest fraction of the longitudinal momentum, x Q > x q .
Note that only this 5-particle Fock state of the proton is considered here since it gives the most forward x F production of J/ψ from intrinsic charm. One can also consider higher Fock components such as |uudccqq but these will reduce the average momentum fraction of the J/ψ and also enter at lower probabilities, see e.g. Refs. [67] for examples of charm hadron distributions from higher Fock states.
Heretofore, only the x F dependence of the J/ψ from the 5-particle Fock state has been considered. In this case average values for the transverse masses of the constituent quarks, m 2 i = m 2 i +k 2 T i , can be used. Previously, m q = 0.45 GeV and m c = 1.8 GeV were chosen [31]. In this case, the J/ψ x F distribution can be calculated assuming simple coalescence of the c and c in a single state, represented in Eq. (18) by the addition of a delta function, δ(x F − x c − x c , for the longitudinal, z, direction. We can related the summed x c and x c momentum fractions to the x F of the J/ψ assuming that the it is brought on-shell by a soft scattering with the target. In this work, the p T distribution of the J/ψ from such a Fock state is calculated explicitly for the first time. Thus the k x and k y integrals are also represented in Eq. (18). To simplify the calculations, it is assumed that the transverse momentum of the J/ψ, p T , is in the x direction only. This introduces two delta functions to Eq. (18), δ(p T −k x c −k x c ) and δ(k y c + k y c ).
Employing these delta functions, the longitudinal momentum fraction of the c can be integrated, as well as the transverse momentum of the c quark in the proton, leaving the light quark x and k as well as that of the c quark to be integrated numerically.
The J/ψ x F and p T probability distributions from intrinsic charm are shown in Fig. 4. The distributions are normalized to unity. How this normalized probability distribution is converted into a cross section for J/ψ production by intrinsic charm is explained shortly. Results are shown for different values of k T for the light and charm quarks in the proton.
As shown in Fig. 4(a), the x F dependence is effectively independent of the chosen k T limits. The average x F of the distribution is 0.53, almost exactly the same as previous calculations assuming constant values ofm q = 0.45 GeV andm 2 c = 1.8 GeV [31]. The charm quark x c distribution obtained by integrating over the momenta of the other four partons in the Fock state, assuming the quark is liberated from the projectile proton by a soft interaction. It has an average value of x c = 0.34, slightly larger than half that of the J/ψ.
The p T distributions, on the other hand, show a significant dependence on the range of k T integration. The chosen default values are k max q = 0.2 GeV and k max c = 1.0 GeV, shown by the red curve in Fig. 4(b). Note that these limits apply to both the x and y directions of the transverse momentum. The charm quark distribution, with these same default values for the range of integration, is given by the black curve.
If the limits of the integration range are doubled to k max q = 0.4 GeV and k max c = 2.0 GeV, the p T distribution is broadened, see the blue dashed curve in Fig. 4(b), while halving the integration range to 0.1 GeV and 0.5 GeV for the light and charm quarks respectively, results in the dot-dashed magenta curve. The average J/ψ p T resulting from these values of k max i are 2.06, 2.48, and 1.90 GeV for the default range, doubling the range, and halving the range, respectively. The average p T of the charm quark in the default integration range is 2.00 GeV.
It is notable that the J/ψ p T distribution from intrinsic charm is considerably broader than that of the J/ψ from the CEM shown in Fig. 1(b). This can be easily understood, even when an average k T 2 ∼ 1 GeV 2 is employed with the CEM, when one considers the two sources.
In perturbative QCD, the p T range depends on the center of mass energy with a limit of p T ∼ √ s N N /2 at x F = 0 for massless partons and less for massive quarks, where p T is replaced by m T . Leading order 2 → 2 scattering is also assumed for this estimate, leading to p T ≤ 7.2 GeV in the massless case and p T ≤ 6.5 GeV for the J/ψ in the CEM. In addition, since the initial partons taking part in the interaction come from two different hadrons, one from the projectile and the other from the target, at least one of them will carry a much smaller fraction of the momentum that the intrinsic charm quarks in the proton Fock state.
On the other hand, when the J/ψ arises from an intrinsic charm state of the proton, according to Eq. (18), there is no energy limit on the p T distribution other than that imposed by momentum conservation. The J/ψ kinematics all come from the incident proton. Only a soft interaction with the target is sufficient to disrupt the Fock state and bring the J/ψ on mass shell.
Even though intrinsic charm may dominate at higher p T for the SeaQuest energy and kinematics, at higher energies, the J/ψ cross section from perturbative QCD will be significantly higher and the average p T of the CEM calculation grows with energy, along with the intrinsic k T kick imparted to the J/ψ at low p T , see Eq. (7). In particular, at colliders, the acceptance is generally around midrapidity (low x F ) and the largest part of the intrinsic charm contribution is thus outside the detector acceptance and detection of intrinsic charm is less favorable.
The intrinsic charm production cross section from the |uudcc component of the proton can be written as The factor of µ 2 /4 m 2 c arises from the soft interaction which breaks the coherence of the Fock state where µ 2 = 0.1 GeV 2 is assumed, see Ref. [31], and σ in pN = 30 mb is appropriate for the SeaQuest energy. Employing this cross section for charm production from the Fock state, the J/ψ contribution is assumed here to be calculated by scaling Eq. (19) by the same factor, F C , used to scale the CEM calculation to obtain the inclusive J/ψ cross section from the cc cross section in Eq. (3), The nuclear dependence of the intrinsic charm contribution is assumed to be the same as that extracted for the nuclear surface-like component of J/ψ dependence by the NA3 Collaboration [1] in Eq. (2), so that with β = 0.71 [1] for a proton beam.
To represent the uncertainties on intrinsic charm, several values of the intrinsic charm probability, P 0 ic 5 , are employed. The EMC charm structure function data is consistent with P 0 ic 5 = 0.31% for low energy virtual photons but P 0 ic 5 could be as large as 1% for the highest virtual photon energies [19,68]. For a lower limit, a probability of 0.1% is used. All three results are generally shown in Sec. IV unless otherwise noted.
In this work, the formulation for intrinsic charm in the proton wavefunction postulated by Brodsky and collaborators in Refs. [16,17] has been adapted. Other variants of the intrinsic charm distribution in the proton exist, including meson-cloud models where the proton fluctuates into a D(uc)Λ c (udc) state [69][70][71][72] and a sea-like distribution [73,74]. Intrinsic charm has also been included in global analyses of the parton densities [73][74][75][76][77]. The range of P 0 ic 5 explored here is consistent with the results of the global analyses. For more details of these other works, see the review of Ref. [78]. New constraints from lattice gauge theory are also available in Ref. [79]. See also the recent review of Ref. [80] for more applications of intrinsic heavy quark states.

IV. RESULTS
This section presents the results for the nuclear suppression factor, where the cross section per nucleon in a nuclear target is compared to that of the deuteron target. The relative per nucleon cross sections are a more straightforward way to present the nuclear dependence than using the exponent α as was done previously. Using the deuteron target reduces the potential, albeit small, isospin dependence of the ratios. The calculations shown in this section are all in the x F and p T bins defined by SeaQuest [33]. The calculated points shown in the figures are placed at the arithmetic center of the bin rather than a cross-section weighted center. Changing the location of the calculated points within the bins could slightly modifiy the shape of the x F and p T dependence in this section but not enough to alter any conclusions.
Section IV A shows the results with cold nuclear matter effects on the perturbative QCD contribution alone while Sec. IV B includes intrinsic charm, both without and with nuclear absorption on the pQCD contribution. Results are calculated both as a function of x F and p T and are shown for all three nuclear targets: A = 12 (carbon); A = 56 (iron); and A = 184 (tungsten).

A. RpA From Perturbative QCD Effects Alone
The focus of the discussion in this section is on only the nuclear effects described in Sec. II. In this case, the cross sections employed in the calculation of R pA are where In the case of the deuteron target, R j ≡ 1 in EPPS16, and δk 2 T = 0 in Eq. (14) so that g A (k T ) = g p (k T ). It is also assumed that σ abs = 0 since the proton and neutron in the deuteron are far separated so the incident proton is likely to interact with only one of the two nucleons in the deuteron. Thus σ pd = 2σ CEM (pp). As was the case in Eq. (5), here the broadening is applied as discussed in Sec. II D. Results are shown in Fig. 5 as a function of x F and in Fig. 6 as a function of p T for effects on the three nuclear targets. The effects are added sequentially in these figures. Only nPDF effects are common to all the calculated ratios. These nPDF results are presented with the central, best fit, set given by the solid curves while the uncertainties added in quadrature are outlined by the dashed curves. First nPDF effects alone are shown by the red points and curves. In this case, g A (k T ) = g p (k T ), no additional broadening is taken into account due to the nuclear target. The magenta points and curves show the effect of enhanced k T broadening in the nuclear target. Next, absorption is added in the two cases with nPDF effects plus absorption shown in the blue points and curves while all three effects: nPDF, absorption and enhanced k T broadening are shown by the cyan points and curves.
The results with nPDF effects alone are generally independent of x F and p T and all show R pA > 1. As seen in Fig. 3, the x F and p T acceptance of the SeaQuest experiment restricts the x 2 range in the nuclear targets to a rather narrow window close to x 2 ∼ 0.1 where antishadowing is dominant. The low scale of the calculation, not much greater than the minimum scale of EPPS16, indicates that the uncertainties will be larger than at higher scales. The EPPS16 nPDF uncertainties are large. The effects grow in magnitude as the target mass is increased from A = 12 to A = 184, as do the widths of the uncertainty bands. It is worth noting that only the lower limit of EPPS16 band gives R pA < 1 in this case.
The effect of k T broadening in the nucleus relative to the deuteron changes the shape of the p T distribution, particularly at low p T . By adding the δk 2 T due to multiple scattering to the intrinsic k T needed for J/ψ production in p + p collisions, the average p T is increased in p + A collisions. The tails of the p T distributions are typically less affected, particularly at collider energies where the enhanced k T broadening in a nuclear has to be large to produce a substantial effect because the added δk 2 T is small compared to p T of the J/ψ. In addition, if k 2 T

1/2
A is much less than the mass of the particle to which it is applied the effect is reduced. For example, the effect of k T broadening, even in p + p, is much smaller for bottom quarks than for charm quarks, see Refs. [81,82].
Because the mass scale of the J/ψ is larger than the applied k T kick, the low center of mass energy of the SeaQuest experiment results in a larger effect than at collider energies and even than at other fixed-target energies. This effect is particularly visible on the p T dependence of the suppression factor shown in Fig. 6. The broadening results in a substantial increase in R pA with p T , especially for heavier targets where δk 2 T is a factor of 2-3 larger than the change in carbon.
It is notable, however, that there is also a slight change in the x F dependence of R pA due to broadening since x F = 2m T sinh y/ √ s N N . This change is significantly smaller than as a function of p T because x F depends instead on m T = m 2 + p 2 T and m and p T are of comparable magnitude over the measured p T range. Finally, there would be no change due to k T broadening if the results were given as a function of rapidity instead.
Rather strong absorption is required to nullify the effects of antishadowing and produce R pA < 1, as seen when absorption is added to the calculation. Since a constant 9 mb cross section has been assumed, as inferred for midrapidity (x F ∼ 0) in Ref. [54], absorption does not change the dependence of R pA on x F and p T . The effective exponent α obtained with σ abs = 9 mb is 0.88.
Reference [54] also showed that σ abs could increase significantly with x F for x F > 0.25, resulting in a decrease of R pA with x F . In that reference and in other works, see for example Ref. [28,30,32,83], the decrease with x F was associated with energy loss, either in the initial state or the final state. Here it is argued that this behavior could be attributed to intrinsic charm. This possibility will be tested in the next section where the effects on an intrinsic charm contribution is studied.
In Ref. [54], away from midrapidity but for x F ≤ 0.25, the effective absorption cross section could be interpreted as starting to decrease with x F , especially for higher fixed-target energies. This could be attributed, for example, to the J/ψ being formed outside the target. Given the forward acceptance of the SeaQuest experiment, well forward of the acceptances of most of the experiments included in Ref. [54], it is possible that a lower value of σ abs could be extracted from comparison with the SeaQuest data.

B. RpA Including Intrinsic Charm
Intrinsic charm is added to the calculations of R pA here. In this case, Here σ CEM (pA) was defined in Eq. (23) while σ J/ψ ic (pA) is given in Eq. (21). Note that, unlike the nuclear volume-like mass dependence of the perturbative QCD effects, the nuclear surface-like A dependence is assumed to apply to the deuteron. Later in this section, this assumption will be relaxed to check the dependence of the results on this assumption. Three values of P 0 ic 5 in Eq. (18) are shown in the following four figures: 0.1%, 0.31% and 1%. This range can be taken as an uncertainty band on intrinsic charm.
As in the previous section, the effects are added sequentially in these figures. Now nPDF effects and intrinsic charm are common to all the calculated ratios. The EPPS16 results are presented with the central, best fit, set given by the solid curves while the uncertainties added in quadrature are outlined by the dashed curves. Those calculations, affecting σ CEM are added to the intrinsic charm cross section, as in Eqs. (27) and (28). Results are shown without absorption, σ abs = 0, in Figs. 7 and 8 while absorption is included with σ abs = 9 mb in Figs. 9 and 10.
In these figures, the red, blue and black curves show the nPDF effects with P 0 ic 5 = 0.1%, 0.31% and 1% respectively in Figs. 7 and 8. The solid curves show the EPPS16 central value while the dashed curves outline the uncertainty band. For these calculations, g A (k T ) = g p (k T ). The magenta, cyan and green solid and dashed curves include enhanced k T broadening in the nuclear target for P 0 ic 5 = 0.1%, 0.31% and 1% respectively. The same color scheme is used in Figs. 9 and 10 but with absorption included in σ CEM (pA).
The nuclear dependence of intrinsic charm is not shown separately because no other nuclear effects are active with this component, only the A β dependence shown in Eq. (21). Thus R pA for intrinsic charm alone would be independent of x F and p T . Since intrinsic charm is included only in the initial proton and the J/ψ comes on shell without a hard interaction with the target, there are no nPDF effects. For the same reason, the intrinsic charm contribution is unaffected by multiple scattering of the proton in the target nucleus.
The effect of the intrinsic charm contribution is immediately clear in Fig. 7: R pA < 1 without nuclear absorption and is now a decreasing function of x F . It is also clear that, at least for the SeaQuest energy and in its forward acceptance, the intrinsic charm cross section is comparable to the perturbative QCD cross section. As P 0 ic 5 increases from 0.1% to 1%, the intrinsic charm contribution comes to dominate the nuclear dependence as a function of x F . The most immediate qualitative measure of this dominance is apparent in the width of the EPPS16 uncertainty band. These bands, are wide in Figs. 5 and 6 without intrinsic charm, narrower but still relatively far apart for P 0 ic 5 = 0.1%, and barely distinguishable for P 0 ic 5 = 1% because, in this case, the nuclear dependence is wholly dominated by intrinsic charm.
There is some notable separation between the different assumed values of P 0 ic 5 without absorption included. The difference due to the k T broadening, apparent for the red and magenta curves in Fig. 7, especially for the iron and tungsten targets, is almost indistinguishable for P 0 ic = 1%, in the black and green curves. The dominance of the intrinsic charm contribution as P 0 ic 5 increases also becomes obvious in the weakening of the x F dependence with increasing P 0 ic 5 . When a 1% probability is assumed, the results are almost independent of x F except in the lowest x F bins, as would be expected if J/ψ production was via the intrinsic charm contribution alone.
Similar behavior is observed as a function of p T in Fig. 8. Even without an enhanced k T broadening in the nucleus, a rise in R pA with p T is still observed, although it is not as strong a function of p T without broadening as it is with it. In this case, with the same k T kick in p + A as in p + p, there is an increase compared to p T → 0 which then levels off in the last two, larger p T , bins. On the other hand, when broadening is included, the rise in R pA with p T continues over the entire p T range shown. The perturbative QCD and intrinsic charm contributions are thus competitive with each other as long as P 0 ic 5 is low, 0.31% or less, but not for P 0 ic 5 = 1% where the distinction between calculations with and without broadening is negligible.
Results including absorption are given in Figs. 9 and 10. The strong absorption cross section employed in these calculations, σ abs = 9 mb, results in dominance of the intrinsic charm contribution even when P 0 ic 5 is as low as 0.1%. Note, however, that the dependence of R pA with x F and p T is not significantly affected by an assumed 1% probability for intrinsic charm because, for this value of P 0 ic 5 , the intrinsic charm contribution is already dominant with σ abs = 0, the higher probability overcomes the stronger surface nuclear target dependence of intrinsic charm. Now, with absorption included, the A dependence of the two contributions are similar in magnitude since S abs A ≈ A α = A 0.88 , see the discussion under Eq. (12) in Sec. II C. Even though α is still greater than β = 0.71, assumed for intrinsic charm [1,31], the two values are more comparable.
The calculations including absorption further compresses the pQCD uncertainties due to nPDF effects and k T broadening. The differences between the calculations with nPDF effects only and nPDFs with k T broadening on the FIG. 7: (Color online) The nuclear modification factors for J/ψ production in SeaQuest as a function of xF for the combined pQCD and intrinsic charm cross section ratios for carbon (a), iron (b) and tungsten (c) targets relative to deuterium. All calculations are shown with σ abs = 0. Results with nPDF effects and the same kT in p + d and p + A are shown in the red, blue and black curves while nPDF effects with an enhanced kT kick in the nucleus are shown in the magenta, cyan and green curves. The probability for IC production is 0.1% in the red and magenta curves; 0.31% in the blue and cyan curves; and 1% in the black and green curves. The solid lines show the results with the central EPPS16 set while the dashed curves denote the limits of adding the EPPS16 uncertainities in quadrature.
perturbative part are almost indistinguishable as a function of x F except for the heaviest targets, as seen in Fig. 9. The same can be seen as a function of p T in Fig. 10. In this figure, the differences between the results without and with k T broadening are only visible for P 0 ic 5 = 0.1%. The last part of this section further tests the assumptions made about how intrinsic charm is implemented. First, the assumption regarding whether one treats the deuterium target as a nucleus or like a proton is tested, with results shown in Fig. 11. Next, the influence of the range of k T integration on the intrinsic charm p T distribution is checked in Fig. 12.
When calculating R pA with a deuteron target, it has so far been assumed that the surface-like A dependence also applies to the deuteron, unlike the assumption that σ abs = 0, a pure volume-like A dependence, for perturbative QCD effects. This has the effect of reducing the intrinsic charm contribution in p + d interactions relative to p + A for heavier nuclei. In Fig. 11, it is assumed that, also in the case of intrinsic charm, the deutron target is effectively the same as a proton target, with β = 1 in Eq. (21). Instead of showing all values of P 0 ic 5 , only calculations with P 0 ic 5 = 0.31% are presented.
Results are given for all three targets with EPPS16 and absorption. Calculations both with and without an intrinsic k T enhancement are presented. The R pA ratios are, from top to bottom, p + C (red without k T enhancement and magenta with); p + Fe (blue and cyan); and p + W (black and green). Making the intrinsic charm contribution in p + d the same as in p + p has the effect of reducing the magnitude of the ratios by a few percent without changing the shape of R pA as a function of x F or p T . The difference would be slightly larger for a smaller value of P 0 ic 5 but would not change the overall result significantly. Thus, at least from the point of view of statistics with a proton or FIG. 8: (Color online) The nuclear modification factors for J/ψ production in SeaQuest as a function of pT for the combined pQCD and intrinsic charm cross section ratios for carbon (a), iron (b) and tungsten (c) targets relative to deuterium. All calculations are shown with σ abs = 0. Results with nPDF effects and the same kT in p + d and p + A are shown in the red, blue and black curves while nPDFs with an enhanced kT kick in the nucleus are shown in the magenta, cyan and green curves. The probability for IC production is 0.1% in the red and magenta curves; 0.31% in the blue and cyan curves; and 1% in the black and green curves. The solid lines shown the results with the central EPPS16 set while the dashed curves denote the limits of adding the EPPS16 uncertainities in quadrature. deuteron target, there is no reason, calculationally, to prefer one over the other.
To end this section, Fig. 12 shows how R pA could be modified as a function of p T if a wider range of k T integration for intrinsic charm had been chosen. In particular, k max q = 0.4 GeV and k max c = 2 GeV are assumed since this choice resulted in a broader J/ψ distribution from intrinsic charm with a correspondingly reduced peak of the p T distribution at low p T , in the range of the SeaQuest acceptance. Results are not also shown for the lower values because, although the p T distribution becomes somewhat narrower, the effect is not large enough to modify the results. The nuclear suppression factor is only shown as a function of p T because the range of k T integration has no effect on the x F distribution of the produced J/ψ, see Fig. 4.
Similar to the results shown in Fig. 11, only P 0 ic 5 = 0.31% is presented for all three targets, without and with enhanced k T broadening, in addition to nPDF effects and absorption. The broader intrinsic charm p T distribution does not change the overall suppression of R pA significantly but it does show a larger difference between calculations with additional k T broadening and those without.

V. CONCLUSIONS
The low center of mass energy of the SeaQuest experiment, as well as its forward acceptance, make it an ideal environment for probing the existence of an intrinsic charm contribution to J/ψ production. The center of mass energy is a factor of 3.1 above the J/ψ production threshold. This relatively low energy makes the perturbative QCD FIG. 9: (Color online) The nuclear modification factors for J/ψ production in SeaQuest as a function of xF for the combined pQCD and intrinsic charm cross section ratios for carbon (a), iron (b) and tungsten (c) targets relative to deuterium. All calculations are shown with nuclear absorption included. Results with EPPS16 and the same kT in p + d and p + A are shown in the red, blue and black curves while EPPS16 with an enhanced kT kick in the nucleus are shown in the magenta, cyan and green curves. The probability for IC production is 0.1% in the red and magenta curves; 0.31% in the blue and cyan curves; and 1% in the black and green curves. The solid lines shown the results with the central EPPS16 set while the dashed curves denote the limits of adding the EPPS16 uncertainities in quadrature. cross section compatible with or less than the J/ψ cross section from intrinsic charm, depending on the experimental x F range. Higher center of mass energies increase the J/ψ cross section in the CEM dramatically, see, e.g. Ref. [39], while the intrinsic charm contribution grows more slowly, depending only on σ in pN , see Eq. (19). In addition, the high x F range covered by SeaQuest is exactly the region where intrinsic charm should dominate production. As seen in a comparison of the x F distributions in Figs. 1(a) and 4(a), the CEM cross section decreases by an order of magnitude over the SeaQuest x F acceptance, with a maximum at x F = 0, outside the SeaQuest acceptance. On the other hand, the peak of intrinsic charm probability distribution is at x F ≈ 0.53, within the range of the SeaQuest measurement.
A comparison of the SeaQuest J/ψ production data on its nuclear targets could set limits on σ abs in the perturbative QCD contribution and P 0 ic 5 , the probability of the intrinsic charm contribution in the proton.
FIG . 11: (Color online) The nuclear modification factors for J/ψ production in SeaQuest as a function of xF (a) and pT (b) for the combined pQCD and intrinsic charm cross section ratios for carbon (red and magenta curves), iron (blue and cyan curves) and tungsten (black and green curves) targets relative to a proton target. All calculations are shown with nuclear absorption included. Results with EPPS16 and the same kT in p + p and p + A are shown in the red, blue and black curves while EPPS16 with an enhanced kT kick in the nucleus are shown in the magenta, cyan and green curves. The probability for IC production is 0.31% in all cases. The solid lines shown the results with the central EPPS16 set while the dashed curves denote the limits of adding the EPPS16 uncertainities in quadrature.
FIG. 12: (Color online) The nuclear modification factors for J/ψ production in SeaQuest as a function of pT for the combined pQCD and intrinsic charm cross section ratios for carbon (red and magenta curves), iron (blue and cyan curves) and tungsten (black and green curves) targets relative to a proton target. All calculations are shown with nuclear absorption included. Results with EPPS16 and the same kT in p + d and p + A are shown in the red, blue and black curves while EPPS16 with an enhanced kT kick in the nucleus are shown in the magenta, cyan and green curves. The probability for IC production is 0.31% in all cases. A broader IC distribution is assumed here, corresponding to the blue dashed curve in Fig. 4(b). The solid lines shown the results with the central EPPS16 set while the dashed curves denote the limits of adding the EPPS16 uncertainities in quadrature.