Nuclear-modification factor of charged hadrons at forward and backward rapidity in $p$$+$Al and $p$$+$Au collisions at $\sqrt{s_{_{NN}}}=200$ GeV

The PHENIX experiment has studied nuclear effects in $p$$+$Al and $p$$+$Au collisions at $\sqrt{s_{_{NN}}}=200$ GeV on charged hadron production at forward rapidity ($1.4<\eta<2.4$, $p$-going direction) and backward rapidity ($-2.2<\eta<-1.2$, $A$-going direction). Such effects are quantified by measuring nuclear modification factors as a function of transverse momentum and pseudorapidity in various collision multiplicity selections. In central $p$$+$Al and $p$$+$Au collisions, a suppression (enhancement) is observed at forward (backward) rapidity compared to the binary scaled yields in $p$+$p$ collisions. The magnitude of enhancement at backward rapidity is larger in $p$$+$Au collisions than in $p$$+$Al collisions, which have a smaller number of participating nucleons. However, the results at forward rapidity show a similar suppression within uncertainties. The results in the integrated centrality are compared with calculations using nuclear parton distribution functions, which show a reasonable agreement at the forward rapidity but fail to describe the backward rapidity enhancement.

is observed at forward (backward) rapidity compared to the binary scaled yields in p+p collisions. The magnitude of enhancement at backward rapidity is larger in p+Au collisions than in p+Al collisions, which have a smaller number of participating nucleons. However, the results at forward rapidity show a similar suppression within uncertainties. The results in the integrated centrality are compared with calculations using nuclear parton distribution functions, which show a reasonable agreement at the forward rapidity but fail to describe the backward rapidity enhancement.

I. INTRODUCTION
Measurements of particle production in heavy-ion collisions enable the study of properties of a hot and dense nuclear medium called the quark-gluon plasma (QGP) [1][2][3][4]. An initial striking observation at the Relativistic Heavy Ion Collider (RHIC) was that production of high transverse momentum (p T ) hadrons in Au+Au collisions is strongly suppressed compared to that in p+p collisions scaled by the number of binary collisions. This suppression indicates that partons experience substantial energy loss as they traverse the QGP, a phenomenon called jetquenching [5]. A control experiment involving a deuteron projectile on a heavy-ion target, d+Au, was carried out to test whether the feature of strong energy loss is still present in a collision system of much smaller size. The results in d+Au collisions at midrapidity presented in Ref. [6] showed no suppression at high p T , initially leading to the conclusion that QGP itself-and associated jet quenching-were unique to collisions of larger heavy ions. In the ten years because these initial measurements, indications of QGP formation in smaller collision systems including d+Au have been found, though without evidence of jet quenching phenomena [7].
Although there were no indications of strong suppression of high p T particles in d+Au collisions, detailed measurements do indicate other particle-production modifications relative to p+p collisions [8][9][10][11][12]. At midrapidity, a centrality-dependent enhancement of charged hadron production was observed at intermediate p T (2 < p T < 5 GeV/c) [11] in d+Au collisions at √ s N N = 200 GeV.
These nuclear effects may be due to initial-and/or finalstate multiple scatterings of incoming and outgoing partons [13,14]. Processes such as radial flow [15] and recombination [16] developed for heavy-ion collisions were also investigated to explain a stronger enhancement of p andp over π ± and K ± [11]. Recent results of collectivity amongst identified particles in small collision systems at RHIC and the Large Hadron Collider [7] have been also explained within the hydrodynamic evolution model [17,18]. The study of particle production at forward and backward rapidity can provide additional information on nuclear effects such as initial-state energy loss [19] and modification of nuclear parton distribution functions (nPDF) [20][21][22][23][24]. Of particular interest are gluons at * PHENIX Spokesperson: akiba@rcf.rhic.bnl.gov † Deceased small Bjorken x Bj (fraction of the proton's longitudinal momentum carried by the parton), where the dramatic increase of gluon density leads to expectation of saturation. This is often described within the color glass condensate (CGC) framework [25]. A strong centrality dependent suppression of single and dihadron production has been observed at forward rapidity in d+Au collisions at √ s N N = 200 GeV [8][9][10]. A CGC calculation provides a good description of the experimental data [26,27]. Also, a perturbative quantum chromodynamics (pQCD) calculation considering coherent multiple scattering with small-x Bj gluons reproduces the suppression of particle production at forward rapidity [19,28]. Another very different explanation for the suppression at forward rapidity is that color fluctuation effects modify the size of the high-x Bj partons in the proton [29,30]. Accessible quark and gluon x Bj ranges depend on the pseudorapidity (η) and transverse momentum of final state hadrons or jets. Therefore, measurements over a wide kinematic range are quite useful to further understand nuclear effects in small collision systems. PHENIX experiment has two muon spectrometers that provide wide coverage at forward (x Bj ≈ 0.02, shadowing region) and backward rapidity (x Bj ≈ 0.1, antishadowing region). In the previous study of nuclear effects on charged hadron production in d+Au collisions at √ s N N = 200 GeV [31], a significant suppression was observed at forward rapidity in high multiplicity collisions compared to that in low multiplicity collisions, whereas a moderate enhancement is seen at backward rapidity. Although the direction of modification is consistent with the expectation from nPDF modification, no specific model comparison was presented.
High statistics data samples of p+p, p+Al, and p+Au collisions at √ s N N = 200 GeV were collected in 2015 by PHENIX. These data samples combined with the availability of a new forward silicon vertex tracking detectors, which enable the selection of particle tracks coming from the collision point, significantly improved p T and η resolutions. The charged hadron analysis with these data sets can extend the previous study in d+Au collisions [31], and a comparison between p+Al and p+Au of very different size of nuclei can provide new information on nuclear effects on charged hadron production in p+A collisions.
In this paper, we present nuclear modification factors of charged hadron production at forward and backward rapidity in p+Al and p+Au collisions at √ s N N = 200 GeV of various multiplicities. Section II describes the experimental setup and the data sets used in this analysis. Section III details the analysis methods.
Section IV discusses systematic uncertainties. Section V presents results and discussion. Section VI gives the summary and conclusions. 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 2015 is shown in Fig. 1. 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 muon spectrometer is composed of a muon tracker (MuTr) embedded in a magnetic field followed by a muon identifier (MuID).

II. EXPERIMENTAL SETUP
The FVTX is a silicon detector with four stations in each arm. Each station comprises 96 sensors along the φ direction. Each silicon sensor is finely segmented along the radial direction, with a strip pitch of 75 µm. The primary purpose of the FVTX is to measure a precise collision vertex also constrained by the silicon vertex tracker (VTX) at midrapidity. The FVTX was also designed to measure precise momentum vector information of charged particles entering the muon spectrometer before suffering large multiple scattering in the hadron absorber. More technical details on the FVTX are available in Ref. [33]. Following the FVTX is the hadron absorber, composed of layers of copper, iron, and stainless steel, corresponding to 7.2 nuclear interaction lengths (λ I ). Hadrons entering the absorber are suppressed by a factor of approximately 1000, thus significantly reducing hadronic background for muon-based measurements.
The MuTr has two arms each consisting of three stations of cathode strip chambers, which are inside a magnet with a radial field integral of B · dl = 0.72 T ·m. The MuTr provides a momentum measurement for charged particles. The 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 MuTr and MuID are also used to trigger events containing at least one muon or hadron candidate. The MuID trigger is designed to enrich events with muons by requiring at least one hit in either gap 3 or 4. Hadrons that stop only after partially penetrating the MuID can be enhanced by requiring no hit in gap 4. The MuTr trigger is used to sample high momentum tracks by requiring a track sagitta less than three MuTr cathode strips wide at the middle station of the MuTr. A more detailed discussion of the PHENIX muon arms can be found in Ref. [34,35].
The beam-beam counters (BBC) [36] comprise two arrays of 64 quartzČerenkov detectors located at z = ±144 cm from the nominal interaction point. Each BBC has an acceptance covering the full azimuth and 3.1 < |η| < 3.9. The BBCs 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. 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 [37], is 55% for inelastic p+p events and 79% for events with midrapidity particle production [38]. In p+Al and p+Au collisions, charged particle multiplicity in BBC in the Aland Au-going direction (−4.9 < η < −3.1) is used to categorize the event centrality. The BBC trigger is for 72% (84%) of inelastic p+Al (p+Au) collisions. Centrality dependent bias factors to account for the efficiency for MB triggered events and hard scattering events have been obtained based on the method developed in Ref. [39].

B. Hadron selection
The majority of hadrons emitted from the collision are stopped inside the hadron absorber. Hadrons which pass through the hadron absorber enter the MuTr and can still be stopped in the middle of the MuID by producing hadronic showers in the additional steel absorber planes. Low momentum muons can also be stopped due to ionization energy loss, but the momentum distribution measured in the MuTr is very different for these muons and hadrons which are stopped in the MuID. Figure 2 shows the longitudinal momentum (p z ) distributions of reconstructed tracks at the north arm MuID gaps 2 and 3 from a full geant4 detector simulation of charged hadrons (see Sec. III D). Muon tracks from light hadron decays show a narrow p z distribution in 2.5 < p z < 3.0 GeV/c, whereas tracks from hadrons show a much broader distribution. Therefore, tracks from hadrons can be enriched with a proper p z cut (3.5 GeV/c for gap 2 and 4 GeV/c for gap 3). The inset plots show the hadron fraction as a function of p T with the p z cuts. The hadron purity is > 98% (> 90%) at MuID gap 2 (gap 3) for p T > 1.5 GeV/c. The contamination of muons in the combined sample for both MuID gap 2 and gap 3 is less than 5% based on this simulation study.
One benefit from the FVTX is that the initial momentum vector of hadrons can be measured precisely before they undergo significant multiple scattering inside the absorber. In particular, the FVTX has very fine segmentation in the radial direction which can improve the p T and η resolution of measured tracks, both of which are important for this analysis. Figure 3 shows the ∆η distribution between reconstructed tracks (η Reco ) and true tracks (η Gen ) as a function of p T for hadron candidates from the geant4 simulation. In the case where momen-tum information from only the MuTr is used, shown in Fig. 3 (a), the smearing in η is quite large. This is significantly improved by requiring association with FVTX tracks, shown in Fig. 3

C. Trigger efficiency
One consideration with the FVTX association requirement is the possibility of multiple FVTX tracks within the search window of a projected MuTr track, due to the higher FVTX track multiplicity and the smeared momentum information from the MuTr as shown in Fig. 3 (a). In this case, a MuTr track can be associated with a wrong FVTX track. This is referred to as a mis-association. Such mis-associations result in further smearing of the reconstructed p T and η. The FVTX-MuTr association efficiency depends on the event multiplicity. The probability of mis-association can be evaluated with a data driven method developed in [40,41] by associating a MuTr track with FVTX tracks from another event of similar FVTX track multiplicity. The same method has been used in this analysis, and the estimated fraction of misassociations in the p+p data is ∼ 1.5% at p T ≈ 1.5 GeV/c and decreases down to ∼ 0.5% at p T ≈ 5 GeV/c. In the 0%-5% highest multiplicity p+Au collisions, the estimated fraction of mis-associations in the south arm (Augoing direction) is ∼ 3% at p T ≈ 1.5 GeV/c and ∼ 1% at p T ≈ 5 GeV/c, which is a factor of two higher than the estimate for p+p collisions. The mis-association fraction is also checked with hadron simulation events embedded into real data events, and is consistent with the data driven values. The embedding simulation described in Sec. III D is used to take into account the multiplicity dependent FVTX-MuTr association efficiency.  In addition to the requirements on p z and FVTX-MuTr association, track quality cuts are applied. MuTr tracks are required to have at least 11 hits out of a maximum of 16 hits, and a 3σ MuTr track fit quality cut is applied. For association between MuTr and MuID tracks, three standard deviation cuts are applied to the angle and distance between MuTr and MuID tracks projected to the MuID gap 0. The associated FVTX track is required to have hits in at least three of the four stations, and an additional 3σ fit quality cut is applied. Momentum-dependent cuts are applied to the angle difference in the radial and azimuthal directions between FVTX and MuTr tracks projected to the middle of the absorber (z = 70 cm). These selections help reject tracks Acceptance and reconstruction efficiency for different charged hadrons as a function of pT in the south arm evaluated in p+p collisions. from decay muons, secondary hadrons, and FVTX-MuTr mis-associations.
The trigger efficiency is evaluated using hadron candidates from MB triggered events by measuring the fraction of hadron candidates satisfying the trigger requirements. Figure 4 shows the trigger efficiency for hadrons as a function of p T at MuID gap 2 and gap 3 of the south arm in the p+p data. The trigger efficiency for hadrons at MuID gap 3 is higher than that for hadrons at MuID gap 2. The efficiency at the north arm in the p+p data is similar. Due to the larger statistical fluctuations at p T > 5 GeV/c, a fit function is used to obtain the p Tdependent trigger efficiency correction factors. The trig-   ger efficiency is separately evaluated for each muon arm as well as each centrality bin of p+Al and p+Au collisions to account for possible multiplicity effects and detector performance variation during the data taking period. The relative variation of the trigger efficiency over the data taking period is less than 10%. Because this variation of the trigger efficiency is accounted for by the detector performance variation described in Sec. III D, no additional systematic uncertainty is assigned.

D. Acceptance and reconstruction efficiency
Calculation of the absolute acceptance and efficiency for hadrons requires a detailed simulation of the hadronic interactions in the thick absorber material. There are significant uncertainties as observed from various geant4 implementations of such interactions. However, the response of hadrons inside the absorber is independent of collision systems, and hence this uncertainty will cancel out when comparing hadron yields between two collision systems. Therefore, nuclear effects on hadron production can be studied by taking into account only the additional multiplicity-dependent efficiency corrections. To obtain the multiplicity-dependent efficiency corrections, a full geant4 detector simulation was developed as follows: 1. Generate a mixture of hadrons (π ± , K ± , K 0 S , K 0 L , p, andp) based on initial p T and η distributions studied in [12,42]. Based on measurements of identified charged hadrons at midrapidity [11,43,44], an extrapolation to forward and backward rapidity is done by multiplying the ratio of p T spectra between mid and forward/backward rapidity from event generators [45,46]. These simulated hadrons originate from a z distribution which matches the measured z BBC data.
2. Run a full geant4 simulation for the detector response of hadrons.
3. Reconstruct simulated detector hits embedded on top of background hits from real data for each centrality bin in each collision system. Apply the datadriven detector dead channel maps to account for variations in detector performance. Figure 5 shows an example of acceptance and efficiency result as a function of p T for different species of hadrons at MuID gaps 2 and 3 of the south arm in p+p collisions. The acceptance and efficiency for π ± and K − is comparable, and K + has the highest acceptance and efficiency due to its longer nuclear interaction length. The acceptance and efficiency for p andp is much smaller than other charged hadrons.
Due to these species-dependent corrections, the overall acceptance and efficiency will depend on the relative production of these hadrons. In order to correctly account for the species dependence, an initial K ± /π ± ratio for each collision system is estimated separately. The contribution of p andp to reconstructed tracks based on this hadron simulation is less than 5%, and thus we do not include them in the overall result. Figure 6 shows the combined acceptance and efficiency for π ± and K ± as a function of p T in p+p collisions for various η ranges. The acceptance and efficiency is higher at more forward rapidity where path length through the absorber is shorter, and the total momentum of tracks for a given p T range is also larger. To have a more accurate correction, the full p T and η dependent correction is applied.

E. Nuclear modification factor
Nuclear effects on charged hadron production in p+Al and p+Au collisions are quantified with the nuclear modification factor, where dY pA /dp T dη is the charged hadron yield in a certain centrality bin of p+Al and p+Au collisions. These yields are corrected for the trigger efficiency, acceptance and reconstruction efficiency, and centrality bias factor introduced in Sec. II. dY pp /dp T dη is the hadron yield in p+p collisions corrected for the trigger efficiency, acceptance and reconstruction efficiency, and BBC efficiency. Finally N coll is the mean number of binary collisions for the corresponding centrality bin as calculated with the MC Glauber framework [47]. The N coll values, bias correction factors, and related systematic uncertainties for each centrality bin of p+Al and p+Au collisions appear in Table I.

IV. SYSTEMATIC UNCERTAINTIES
In this section, sources of systematic uncertainty in the nuclear modification factor are described, and the procedure used to determine each systematic uncertainty is discussed.
A. Acceptance and efficiency

Initial hadron distribution
Because there are limited measurements of identified charged hadrons at forward and backward rapidity (1.2 < |η| < 2.4), some model assumptions are necessary. Such forward rapidity particle yields have previously been estimated for use in earlier PHENIX p+p and d+Au collisions studies-see Ref. [12] for details. Here we follow that previous work as input for our simulation studies. To account for uncertainties on the estimated p T and η distributions, weight factors in p T and η for each collision system are extracted by comparing reconstructed p T and η distributions between data and simulation. The variation of acceptance and efficiency with modified initial p T and η distributions based on the weighting factors is less than 3% for p+p data. For p+Al and p+Au data, the variation at forward (backward) rapidity is less than 3% (5%). The variation is included in the systematic uncertainty.
In addition, there is an uncertainty in the K ± /π ± ratio which influences the combined acceptance and efficiency due to the longer nuclear interaction length of K + . Based on the uncertainties of measurements at midrapidity [11,43,44] used as an input for extrapolation to forward and backward rapidity and a possible extrapolation uncertainty estimated by comparing with the data at more forward rapidity [48], an effect of a ±30% variation of K ± /π ± on the acceptance and efficiency has been evaluated.
The K ± /π ± at midrapidity in various centrality bins of d+Au collisions are compatible with each other [11], and the difference of K ± /π ± between d+Au and p+Al and p+Au collisions in hijing [46] is less than 10%. These additional sources of uncertainty are covered by the 30% variation of K ± /π ± . The variation of acceptance and efficiency due to the 30% K ± /π ± change is less than 5% (7%) in p+p (p+Al and p+Au) collisions.

Proton contamination
As described in Sec. III D, the acceptance and efficiency is calculated for π ± and K ± . There is an ∼ 5% proton contamination where the fraction may vary with the initial p/(π + K) ratio. Based on the results in p+p and d+Au collisions at midrapidity [11,43,44], the p/π ratio at p T ≈ 2 GeV/c in 0%-20% central d+Au collisions is about 30% larger than in p+p collisions, which results in an increase of the contamination to 6.5% in 0%-20% central d+Au collisions as compared with 5% in p+p collisions. However, there is a lack of p/π measurements in a broader p T range in various centrality ranges of p+Al and p+Au collisions. Therefore, a conservative uncertainty of 5% is assigned corresponding to a factor of two difference in p/(π + K) ratios between p+Al, p+Au, and p+p collisions.

Hadron simulation
Although hadron response inside the absorber will not vary between different collision systems, the variation of acceptance and efficiency among three hadron interaction models (qgsp bert, qgsp bic, and ftfp bert) in geant4 has been checked. A detailed description of the three models and a previous study for muons can be found in Refs. [49,50]. The variation of the combined acceptance and efficiency for π ± and K ± between the three models is less than 2% in p T and η.

Variation of detector efficiency
During the data taking period, the detector performance varied due to temporary dead channels, changes in the instantaneous beam luminosity, and other experimental factors. The average detector efficiency for each collision system is included in the hadron simulation. The raw yield variation in FVTX and muon tracks is considered as a source of systematic uncertainty. The level of variation appears in Table II. The FVTX performance is quite stable during the entire data taking period, and the variation of the muon arm is observed to be larger in the south arm in the p+Au data due to a larger sensitivity of the MuID efficiency to the instantaneous beam luminosity of Au ions. A 1σ variation of the raw yield is assigned as a systematic uncertainty for each detector, and two systematic uncertainties are added in quadrature.

FVTX-MuTr mis-association
The probability of FVTX-MuTr mis-association depends on the FVTX track multiplicity, and the misassociation may artificially increase the acceptance and efficiency when requiring FVTX track association. The procedure for calculating the acceptance and efficiency using embedded simulations takes into account the multiplicity dependent FVTX-MuTr mis-association. The primary method to estimate the fraction of FVTX-MuTr is the data driven method described in Sec. III B, and the systematic uncertainty is evaluated by comparing with the estimated fraction from the embedded simulation. The difference is less than 1% of the maximum ∼ 3% of FVTX-MuTr mis-association contamination in 0%-5% p+Au collisions. A 1% systematic uncertainty is assigned for the estimation of FVTX-MuTr mis-association.

Vertex resolution
Because the location of the FVTX is close to the interaction point, the η acceptance of the FVTX depends on the z position of collisions. In the hadron simulation for acceptance and efficiency calculation, the measured z BBC distribution for each collision system is used, but there is uncertainty due to the resolution of z BBC . When considering the 2 cm of z BBC resolution, the variation of acceptance and efficiency is less than 0.5% in all three collision systems. A 0.5% systematic uncertainty is assigned due to the z BBC resolution.

B. Contamination from secondary hadrons
Remaining secondary hadrons can introduce a smearing of kinematic variables (p T and η) used in this analysis. The hadron simulation for calculating acceptance and efficiency already includes this component, however there can be a discrepancy in the relative contribution of secondary hadrons between the data and simulation. The systematic uncertainty on R pA is estimated by varying the FVTX-MuTr matching quality cuts (projection angles between FVTX and MuTr tracks) which affect the remaining fraction of secondary hadrons. Based on the hadron simulation, a tighter or looser FVTX-MuTr matching quality cut changes the relative fraction of secondary hadrons by ∼ 25%; the variation on R pA is less than 3%.

C. Multiple collisions
Due to the high instantaneous beam luminosity particularly in p+p and p+Al collisions, there is a chance of having multiple collisions in a single bunch crossing. This can introduce a bias in the yield calculation as well as centrality determination. The effect has been checked by analyzing two data groups with low and high instantaneous beam luminosity, and the difference in R pA is less than 5%. The variation due to multiple collisions is already considered in the systematic uncertainty from the variations in detector efficiency with data-taking period. Therefore, no additional systematic uncertainty is assigned.

D. BBC efficiency and centrality selection
The BBC efficiency in p+p collisions is ∼ 55% for MB events and ∼ 79% for hard scattering events, and a 10% systematic uncertainty is assigned based on previous studies [38]. This uncertainty is a global scale uncertainty.
As described in Table I, there are systematic uncertainties on N coll and bias correction factor calculations. The procedure to estimate these systematic uncertainties has been studied for d+Au collisions [39], and the same procedure is used for p+Al and p+Au collisions. Table III shows the summary of systematic uncertainties. All systematic uncertainties are point-to-point correlated. Because most of sources on the acceptance and efficiency are independent in each collision system, there is no cancellation of systematic uncertainty for R pA calculation.  0%-100% centrality are obtained by integrating over all centrality and applying the bias correction factors. Bars (boxes) around the data points represent statistical (systematic) uncertainties, and boxes around unity represent the global systematic uncertainty due to uncertainties in the BBC efficiency and the calculated N coll . The results for p+Al indicate that there is little modification at forward rapidity (i.e. in the p-going direction), whereas a small enhancement is observed in p T < 2 GeV/c at backward rapidity (i.e. in the Al-going direction). In p+Au results, a suppression is seen in p T < 3 GeV/c at forward rapidity unlike the p+Al results. At backward rapidity, a similar trend of enhancement is observed in the p+Au data, though with larger magnitude. GeV. Also shown are comparisons to a pQCD calculation [14] and calculations based on the nPDF sets [22,23].

E. Summary of systematic uncertainty
Comparisons with estimated R pA based on nuclear modified PDFs are shown from the ncteq15 nPDF [22] and the epps16 nPDF [23] interfaced with pythia v8.235 [51]; the parameters used in the event generation of pythia are listed in Table IV. Note that the multiplication factor for multiparton interactions is determined by comparing the η-dependent multiplicity distribution in p+p collisions at √ s = 200 GeV [52]. The calculations indicate a modest expected suppression at forward rapidity from shadowing of low-x Bj partons in the Au nucleus, and are in agreement with the data within uncertainties. However, at backward rapidity, sensitive to potential anti-shadowing of higher-x Bj partons in the Au nucleus, the calculations result in no modification in contradistinction from the data. pQCD calculations considering incoherent multiple scatterings inside the nucleus  RpA of charged hadrons as a function of pT at (a) forward and (b) backward rapidity in p+Au 0%-100% centrality selected collisions at √ s N N = 200 GeV. Also shown are comparisons to a pQCD calculation [14] and calculations based on the nPDF sets [22,23].
before and after hard scattering [14] at backward rapidity are also compared with the data, and it agrees with the both p+Al and p+Au data. Figure 9 shows R pA of charged hadrons integrated over the interval 2.5 < p T < 5 GeV/c as a function of η in the 0%-100% centrality selection of (a) p+Al and (b) p+Au collisions at √ s N N = 200 GeV. Again the data are compared with pQCD calculations at backward rapidity and GeV. Also shown are comparisons to a pQCD calculation [14] and calculations based on the nPDF sets [22,23].
calculations based on two nPDF sets. In p+Au collisions, there is a modest hint that enhancement at backward rapidity becomes larger as η approaches midrapidity, while the suppression at forward rapidity becomes stronger. In p+Al collisions, R pA at forward rapidity is quite similar to what is observed in p+Au collisions, whereas it shows a smaller enhancement at backward rapidity than the results in p+Au collisions. The comparison with ncteq15 and epps16 nPDF calculations indicates that the R pA at forward rapidity agrees in both p+Al and p+Au collisions, but the enhancement at backward rapidity in p+Au collisions is not reproduced by the both calculations. In case of the comparison with the pQCD calculations at backward rapidity, the magnitude of enhancement is similar. However, the pQCD calculations show a stronger enhancement at more backward rapidity which is different from the trend in the data. Because initial and final-state nuclear effects on hadron production may depend on the density of initial partons in the nucleus and on the density of final-state produced particles, R pA has been measured in various centrality bins of p+Al and p+Au collisions. Figures 10 and 11 show R pA of charged hadrons as a function of p T or η at forward and backward rapidity from the most central bin (0%-5%) to the most peripheral bin (40%-72%) for p+Al collisions at √ s N N = 200 GeV. The results at forward and backward rapidity are plotted together in each plot. First, there is a clear centrality dependence both at forward and backward rapidity. The magnitude of the modification, which shows enhancement at backward rapidity and suppression at forward rapidity, becomes stronger in more central p+Al collisions. The observed R pA in the most peripheral (40%-72%) p+Al collisions is consistent with unity in both rapidity regions, indicating little modification of charged hadron production compared to the p+p data. Both the magnitude of the modification and the p T dependence are larger in central collisions. at forward and backward rapidity in central p+Au collisions. The centrality dependence of R pA as a function of η shown in Fig. 11 is consistent with what is seen in R pA as a function of p T . The η dependence at backward rapidity is weakly centrality dependent, but there is a clear η dependence at forward rapidity in the most central collisions. Figures 12 and 13 show R pA of charged hadrons as a function of p T and η in various centrality classes of p+Au collisions. Similar to the results in p+Al collisions, the magnitude of modification becomes larger in more central collisions both at forward and backward rapidity, and the R pA values in the most peripheral p+Au collisions are consistent with unity. When comparing p+Al and p+Au results in the 0%-5% central collisions shown in the panel (a) of Figs. 10, 11, 12, and 13, R pA at forward rapidity is comparable between the two collision systems. However, the enhancement at backward rapidity is much stronger in p+Au collisions. Figure 12 compares pQCD calculations with the p+Au data at backward rapidity. Similarly with the comparison in the integrated centrality, the calculation can reproduce the p T and centrality dependent enhancement. Figure 14 shows R pA as a function of N part for charged hadrons in the range 2.5 < p T < 5 GeV/c at (a) forward and (b) backward rapidity in p+Al and p+Au collisions at √ s N N = 200 GeV. Unlike the previous results, the systematic uncertainty on N coll is included in boxes around data points. The data show that R pA at backward rapidity (filled [black] circles), i.e. in the A-going direction, increases monotonically with N part , and the trend is reproduced by the pQCD calculation. However, R pA at forward rapidity (open [red] circles), i.e. in the p-going direction, reveal that each collision system has its own decreasing trend as N part becomes larger. R pA at forward rapidity in 0%-5% of p+Al and p+Au collisions are consistent (R pA ∼ 0.7), although N part (9.7 in p+Au and 4.1 in p+Al collisions) are quite different. The trend of a larger enhancement (suppression) at backward (forward) rapidity in more central collisions is consistent with the previous results of charged hadrons and muons from heavy flavor decay in d+Au collisions [12,31]. A closer look on η-dependent R pA in 0%-5% p+Al and 40%-60% p+Au collisions of similar N part is shown in Fig. 15. At backward rapidity, it shows not only a consistent magnitude of R pA but also a quite similar trend of R pA in η. In case of the comparison at forward rapidity, R pA of the 40%-60% p+Au centrality bin is consistent with unity in all η bins, whereas a η-dependent suppression is seen in 0%-5% p+Al collisions.
The suppression of charged hadron production at forward rapidity in integrated centrality of p+Al and p+Au collisions can be explained by the nPDF modification based on the comparison with the ncteq15 and epps16 calculations shown in Figs 7, 8, and 9. It would be useful to extend another calculation within the CGC framework [27], which successfully describes the suppression of charged hadron production at forward rapidity in d+Au collisions [8,9]. More differential calculations from these various frameworks are needed to compare to the systematic trends found in our new results. In addition to these models which consider modification of the parton distribution functions inside the nucleus, the pQCD calculation of dynamic shadowing considering coherent multiple scatterings inside the nucleus [19] also predicts a rapidity and impact parameter dependent suppression of hadron production at forward rapidity. The centrality dependent suppression at forward rapidity shown in both p+Al and p+Au collisions also can be described by the color fluctuation effects expecting a stronger centrality dependence in p+Au collisions than d+Au collisions [30]. It will be quite useful to have theoretical calculations for detailed comparison with the data in p T , rapidity, and centrality.
For the enhancement of charged hadron production observed at backward rapidity, estimates from the nPDF sets clearly fail to describe the data. A pQCD calculation considering incoherent multiple scatterings both before and after hard scattering [14], which can describe the enhancement of heavy quark production at backward rapidity in d+Au collisions [12], successfully explains the centrality and A-dependent enhancement. In addition, there is also a possibility of hydrodynamic behavior showing a larger elliptic flow of charged particles at backward rapidity where the multiplicity is also larger than other rapidity ranges [53].   are comparisons to a pQCD calculation [14].

VI. SUMMARY
PHENIX has measured the nuclear modification factor R pA of charged hadrons as a function of p T and η at forward and backward rapidity in various centrality ranges of p+Al and p+Au collisions at √ s N N = 200 GeV.
The results in central p+Al and p+Au collisions show a suppression (enhancement) in the forward p-going (backward, A-going) rapidity region compared to the binary scaled p+p results of 0.7 (2.0) for p+Au and 0.9 (1.2) for p+Al in 2.5 < p T < 5 GeV/c at a level of significance 3.3σ (3.2σ) for p+Au and 2.7σ (1.1σ) for p+Al. In contrast, there is no significant modification of charged hadron production observed in peripheral p+Al and p+Au collisions in either rapidity region. The enhancement at backward rapidity shows a clear A-dependence, but the suppression at forward rapidity is comparable between the two collision systems despite more than a factor two larger N part in p+Au collisions. The results integrated over centrality are compared to a calculation with the ncteq15 and epps16 nPDF sets. The calculation agrees with the data at forward rapidity both in the integrated centrality of p+Al and p+Au collisions, but it fails to describe the enhancement observed at backward rapidity in p+Au collisions. Because the nPDF sets does not yet provide an impact parameter dependent nPDF, the comparison is limited to the case of integrated centrality. These data measured in various centrality ranges can be useful to test impact parameter dependent nPDFs in different nuclei in the future. The pQCD calculation considering incoherent multiple scatterings inside the nucleus can describe the data at backward rapidity. In addition, a comparison with different models can help to improve the understanding of nuclear effects in small collision systems.