Measurement of $J/\psi$ at forward and backward rapidity in $p$+$p$, $p$$+A$l, $p$$+A$u, and $^3$He+Au collisions at $\sqrt{s_{_{NN}}}=200~{\rm GeV}$

Charmonium is a valuable probe in heavy-ion collisions to study the properties of the quark gluon plasma, and is also an interesting probe in small collision systems to study cold nuclear matter effects, which are also present in large collision systems. With the recent observations of collective behavior of produced particles in small system collisions, measurements of the modification of charmonium in small systems have become increasingly relevant. We present the results of $J/\psi$ measurements at forward and backward rapidity in various small collision systems, $p$$+$$p$, $p$$+$Al, $p$$+$Au and $^3$He$+$Au, at $\sqrt{s_{_{NN}}}$=200 GeV. The results are presented in the form of the observable $R_{AB}$, the nuclear modification factor, a measure of the ratio of the $J/\psi$ invariant yield compared to the scaled yield in $p$$+$$p$ collisions. We examine the rapidity, transverse momentum, and collision centrality dependence of nuclear effects on $J/\psi$ production with different projectile sizes $p$ and $^3$He, and different target sizes Al and Au. The modification is found to be strongly dependent on the target size, but to be very similar for $p$$+$Au and $^{3}$He$+$Au. However, for 0%--20% central collisions at backward rapidity, the modification for $^{3}$He$+$Au is found to be smaller than that for $p$$+$Au, with a mean fit to the ratio of $0.89\pm0.03$(stat)${\pm}0.08$(syst), possibly indicating final state effects due to the larger projectile size.

Charmonium is a valuable probe in heavy-ion collisions to study the properties of the quark gluon plasma, and is also an interesting probe in small collision systems to study cold nuclear matter effects, which are also present in large collision systems. With the recent observations of collective behavior of produced particles in small system collisions, measurements of the modification of charmonium in small systems have become increasingly relevant. We present the results of J/ψ measurements at forward and backward rapidity in various small collision systems, p+p, p+Al, p+Au and 3 He+Au, at √ s N N =200 GeV. The results are presented in the form of the observable RAB, the nuclear modification factor, a measure of the ratio of the J/ψ invariant yield compared to the scaled yield in p+p collisions. We examine the rapidity, transverse momentum, and collision centrality dependence of nuclear effects on J/ψ production with different projectile sizes p and 3 He, and different target sizes Al and Au. The modification is found to be strongly dependent on the target size, but to be very similar for p+Au and 3 He+Au. However, for 0%-20% central collisions at backward rapidity, the modification for 3 He+Au is found to be smaller than that for p+Au, with a mean fit to the ratio of 0.89 ± 0.03(stat)±0.08(syst), possibly indicating final state effects due to the larger projectile size.

I. INTRODUCTION
The cross section for production of charmonium in proton collisions with heavy nuclei is strongly modified relative to that in p+p collisions. The effects that cause this modification are often referred to as cold nuclear matter (CNM) effects because of the long-time presumption that the energy density and temperature produced in the collision of a single proton with a nucleus were not sufficient to form a deconfined quark-gluon plasma. The CNM effects that can modify charm production in p+A collisions include modification of the nuclear-parton-distribution functions (nPDFs) in a nucleus [1,2], initial state parton energy loss [3], breakup of the forming charmonium in collisions with target nucleons [4,5], coherent gluon saturation [6,7], and transverse momentum broadening [8]. These mechanisms are generally expected to act in the early stages of the collision, and effect either the production rates of charm quarks or their propagation through the nucleus. All of these processes are strongly (and differently) dependent on the rapidity and transverse momentum of the produced charmonium, and the collision energy. They are therefore best studied using p+A data covering the broadest possible range of collision energy, rapidity and transverse momentum.
The assumption that effects due to soft particles produced in the collision are not important in p or d+A collision at colliders was called into question by the observation of strong suppression of the ψ(2S) relative to the J/ψ in central d+Au collisions [11], and then in p+Pb collisions [16]. Because CNM effects on the production of charm quarks and their transport through the nucleus are expected to affect both states similarly, they do not appear to be able to explain this observation. However, it can be reproduced by the co-mover break up model [23], where charmonium is dissociated by interactions with produced particles in the final state, which naturally gives a larger suppression effect on the much more weakly bound ψ(2S). The observation of flow-like behavior in p+Pb collisions at LHC (see for example [24]) and later in d+Au collisions at RHIC [25,26] suggested that a quark-gluon plasma of small size may be formed in high energy collisions of these light systems. This led to the application of transport models to p+Pb and d+Au data, which were originally developed for charmonium produc-tion in heavy ion collisions [27,28]. A plasma phase in these small collision systems gives different suppression between the charmonia states and allows a description of the data. In the case of most central midrapidity d+Au collisions at √ s N N = 200 GeV, additional suppression beyond CNM effects has been predicted of approximately 20% for the J/ψ, and 55% for the ψ(2S) [27], in good agreement with the data [9,11].
In 2014 and 2015, the Relativistic Heavy Ion Collider (RHIC) provided collisions of p+Al, p+Au and 3 He+Au for a systematic study of small systems. A comparison of flow data from p+Au, d+Au and 3 He+Au with hydrodynamic models found that the data were all consistent with hydrodynamic flow in the most central collisions [29][30][31]. An obvious question is whether increased energy density provided by the 3 He projectile in comparison to the proton produces any observable effect on charmonium modification in collisions with a Au target.
In this paper we present PHENIX measurements of inclusive J/ψ production in p+Al, p+Au and 3 He+Au collisions at √ s N N = 200 GeV. The inclusive J/ψ cross section includes feed-down from ψ(2S) and χ c states, and a smaller contribution from B-meson decays. The results are directly compared to p+p collisions at the same center of mass energy by calculating the nuclear modification factor R AB . The data are presented as a function of J/ψ p T , rapidity, and centrality, and compared to theoretical models.

II. EXPERIMENTAL SETUP
The PHENIX detector [32] comprises two central arm spectrometers at midrapidity and two muon arm spectrometers at forward and backward rapidity. The detector configuration during the data taking in 2014 and 2015 is shown in Fig. 1. The data presented here are from J/ψ → µ + µ − decays recorded with the muon arm spectrometers. The muon spectrometers have full azimuthal acceptance, covering −2.2 < η < −1.2 (south arm) and 1.2 < η < 2.4 (north arm). Each muon arm comprises a Forward Silicon Vertex Tracker (FVTX), followed by a hadron absorber and a muon spectrometer.
The FVTX [33] is a silicon detector designed to measure a precise collision vertex (also constrained by the Silicon Vertex Tracker (VTX) at midrapidity), and to provide precise tracking for charged particles entering the muon spectrometer before undergoing multiple scattering in the hadron absorber. The FVTX was not used in this inclusive J/ψ analysis, because the acceptance is reduced when requiring muon arm tracks that match tracks in the FVTX. Following the FVTX is the hadron absorber, composed of layers of copper, iron, and stainless steel, corresponding to 7.2 nuclear interaction lengths (λ I ). The absorber suppresses hadrons in front of the muon arm by a factor of approximately 1000, thus significantly reducing hadronic background for muon based measurements.
Each of the muon spectrometers is composed of a muon tracker (MuTr) embedded in a magnetic field followed by a muon identifier (MuID). Each MuTr comprises three stations of cathode strip chambers, inside a magnet with a radial field integral of B · dl = 0.72 T · m. It provides a momentum measurement for charged particles. Each MuID is composed of five layers (referred to as gap 0-4) of steel absorber (4.8 (5.4) λ I for south (north) arm) and two planes of Iarocci tubes. This enables the separation of muons and hadrons based on their penetration depth at a given reconstructed momentum. The MuID in each arm is also used to trigger events containing two or more muon tracks per event, called a dimuon trigger, and each muon track is required to have at least one hit in either gap 3 or gap 4. A more detailed discussion of the PHENIX muon arms can be found in Ref. [34,35].
The beam-beam counters (BBC) are used to determine the collision vertex position along the beam axis (z BBC ) with a resolution of roughly 2 cm in p+p collisions. Each BBC comprises two arrays of 64 quartzČerenkov detectors located at z = ±144 cm from the nominal interaction point, and has an acceptance covering the full azimuth and 3.1 < |η| < 3.9. They also provide a minimum bias (MB) trigger by requiring at least one hit in each BBC. The BBC trigger efficiency, determined from the Van der Meer scan technique [36], is 55% ±5% for inelastic p+p events and 79% ± 2% for events with midrapidity particle production [37,38]. In p+Al, p+Au, and 3 He+Au collisions, charged particle multiplicity in the BBC in the Au/Al-going direction (−3.9 < η < −3.1) is used to categorize the event centrality. The BBC trigger efficiency is 72% ± 4%, 84% ± 3%, and 88% ± 4% of inelastic p+Al, p+Au, and 3 He+Au collisions, respectively.
A Glauber model, combined with a simulation of the BBC response, is used to relate charged particle multiplicity in the BBC to parameters that characterize the collision centrality, as described in [38]. The analysis produces the average number of nucleon-nucleon collisions in each centrality category. It also produces centrality dependent BBC bias correction factors which account for the correlation between BBC charge and the presence of a hard scattering in the event, and are applied as a multiplicative correction on invariant yields. Table I shows the values of N coll and BBC bias correction factor from this analysis.

A. Data set
The data sets used in this analysis are 3 He+Au data collected in 2014, and p+p, p+Al, and p+Au data collected in 2015. All data sets were recorded at a center of mass energy √ s N N =200 GeV. The events considered here are triggered by the dimuon trigger and are required to have a vertex within ±30 cm of the center of the interaction region. The corresponding integrated luminosity  Yields of J/ψ mesons were extracted from the invariant mass spectra constructed from combinations of unlikesign tracks that are identified as muons (see Fig. 2). The mass spectra contain muon pairs from J/ψ decays, as well as significant contributions from combinations of real muons not from a J/ψ, as well as misidentified hadrons. Details about the dimuon selection to reduce the background contributions are described in [39,40].
The mass spectrum constructed from like-sign tracks was used to estimate the background due to random combinations of kinematically unrelated tracks. A modified Hagedorn function was used to represent the correlated background due to kinematically related tracks. For J/ψ signal extraction, Crystal-ball functions [41] were used to describe the J/ψ and ψ(2S) peaks, similar to the previous analysis in small collision systems [12]: where σ andm are the width and mass centroid of the Gaussian component of the line shape and α and n are parameters describing the tail.
The crystal-ball shape and tail parameters for the ψ(2S) were fixed with respect to the J/ψ parameters, using the PDG database value [42] for the energy difference and a width broadening factor taken from simulations. In cases where the statistical precision of the data led to poor definition of the J/ψ signal shape, the mass and width of the J/ψ peak were fixed and a systematic uncertainty was assigned to the yield based on tests made with higher statistics cases. The statistical uncertainties related to the extraction of the J/ψ yields were determined from a covariance matrix in the fitting procedure.

C. Background estimation
The random combinatorial background in the unlikesign mass spectrum was approximated by combining all like-sign tracks from the same events. There is a small correlated contribution to the like-sign pairs from jets and open bottom; however, compared to the other background sources, this is small.
The correlated background comprises unlike-sign muon pairs from charm, bottom, jets, and Drell-Yan. Because the correlated background cannot be estimated independently from the data, it must be fitted to the mass spectrum when the J/ψ yield is extracted. Fitting the correlated background effectively compensates for the small correlated component included in the like-sign estimation of the combinatorial background.
We describe the correlated background using a modified Hagedorn function [39,43,44]: where m µµ is the reconstructed J/ψ mass, p 0 is a nor-malization parameter, p 4 is the high mass tail parameter, and p 1 , p 2 and p 3 are additional fit parameters. It was found during the analysis that when fitting mass spectra with poor statistical precision, the shape of the correlated background was not well defined. This led to a contribution of less than 10% to the point-to-point uncertainty in the J/ψ yields. Therefore, the shape of the correlated background as a function of p T (determined by p 1 , p 2 and p 3 ) was constrained using simulation results based on a detailed study of dimuon mass spectra [9,39,45,46]. A systematic uncertainty on the J/ψ yield was assigned for this procedure by refitting the data with various combinations of correlated background parameters left free.

Acceptance and Reconstruction Efficiency
The study of acceptance and reconstruction efficiency of dimuons from J/ψ decays has been performed using a geant4-based full detector simulation [47]. In this simulation, the MuTr and MuID detector efficiencies are set to values determined from the data. An emulator of the dimuon trigger response is included in the simulation to account for the trigger efficiency. As these efficiencies depend on the instantaneous luminosity being sampled, each data set is divided into three groups with different beam interaction rates, and corrected yields with separate corrections are compared. A systematic uncertainty is assigned to the extracted J/ψ cross sections to reflect the differences, see Sec. III G for details.
The pythia8 event generator package [48] is utilized to generate full J/ψ events used for the full Geant4 detector simulation. To take into account effects from background hits, the simulated hits of pythia8 events are embedded into real data events, separated into centrality classes of the collision system. The track reconstruction is then run on the data with embedded simulated hits to examine the effects of the underlying event on the reconstruction efficiency. Figure 3 shows the acceptance and reconstruction efficiency for the J/ψ as a function of p T in p+p collisions. The difference between the two muon arms is mainly from different inefficient detector areas.

E. Invariant yield and nuclear modification factor
The invariant yield of dimuons from J/ψ decays in a given rapidity and centrality bin for the integrated p T range is where B ll is the branching ratio of J/ψ to dimuons, ∆y is the width of the rapidity bin, N J/ψ is the number of J/ψ obtained from the fit procedure, c BBC is the BBC bias correction factor described in Table I, N evt is the number of sampled MB events in the given centrality bin, ε Ae is the J/ψ acceptance and reconstruction efficiency, and ε trig is the dimuon trigger efficiency. The invariant yield in a y, p T , and centrality bin is where ∆p T is the width of the p T bin, and in this case N evt is the number of events in the centrality bin. Based on the invariant yields calculated with Eq. 4, the J/ψ nuclear modification factor R AB for a given y, p T , and centrality bin is formed to quantify nuclear effects in p+Al, p+Au, and 3 He+Au collisions. The R AB is defined as where d 2 N AB /dydp T is the J/ψ invariant yield for a certain centrality bin of A+B collisions, d 2 N pp /dydp T is the corresponding J/ψ invariant yield for p+p collisions, and N coll is the mean number of binary collisions for that centrality bin in A+B collisions.
The p 2 T values for various centrality bins in all collision systems have been calculated over the full measured p T range (0 < p T < 7 GeV/c). We do not extrapolate the p T distribution beyond 7 GeV/c. A previous study [10] determined that extrapolating to infinite p T increased the p 2 T values by 3%. The value of p 2 T is calculated numerically using the following formula: where p T,i is the center of the i-th p T bin, and w i is the weight factor proportional to the J/ψ invariant yield in the p T bin: where dp T,i is the width of the bin.

G. Systematic Uncertainties
In the measurements we present in the next section, Type A uncertainties are random point to point uncertainties, and are dominated by the statistical precision of the data. Type B systematic uncertainties are correlated point to point uncertainties. Type C global uncertainties are fractional uncertainties that apply to all measurements uniformly.

Signal extraction
As discussed in Sec. III C, the modified Hagedorn function in Eq. 2 was used to describe the correlated background. Initial parameters were estimated based on the previous measurement of dimuon mass spectra [39,45], and two parameters, p 0 and p 4 , were left free to describe dimuon mass distributions in the data more properly. For the systematic uncertainty study, additional parameters, p 1 , p 2 , and p 3 , in the modified Hagedorn function were also freed in the fit procedure. We observe 1.4%-2.8% variations of J/ψ counts depending on rapidity, p T , and centrality.
To describe the combinatorial background shape, the modified Hagedorn function in Eq. 2, used for the correlated background component, was also used to fit likesign dimuon mass distributions. The effect of statistical fluctuations in the like-sign dimuon mass distributions was studied by varying the shape based on the statistical uncertainties of the fit parameters. We observe 1.0%-4.4% variations of J/ψ counts depending on rapidity, p T , and centrality.
The uncertainty related to fixing the J/ψ mass centroid and width was evaluated by directly comparing the difference in yields with the parameters free versus fixed, which ranges from 1.1%-2.9% uncertainty. Table II lists all Type B uncertainties arising from the J/ψ signal extraction.

Acceptance and efficiency correction
The acceptance and reconstruction efficiency correction and trigger efficiency correction are obtained from simulation, so discrepancies between the data and calculations can be a source of systematic uncertainty. The discrepancies can be due to a variation in the detector performance during the data taking period and/or inaccuracy of detector geometry and dead channel maps in the simulation. To quantify these effects, we divide each data set into three groups of different detector efficiency, based on the beam instantaneous luminosity and calculated invariant yields with separate correction factors. In this comparison we observe 1.5%-5.0% variations, depending on rapidity and data set, and assign this variation as a systematic uncertainty. In addition, we compare the azimuthal angle φ distribution of tracks in the MuTr between the data and simulation, and assign a 2.5%-6.0% systematic uncertainty depending on rapidity and data set.
In the simulation procedure, pythia8 was used to generate J/ψ events, and initial J/ψ rapidity and p T shapes in pythia8 are tuned to match the measurements in p+p and d+Au collisions [10,40,49]. These two different assumptions of the distributions are used as bounds to estimate the sensitivity of this analysis to the shapes of these distributions in p+Al, p+Au, and 3 He+Au collisions, which are not known a priori. The variation of acceptance and reconstruction efficiency between two sets of rapidity and p T distributions is less than 2%, so we assigned a 2% conservative systematic uncertainty.
The uncertainty in the dimuon acceptance caused by lack of knowledge of the J/ψ polarization was studied as described in [40]. Because there is no precise measurement of J/ψ polarization, a maximum polarization value (±1 in the helicity frame) was considered to study the systematic uncertainty. The variation of dimuon acceptance becomes larger as J/ψ p T decreases, and 9%-20% systematic uncertainties are assigned depending on p T . We assumed that the J/ψ polarization is not significantly modified in p+Al, p+Au, and 3 He+Au collisions, and this uncertainty is canceled in the R AB calculation. This assumption was also made in a similar PHENIX analysis for J/ψ nuclear modification in d+Au collisions [10].
To evaluate a systematic uncertainty on the dimuon trigger efficiency, the single muon trigger efficiency in the MB triggered data obtained with a large number of muon samples was compared with the emulated single muon trigger efficiency determined from simulation. This difference was propagated to the uncertainty in the dimuon trigger efficiency based on a previous study [39], and a 1.0%-4.8% systematic uncertainty was assigned. The Type B systematic uncertainties related to acceptance and efficiency correction are shown in Table III.

Multiple interaction
Due to the high instantaneous beam luminosity, particularly in p+p and p+Al runs, it is possible to have multiple inelastic collisions from a single beam crossing, which can affect the invariant yield calculation. To investigate this effect, the variation among invariant yields in three groups of different instantaneous luminosity for each data set was studied, revealing a yield variation smaller than 5%. However, the instantaneous luminosity dependence of the acceptance and efficiency correction is already included as a systematic uncertainty, and so no additional systematic uncertainty is assigned.

N coll and BBC efficiency
The systematic uncertainties on the BBC efficiency and the determination of N coll in p+Al, p+Au, and 3 He+Au collisions described in Table I are evaluated by following the procedure developed in the previous PHENIX analyses of d+Au data [38]. The systematic uncertainty on the BBC efficiency in p+p collisions obtained in [37] is 10.1%.

IV. RESULTS
In this section, we present invariant yield, nuclear modification factor, and p 2 T results at forward and backward rapidity. There have been significant changes to the muon arm configuration and to the simulation framework because the d+Au data set was recorded. Figure 4 shows the J/ψ invariant yield as a function of p T in p+p collisions at √ s = 200 GeV at forward and backward rapidity, where bars (boxes) represent point-to-point uncorrelated (correlated) uncertainties. The global systematic uncertainty is 10.1%. The ratio of invariant yields between the forward and backward rapidity regions is presented in the bottom panel, where the systematic uncertainty due to the J/ψ polarization cancels in the ratio. The invariant yields at forward and backward rapidity are consistent within the systematic uncertainties, confirming that the detector efficiency is well understood.
Plots and tables of invariant yield are presented for the other collision systems in the Appendix. We focus here on the nuclear modification factors. Figure 5 shows the rapidity dependence of the nuclear modification factor for 0%-100% centrality in p+Al, p+Au, and 3 He+Au collisions. The rapidity dependence of the nuclear modification for different centrality classes is shown for p+Al in Fig. 6, for p+Au in Fig. 7, and for 3 He+Au in Fig. 8. Figures 9 and 10 show the nuclear modification factor as a function of p T for 0%-100% p+Al, p+Au, and 3 He+Au collisions at backward and forward rapidity. The p T dependence in different centrality classes is presented for p+Al in Fig. 11, for p+Au in Fig. 12 and 13, and for 3 He+Au in Fig. 14. The modification as a function of p T in 0%-20% central collisions is compared between p+Al and p+Au in Fig. 15. Similar comparisons where the target is identical, but the projectile is different are shown for 0%-20% central collisions comparing d+Au and p+Au in Fig. 16 and comparing 3 He+Au and p+Au in Fig. 17.
The p T integrated nuclear modification factor for p+Al, p+Au and 3 He+Au as a function of N coll is shown at both forward and backward rapidity in Figs. 18 and 19. A comparison between p+Al, p+Au and 3 He+Au modifications when plotted as a function of the average nuclear thickness sampled by the charmonium production is presented in Fig. 20. Figure 21 shows the mean p T squared values for the three systems p+Al, p+Au and 3 He+Au as a function of N coll for p T < 7 GeV/c at forward and backward rapidity.

A. Rapidity dependence
The rapidity dependence of the modification for 0%-100% centrality, seen in Fig. 5, shows only weak modification for p+Al collisions. For both p+Au and 3 He+Au significant suppression is seen at forward rapidity, with smaller suppression at backward rapidity. The modifications for p+Au and 3 He+Au are very similar.
The rapidity dependence in three centrality bins for p+Al collisions, seen in Fig. 6, shows only weak modification in all centrality bins, both at forward and backward rapidity.
The p+Au data presented here contain finer centrality binning for central collisions than was previously available from d+Au. The rapidity dependence in six centrality bins for p+Au collisions, seen in Fig. 7, shows a factor of more than two suppression at the most forward rapid-       ity in the 0%-5% centrality bin, and a marked increase in suppression with increasing rapidity in the forward direction. At backward rapidity, the modifications in all centrality bins show little centrality dependence, all being somewhat suppressed.
The rapidity dependence of the 0%-100% centrality data is compared in Fig. 5 with model calculations from R. Vogt [50,51] and Shao et al. [52][53][54][55] showing the effect of nPDF modifications using the Eskola-Paakkinen-Paukkunen-Salgado (EPPS16) [1] next-to-leading order (NLO) and/or nuclear coordinated theoretical and experimental tests of quantum chromodynamics (nCTEQ15) NLO parameterizations [2]. The Vogt EPPS16 NLO shadowing calculations in general follow the methods de-scribed in [50], while the J/ψ mass and scale parameters are discussed in [51]. The Shao, et al. model calculations for p+Au collisions are based on a Bayesian reweighting method which uses tighter J/ψ constraints from p+Pb data at the LHC [52]. The dominant uncertainty in the reweighting method is the energy scale dependence, µ 0 , where µ 2 0 = M 2 + p T 2 for the J/ψ mass and transverse momentum. The reweighting however is not applied for lighter 3 He and Al nuclei, with the predictions for these nuclei based on the original method described in [53][54][55]. For these predictions, the previous PHENIX J/ψ measurement in p+p collisions [40] is used as a baseline. The calculations were performed at three different energy scales (µ 0 , 0.5 µ 0 , and 2 µ 0 ) and provide two different The theory bands are discussed in the text. Note that the theory bands compared with the 0%-5% and 5%-10% centrality data are for 0%-10%. FIG. 13. Nuclear modification factor of inclusive J/ψ as a function of pT at 1.2 < y < 2.2 in six centrality bins for p+Au collisions. The theory bands are discussed in the text. Note that the theory bands compared with the 0%-5% and 5%-10% centrality data are for 0%-10%.  confidence levels (68% and 90% CL). The uncertainty band shown is for the 68% CL, and we have taken the envelope of the uncertainty bands from the calculations at the three energy scales.
The calculations describe the data very well at forward rapidity for all three systems, and for p+Al at backward rapidity. For p+Au and 3 He+Au at backward rapidity the calculated modifications are too large by roughly 40%. However, the calculations do not contain effects of nuclear absorption, which is expected to be important at backward rapidity in PHENIX at √ s N N = 200 GeV [4], where the nuclear crossing time is comparable with the charmonium formation time. That is not expected to be the case at forward rapidity in PHENIX at √ s N N = 200 GeV, or at the rapidities of interest at LHC energies. Because nuclear absorption is not included in the model calculations, they should be expected to overpredict the modification in p+Au and 3 He+Au at backward rapidity.
An estimate of the effect of nuclear absorption at backward rapidity can be obtained from a model [5] fitted to absorption cross sections derived from shadowing corrected data measured at a broad range of beam energies [4]. The model assumes that the cc pair size grows linearly with time until it reaches the size of a fully formed charmonium meson. Then the absorption cross section depends on the proper time before the pair escapes the target. The effect of the modification due to nuclear absorption at backward rapidity from this model is added to Fig. 5, by folding it into the shadowing calculation. The results indicate that the measured modifications are reasonably consistent with shadowing plus nuclear absorption.

B. pT dependence
The p T dependence for 0%-100% centrality, seen at backward rapidity in Fig. 9 and at forward rapidity in Fig. 10, shows little modification for p+Al but shows strong, and similar, p T dependence for p+Au and 3 He+Au. These data are also compared with the calculations of Shao et al. [52]. As for the rapidity dependence, the calculations describe the forward rapidity data well for all three collision systems and for the backward rapidity p+Al. But the backward rapidity modification for p+Au and 3 He+Au is overpredicted. Significant nuclear absorption is expected at backward rapidity and low p T , and calculations that do not include it should overpredict the modification there.
The p+Au modifications vs p T , seen at forward rapidity in Fig. 13 for all centrality bins, shows very strong dependence on centrality. The modification falls to 0.35 at low p T for the 5% most central collisions. At backward rapidity, as shown in Fig. 12 the suppression is considerably weaker at low p T for the most central collisions, but it changes more slowly with centrality. The result is that for collision centralities above 20% the behavior of the modification versus p T becomes rather similar at forward and backward rapidity.
The theory curves shown in Figs. 12 and 13 are adapted transport models provided by X. Du and R. Rapp, based on the original transport model by Zhao & Rapp for A+A collisions [56]. The theory was extended for d+A collisions [27] and most recently for p+A collisions [57]. The transport model includes a fireball generated by a Monte-Carlo Glauber model [58] in addition to shadowing from Eskola-Paukkunen-Salgado (EPS09) [59] NLO and an ab-sorption cross section at backward rapidity. The J/ψ production cross section is described in [57], and charged particle multiplicity [60], hadronic dissociation rates [27], and open charm production cross sections [57] are also considered. The calculations reproduce the data at high p T , but generally underpredict the suppression at low p T at forward rapidity. Because the modification of J/ψ production in the transport model is not very strong at forward rapidity, the suppression there is dominated by the EPS09 shadowing contribution.
The comparison in Fig. 15 of 0%-20% p+Al with p+Au modifications contrasts the weak modification in central p+Al collisions with the strong modification, particularly at forward rapidity, in central p+Au collisions.
In a previous PHENIX measurement of charged particle multiplicity [60], it was found that twice as many par- ticles are produced in 0%-20% central 3 He+Au collisions than in 0%-20% central p+Au collisions. To look for evidence of an effect from this, Fig. 17 shows a direct comparison between the modifications for p+Au and 3 He+Au in the 0%-20% centrality bin. The ratio of 3 He+Au to p+Au is included in Fig. 17. All systematic uncertainties from the p+Au and 3 He+Au systems are included except the initial shape uncertainty, which cancels upon taking the ratio. All systematic uncertainties stemming from the p+p system cancel. A mean value has been fitted to the ratios, and the result is shown on the plot together with the fit uncertainty and the uncertainty from the systematic errors. The systematic uncertainty was determined by remaking the fit with all points moved to the upper or lower limits of their systematic uncertainty. The ratio at forward rapidity is R3 HeAu /R pAu = 0.96 ± 0.03(stat) ± 0.05(syst), which is consistent with unity. At backward rapidity it is The results are consistent with J/ψ production being reduced for the 3 He projectile, with the backward rapidity ratio having a probability of 90% of being less than one.

C.
N coll dependence The p T integrated modifications as a function of N coll in each centrality bin are shown in Fig. 18. No scaling with N coll is expected between p+Au and 3 He+Au, because 3 He+Au will have roughly three times as many collisions as p+Au in the same centrality class. The N coll dependence of the p+Au modification is shown again in Fig. 19, where it is compared with the p T integrated modification predicted by Du and Rapp. The theory calculation shows both the CNM baseline and the effect of the transport model. Again, it is seen that the suppression is underpredicted. At backward rapidity some nuclear absorption is expected. At forward rapidity, it appears that the CNM effects are not strong enough to explain the data. However, the model predicts a suppression beyond CNM effects at backward rapidity for central collisions of approximately 10%.
Modifications that are due to CNM effects (including nuclear absorption) would be expected to depend on the thickness of the target nucleus at the impact parameter of the nucleon that was involved in the hard process. The nuclear thickness can be written where ρ A (z, r T ) is the density distribution of nucleons in nucleus A taken from the Woods-Saxon distribution used in the Glauber model discussed in section II. The parameter z is the location in the nucleus along the beam direction, and r T is the transverse distance from the center of the nucleus. T A (r T ) is the average number of nucleons per unit area at the projectile nucleon impact parameter r T . To get the average value of T A sampled for charmonium production within a given centrality bin, the values of T A (r T ) are weighted by the distribution of r T values within the centrality bin, to reflect the number of projectile nucleons having one or more inelastic collisions at that r T , and additionally by the probability of a hard process at that r T -which is proportional to T A (r T ). Figure 20 shows the p+Al, p+Au and 3 He+Au modifications plotted versus T A , in each centrality bin. The modifications seem to fall on a common curve within uncertainties, as would be expected if they were primarily due to CNM effects.
The p 2 T values versus N coll , shown in Fig. 21, fall on a common curve for all three systems. The N coll dependence is mild, with p 2 T increasing from 3.3 in p+p collisions to approximately 4 in p+Au and 3 He+Au collisions. The p 2 T is very similar between forward and backward rapidity, as was also observed in d+Au collisions [10].

VI. SUMMARY AND CONCLUSIONS
We have presented invariant yields for inclusive J/ψ production in p+p, p+Al, p+Au and 3 He+Au collisions at √ s N N = 200 GeV, and the corresponding nuclear modifications for p+Al, p+Au and 3 He+Au. The new p+Au results are found to agree within uncertainties with the previous PHENIX d+Au results [9]. The p+Al modifications are found to be much weaker at all centralities than those in p+Au. The 0%-100% centrality data for p+Al are found to be well described in rapidity and p T by calculations containing only shadowing effects from the EPPS16 NLO and nCTEQ15 NLO parameterizations, aside from slightly underpredicting the modification at 4-6 GeV/c at forward rapidity.
The 0%-100% centrality p+Au and 3 He+Au data are also compared with calculations based on the EPPS16 NLO and nCTEQ15 NLO shadowing parameterizations. At forward rapidity, the calculations describe the p+Au and 3 He+Au modifications well in both rapidity and p T , again with the exception of slightly underpredicting the modification at 4-6 GeV/c at forward rapidity. At backward rapidity, the calculations overpredict the modifications. We found that adding the predicted nuclear absorption modification taken from previous work to the backward rapidity p T integrated data reduced the modifications to values consistent with the data.
The ratio of the 3 He+Au and p+Au modifications for the 0%-20% centrality bin at forward rapidity is R3 HeAu /R pAu = 0.96 ± 0.03(stat) ± 0.05(syst), which is smaller but consistent with unity. At backward rapidity it is The results are consistent with a reduction in the modification for the heavier projectile case. Given the systematic uncertainty, the backward rapidity ratio has a 90% probability of being less than 1.0.
For p+Au at forward rapidity, the nuclear modification vs p T shows very strong centrality dependence, dropping to approximately 0.35 at low p T in the most central 5% of collisions. At backward rapidity the suppression is weaker for central collisions, but it changes more slowly. Comparison with theory calculations that include EPS09 shadowing and a final state transport model are able to reproduce the general shape of the p T dependence at each centrality, but greatly underpredict the suppression at low p T for central collisions.
The p T integrated modification for p+Au drops steeply with centrality at forward rapidity, reaching approximately 0.5 for the 5% most central collisions. The modification at backward rapidity is found to have weak centrality dependence. Because nuclear absorption is evidently important at backward rapidity, the weak centrality dependence there is likely due to a trade-off between anti-shadowing and nuclear absorption. It was found that plotting the modification vs T A for each centrality bin caused them to fall on a common line for all three systems, as would be expected if CNM effects dominate.

APPENDIX
The invariant yields for all data sets are presented in this appendix. Figure 22 shows inclusive J/ψ invariant yield as a function of rapidity in MB p+p, p+Al, p+Au, and 3 He+Au collisions, and the invariant yields in p+Al, p+Au, and 3 He+Au collisions are scaled with N coll to compare with the invariant yield in p+p collisions. In this and following figures showing results of invariant yield measurement, bars (boxes) around data points represent point-to-point uncorrelated (correlated) uncertainties. Figures 23, 24, and 25 show inclusive J/ψ invariant yield as a function of rapidity in different centrality of p+Al, p+Au, and 3 He+Au collisions, respectively. Invariant yields in p+Al, p+Au, and 3 He+Au collisions are scaled with N coll , and the p+p result is also presented in each panel. Figures 26, 27, and 28 show inclusive J/ψ invariant yield as a function of p T in different centrality of p+Al, p+Au, and 3 He+Au collisions, respectively.  There is also a global uncertainty of 13.6%, 12.2%, and 12.3% corresponding to 0%-20%, 20%-40% and 40%-72% centrality. There is also a global uncertainty of 11.9%, 11.8%, 12.2%, 12.1%, 12.2% and 14.0% corresponding to 0%-5%, 5%-10%, 10%-20%, 20%-40%, 40%-60%, and 60%-84% centrality. J/ψ invariant yield as a function of pT in various centrality bins of 3 He+Au collisions, and the yields in each centrality bin are scaled for better visibility. Bars (boxes) around data points represents point-to-point uncorrelated (correlated) uncertainties. There is also a global uncertainty of 10.7% for 40%-88% centrality and 10.2% for all remaining centralities.