Impact-parameter dependent nuclear parton distribution functions: EPS09s and EKS98s and their applications in nuclear hard processes

We determine the spatial (impact parameter) dependence of nuclear parton distribution functions (nPDFs) using the $A$-dependence of the spatially independent (averaged) global fits EPS09 and EKS98. We work under the assumption that the spatial dependence can be formulated as a power series of the nuclear thickness functions $T_A$. To reproduce the $A$-dependence over the entire $x$ range we need terms up to $[T_A]^4$. As an outcome, we release two sets, EPS09s (LO, NLO, error sets) and EKS98s, of spatially dependent nPDFs for public use. We also discuss the implementation of these into the existing calculations. With our results, the centrality dependence of nuclear hard-process observables can be studied consistently with the globally fitted nPDFs for the first time. As an application, we first calculate the LO nuclear modification factor $R^{1jet}_{AA}$ for primary partonic-jet production in different centrality classes in Au+Au collisions at RHIC and Pb+Pb collisions at LHC. Also the corresponding central-to-peripheral ratios $R_{CP}^{1jet}$ are studied. We also calculate the LO and NLO nuclear modification factors for single inclusive neutral pion production, $R_{dAu}^{\pi^0}$, at mid- and forward rapidities in different centrality classes in d+Au collisions at RHIC. In particular, we show that our results are compatible with the PHENIX mid-rapidity data within the overall normalization uncertainties given by the experiment. Finally, we show our predictions for the corresponding modifications $R_{pPb}^{\pi^0}$ in the forthcoming p+Pb collisions at LHC.


Introduction
In a high-energy hadronic or nuclear collision of particles A and B the inclusive cross sections for hard processes involving a large interaction scale Q 2 ≫ Λ 2 QCD can be computed using the QCD collinear factorization theorem [1,2], where dσ represents the perturbatively computable partonic pieces (cross sections in lowest order), and f A i (f B j ) is the parton distribution function (PDFs) for a given parton flavor i in the colliding particle A (and correspondingly for the flavor j in B). The PDFs are universal, process-independent functions of nonperturbative origin, whose evolution in the scale Q 2 can, however, be obtained from the DGLAP equations [3][4][5][6] derived from perturbative QCD. A precise knowledge of the universal PDFs is thus vital for interpreting any hardprocess results at the present colliders BNL-RHIC and CERN-LHC. This holds as well for proton-proton collisions as for proton-nucleus and nucleus-nucleus collisions. To determine the nonperturbative input in the PDFs, one has developed global analyses which exploit a multitude of experimental hard-process data and the DGLAP evolution. Excellent fits for the free proton PDFs have been obtained, and sets like CT10 [7], MSTW [8], and NNPDF2.0 [9] are nowadays available.
It is well known that the PDFs of nucleons bound to a nucleus, the nuclear PDFs (nPDFs), are modified relative to the free-nucleon PDFs. Analogously to the free-proton case, global DGLAP analyses have been developed also for the nPDFs. The further complication with these is that in addition to the usual x and Q 2 dependences also the nuclear mass-number (A) dependence of the PDFs needs to be dealt with. The global nPDF fits have so far resulted in leading-order (LO) nPDF sets EKS98 [10], HKM [11] and HKN04 [12], and next-to-leading order (NLO) sets nDS [13], HKN07 [14], EPS09 [15], nCTEQ [16,17], and DSZS [18]. Importantly, and similarly to the free proton case, with the error sets of EPS09 (and similar sets in DSZS), one can nowadays quantify how the uncertainties remaining in the nPDFs, illustrated in Fig. 1, are transmitted to the nuclear hard-process cross-sections.
The global analyses mentioned above have all considered only the spatially averaged nPDFs, probed in minimum-bias nuclear collisions with no cuts on the collision centrality (impact parameter). In particular, as the modest amount of available nuclear hard-process data severely limits the number of possible fit parameters, it has so far been impossible to embed the spatial dependence, or the impact-parameter dependence, of the nPDFs directly into the global analysis. An obvious drawback with the globally analysed nPDFs then is that it has not been possible to consistently compute nuclear hard-process cross-sections in different centrality classes.
The purpose of this study is to consider this problem by pinning down the spatial dependence of the nPDFs, i.e. the dependence of the nuclear modifications of the PDFs on the nucleon's position in a nucleus. We do this in a manner which is for the first time fully consistent with the nPDFs from a global analysis. Earlier attempts to this direction, lacking however such a consistency, can be found in Refs. [19][20][21][22]. A further motivation for the current study is the Gribov-Glauber modeling of nuclear shadowing, reviewed lately in Ref. [23], whose output nPDFs are not a result of a global analysis like EPS09 but which have so far been the only ones where the spatial dependence arises in a self-consistent manner from modeling the physics origin of the nuclear effects. On the experimental side, the current study is inspired by e.g. the measurements of single hadron production [24][25][26][27][28][29][30] and J/Ψ production [31,32] in different centrality classes in d+Au collisions at RHIC, as well as by the hard-process measurements in the forthcoming p+Pb collisions at the LHC. Also the theoretical modeling of the J/Ψ production discussed recently in Ref. [33] has motivated our study. Our basic idea for uncovering the spatial dependence in the EKS98 and EPS09 nPDFs is straightforward: We first introduce the spatial dependence of the nuclear modification to the nPDF of each parton type i in each nucleus A at each x and Q 2 in terms of a power series of the standard nuclear thickness functions T A . Then, we determine the coefficients of each power of T A by exploiting the A-dependence of the EPS09 and EKS98 nPDFs (these sets, through the global fits, represent the experimental data here). As an output, we provide the numerical routines named EPS09s (LO and NLO as well as error sets for both) and EKS98s for computing the spatially dependent nPDFs which -simultaneously for all nuclei considered -normalize to the corresponding spatially independent EPS09 and EKS98 nPDFs. These new sets will be downloadable at the link [34].
As concrete examples of how to easily implement our spatially dependent nPDFs and the nuclear collision geometry in the computation of nuclear hard-process cross-sections in different centrality bins, we first discuss the centrality dependence of the LO nuclear modification ratios R 1jet AA (p T ) of primary partonic-jet production in Au+Au collisions at RHIC and Pb+Pb collisions at the LHC. We also study the nuclear modification factors of inclusive π 0 production, R π 0 dAu , in d+Au collisions at RHIC and R π 0 pPb in p+Pb collisions at the LHC, both at mid-and forward rapidities, and considering both the NLO and LO cases. For R π 0 dAu we also make, to our knowledge, a first comparison with the PHENIX centrality dependent data [26] where the overall normalization errors of the data are accounted for in detail. Due to the planned p+Pb program at the LHC, the ratio R h pPb (p T ) for single hadron production has been of growing interest recently [35][36][37][38][39][40], and we will show also here how interesting and useful this ratio would be from the point of view of constraining the nPDFs further.
The paper is organized as follows: In Sec. 2 we define the model framework and explain the fitting procedure. In Sec. 3 we show the results for the spatially dependent nuclear modifications of PDFs. Also a comparison with selected other works is presented here. Applications of our results are discussed in Sec. 4. For clarity, a summary of the standard elements used in the applications here, the formulation of the nuclear collision geometry, different overlap functions and the optical Glauber model, is given in the Appendix A.

Definitions of the Nuclear Modifications
First we need to define how we introduce the spatial dependence to the nPDFs in terms of the hard-process cross-sections. Let us start with the usual spatially averaged nPDFs. The number distribution of an observable k produced in a collision of nuclei A and B at an impact parameter b is given by where T AB (b) is the standard nuclear overlap function normalized to AB (cf. Eq. (A. 20) in App. A.2, see the nuclear collision geometry in Fig. 20), and dσ AB→k+X is the bindependent inclusive hard cross-section of Eq. (1.1) containing the nPDFs and perturbative pieces. The spatially averaged nPDFs in a nucleus A with Z protons and A − Z neutrons are now given by where the nPDFs of a bound neutron, f n/A i , may be (approximately) obtained from those of the bound proton, f p/A i , by using the isospin symmetry (see [15]). As in EKS98 and EPS09, we define the nPDF for each parton flavor in terms of the spatially averaged nuclear modification R A i (x, Q 2 ) and the corresponding free proton PDF f p i (x, Q 2 ), To lighten the notations, we express the nPDFs in Eq. (2. 2) as where the sum runs over all the nucleons N = 1, . . . , A.
Decomposing the T AB into the standard nuclear thickness functions (cf. Eq. (A.20)), and using Eq. (1.1) we may write (2.5) From this, we see that a suitable definition of the spatially dependent nuclear modification r A i (x, Q 2 , s) for the PDF of parton flavor i (per nucleon) is where the thickness function T A is normalized to A and where the case of no nuclear effects corresponds to R A i = r A i = 1. Using these definitions, we can now generalize Eq. (2.5) to include the spatially dependent nuclear modifications, As a consistency check, we note that the definition in Eq. (2.6) guarantees that the minimum-bias cross sections, which are obtained by integrating Eq. (2.7) over the whole b space, become simply AB times the hard cross-section computed with the spatially averaged nPDFs, The key assumption in the present analysis is that the spatial dependence of r A i (x, Q 2 , s) is a function of the nuclear thickness T A (s). The motivation for this comes mainly from the shadowing region at small x, where the partons of sufficiently small values of x may interact with partons from any other nucleon near enough in the transverse direction. Also in the Gribov-Glauber modeling [23] of the initial state nPDFs the nuclear effects become essentially functions of T A . The functional form we choose to use and test here is a simple power series of the thickness functions, Here we would like to emphasize the following points: First, all the A dependence is now in the thickness functions which are fully known, and all the coefficients c i j (x, Q 2 ) which will be our fit parameters, depend on x and Q 2 but not on A. Second, the power series of the form 1 + . . . also fixes by construction that r A i (x, Q 2 , s) → 1 when |s| → ∞, which means that the nucleons at the very edge of the nucleus are essentially regarded as free nucleons. Third, what is known from the EKS98 and EPS09 -types of analyses, are only the spatially averaged nuclear modifications and their A systematics, i.e. T A -weighted integrals of Eq. (2.9) over s for each nucleus. Fourth, since the EKS98 and EPS09 global analyses have not been constructed to reproduce any specific theoretically motivated A dependence of the nPDFs, we can test the validity of the assumption of Eq. (2.9), as well as the number of terms needed, only a posteriori. Using the definitions above, we can see why the simplest 1-parameter approach with n = 1 in Eq. (2.9) (which is used e.g. in [19][20][21][22] as well as in e.g. the HIJING event generator [41]) is not fully consistent with the observed A systematics of the nuclear data. In this case, , and from the definition in Eq. (2.6), one where R A i is given by the globally analysed nPDFs (i.e. nuclear data). The problem then is that the coefficient c i (x, Q 2 ) may depend in fact quite strongly on A, which indicates that the simplest assumption of terminating the power series at the first nontrivial term does not correctly capture the spatial dependence of the measured nuclear structure functions. This redundant A dependence is illustrated in Fig. 2 for gluons in a lead nucleus at x = 0.01 at the initial scales of the sets EKS98, EPS09LO1 and EPS09NLO1. We can see that especially for the NLO set the problem is more serious. One of the driving motivations for the present study is to solve the problem of recovering the A systematics in the spatially dependent nPDFs.

Fitting Procedure
To extract the A-independent coefficients c i j (x, Q 2 ), we need to introduce a fitting procedure, where we utilize the definition (2.6) and the A dependence of the EKS98 and EPS09 nuclear modifications at different values of x and Q 2 for each parton flavor i. To reproduce the A systematics in the spatially independent nuclear modifications with the power-series ansatz of Eq. (2.9), we minimize the χ 2 defined as where the spatially averaged modifications R A i (x, Q 2 ) from EKS98 and EPS09 now represent the "experimental" data. The weight factors W A i (x, Q 2 ) are artificial errors which control the quality of the fit and which are set by hand. Our numerical observation is that for good fits we need 4th-order polynomials in T A , i.e. n = 4 in Eq. (2.9). Furthermore, best fits were obtained with the weight W A i (x, Q 2 ) = R A i (x, Q 2 )− 1 for EKS98 (this corresponds to fitting the deviations from unity within a constant relative error) and W A i (x, Q 2 ) = 1 for EPS09 (corresponds to fitting the modifications within a constant error).
By construction, both EKS98 and EPS09 give no nuclear modifications for deuterium. This cannot be reproduced with the functional form we selected for r A i (x, Q 2 ), and we do not expect the fit form of Eq. (2.9) work for the smallest values of A, either. Consequently, we exclude the nuclei A < 16 from the fit. Thus for EKS98 the sum runs over A = 16, 20, . . . , 300 (i.e. emphasising the large nuclei) and for EPS09 we use all the A ≥ 16 values for which these sets are currently available.

Quality of the fit
First, we demonstrate that our fit framework manages to reproduce the spatially averaged nuclear modifications and especially their A dependence indeed very well. Figure 3 shows the obtained spatially dependent gluon modifications integrated over the transverse plane according to Eq. (2.6), and the corresponding input modifications at different fixed values of Q 2 from the NLO set EPS09NLO1 (left panel), and from the LO sets EPS09LO1 and EKS98 (right panel), for a lead nucleus. In what follows, we refer to these cases as "EPS09sNLO1", "EPS09sLO1" and "EKS98s" where "s" is for "spatial" and "1" for the central sets. As seen in the figure, the match with the input and output distributions is very good; for all parton flavors and the nuclei included in our fits it is within 2 % at x < 0.75 for EPS09NLO, 1 % at x < 0.85 for EPS09LO, and 0.2 % at x < 0.95 for EKS98. Importantly, the keyfeature here, the A dependence of EPS09 and EKS98, is similarly well reproduced, as is demonstrated by Fig. 4 below.
Recall also that in the EPS09 global analysis in addition to the best fit there are also 30 error sets, which enables one to compute how the uncertainties of the nPDFs propagate into physical observables. The above fitting and determination of the spatial dependence are done also for each of these error sets, both in LO and in NLO, and the fit quality is similar as in Figs. 3 and 4. Thus the error propagation calculations (as instructed in [15]) for centrality-dependent nuclear hard cross-sections can now be done as before, using the EPS09s sets.   Figure 3. Left: The spatially averaged nuclear modification R A g (x, Q 2 ) for a lead nucleus (A = 208) from the NLO set EPS09NLO1 (dotted lines) and from the EPS09sNLO1 spatial fit presented here (solid lines) at four different scales. Right: The same with the LO sets EKS98 and EPS09LO1 (dotted) and with the spatial fits EKS98s (dashed) for three different scales and EPS09sLO1 (solid) for four different scales.   The A dependence of the spatially averaged nuclear modification R A g (x, Q 2 ) at fixed values x = 0.001 and Q 2 = 1.69 GeV 2 from the sets EPS09NLO1 (crosses) and EPS09LO1 (pluses) and from the corresponding spatial fits EPS09sNLO1 (solid green line) and EPS09sLO1 (solid blue line). Right: The same but with the LO set EKS98 (circles) and the corresponding spatial fit EKS98s (solid red) at Q 2 = 2.25 GeV 2 . The small nuclei shown with gray markers in both panels were not used in our spatial fits.

Spatial Dependence
After the consistency checks above, let us next discuss the spatial dependence obtained for the nuclear modifications of the PDFs. In Fig. 5 we present the nuclear modification r Pb g (x, Q 2 , s) at the initial scale Q 2 = 1.69 GeV 2 as a function of x and s, as obtained from the fitting to the sets EPS09 NLO and LO, as well as the LO set EKS98. The three main observations are • The overall x-shape of the nuclear modification away from the edge of the nucleus, at |s| < R A , is similar as in the input distribution. This confirms that our fit does not generate any unwanted extra curvature.
• In the center of the nucleus, |s| ≈ 0, the nuclear modification is only slightly larger than the input average modification. This also confirms the earlier similar findings in [19,20,22].
• The nuclear modification dies out as expected, by construction, when |s| > R A . This feature arises from the vanishing T A (s) at the edge of the nucleus.
The observations for the spatial dependence of the sea and valence quarks nuclear modifications are the same. Examples of these can be found in App. B.

Comparison with other approaches
Next, we compare our EPS09s and EKS98s fits with the 1-parameter approach described in the end of Sec. 2.1. [20,22]. 1 This model has been used to study the centrality dependence of the J/Ψ suppression e.g. in Refs. [31,33,42,43] 2 and inclusive hadron production in d+Au collisions at RHIC in Ref. [22]. We also compare with the leading-twist formulation [23,44] of nuclear shadowing which is based on the generalization of Gribov-Glauber theory, QCD factorization and diffractive PDFs measured at HERA. For the spatially averaged nuclear modifications, this model typically predicts a stronger smallest-x shadowing than what is implemented in the parametrizations of EKS98 and EPS09 (see e.g. Ref. [23]). For the comparison, we consider the FGS10 L set [23,45], and choose the value of x not too small, so that the spatially averaged FGS10 L nuclear gluon modification is close to that in EPS09 or EKS98. In Fig. 6 we plot the nuclear modification for gluons at fixed values of x and scale Q 2 = 4 GeV 2 for A = 208 as a function of |s| from our EPS09sNLO1, EPS09sLO1 and EKS98s fits, from the 1-parameter approach using the averaged sets EPS09NLO1, EPS09LO1 and EKS98, and from FGS10 L. Although numerically the differences are not very large, we notice that while both the EPS09sNLO and EKS98s results are close to FGS10 L, the 1-parameter approach leads to a too steep transverse profile for the modifications in all cases.

Applications
Next, we consider some concrete examples of computing the nuclear hard-process crosssections in different centrality classes using the spatially dependent nPDFs. First, we discuss the centrality dependence of primary partonic-jet production in A+A collisions at RHIC and LHC. Then, we consider neutral pion production in d+Au collisions at RHIC and in p+Pb collisions at the LHC. 1 In [20] the spatial dependence enters through the first nontrivial power of the nuclear density ρA(r) or the thickness function TA(s). The latter scenario corresponds to what we refer to as "1-parameter approach" here.
2 In [33] one studies also other types of spatial dependences.

Implementation of EKS98s and EPS09s
For defining the centrality classes we use the optical Glauber model specified in App. A.2.
In this case, each centrality class corresponds simply to a certain impact-parameter interval The generic average number distribution of a hard-process observable k in this centrality class of an A+B collision is  where p inel 23), and dN k AB (b) is obtained from Eq. (2.7). Using the expansion of r A i (x, Q 2 , s) in powers of T A from Eq. (2.9), the integrals over the impact parameter for the spatially dependent parts can be conveniently separated from the spatially independent fit coefficients, free nucleon PDFs and pQCD parts as follows: are the free nucleon PDFs, and we have defined c i,j 0 (x, Q 2 ) ≡ 1 and From Eq. (4.2) we see the most straightforward implementation of the spatially dependent nPDFs. The purely geometric integrals, the T nm AB (b 1 , b 2 ) in Eq. (4.3) for each pair of the powers n and m, can be computed independently of the kinematic variables x 1 , x 2 , and Q 2 and also independently of the parton flavors i, j. Thus, in total we have 25 different geometric integrals to do (or 15 if A = B) but we need to do them only once. In comparison with the spatially averaged case, the fit parameters c i n (x, Q 2 ) thus play the role of the nuclear modifications R i A (x, Q 2 ) for each of the 25 pairs n, m. To arrive at the final b-integrated result for the number distribution of k, we thus need to repeat the computation of the kinematic parts 25 times, each with different sets of the coefficient . The EKS98s and EPS09s routines which we provide in [34], give in addition to the fit coefficients {c i n (x, Q 2 )} also the thickness functions T A (s) (used in the fits here) for the computation of T nm AB (b 1 , b 2 ), as well as the combination T A (s)r A i (x, Q 2 , s) for other possible implementations. Note also that for the b integral in Eq. (4.3) the angular part is trivial, giving just 2π.

The Nuclear Modification Factors R 1jet
AA and R 1jet

CP
Let us now consider the centrality dependence of primary inclusive high-p T parton production in A+A collisions at RHIC and LHC. Following the generic discussion above, we define the nuclear modification ratio R 1jet AA (p T ) relative to the p+p case for each centrality class as where N AA bin b 1 ,b 2 is the average number of binary collisions in this centrality class given by Eq. (A.27) and σ N N inel is the inelastic nucleon-nucleon cross section. Apart from the (small) isospin effect, this ratio yields unity if there are no nuclear effects in the nPDFs. Thus, for peripheral enough centrality bins, this ratio should approach unity. For the details of the partonic cross sections, bookkeeping and kinematics, we refer to [46].
The nuclear mofication factor in the minimum-bias collisions is obtained from above by setting b 1 = 0 and b 2 → ∞, in which case we have where dσ 1jet AA,MB , which contains only the spatially averaged nPDFs, is obtained from Eq. (2.8) by setting B = A, and the p+p baseline dσ 1jet pp from the same equation by setting A = B = p.
In addition to the centrality dependence of R 1jet AA , we are interested in the central-toperipheral ratios, defined as where C and P refer to the central and peripheral bins, correspondingly. The advantage of this ratio (in the experiments) is that the information of the proton-proton baseline is not required. In particular, we would like to see exactly how much R 1jet CP differs from the modification R 1jet AA which is computed with the spatially averaged nPDFs. We perform these example-calculations for both RHIC and LHC but for simplcity only to LO pQCD, since without jet quenching these ratios do not directly correspond to observables. They illustrate, however, the points we wish to make with the spatially dependent nPDFs, and also serve as (LO) pQCD baselines for the observed suppression of high-p T particles.
The two different centrality classes we consider here for Au+Au collisions at RHIC and Pb+Pb collisions at the LHC, are the central 0-5% and peripheral 60-80% bins. The Glauber model input and the resulting impact parameter intervals and average numbers of binary collisions in these centrality classes are summarized in Table 1. In Fig. 7 we plot the ratio R 1jet AA (p T , y = 0) for central, peripheral and minimum-bias collisions, as well as R 1jet CP (p T ) in Au+Au collisions at RHIC. Figure 8 shows the same quantities for the LHC Pb+Pb case. The central and peripheral R 1jet AA and R 1jet CP have been obtained with the spatially dependent nPDFs EPS09sLO1 (left) and EKS98s (right), and the average R 1jet AA in minimum bias collisions with the spatially independent EPS09 and EKS98 nuclear modifications. For the free proton PDFs we have used CTEQ6.1L [47]. The renormalization scale µ and factorization scale Q has been set to be the transverse momentum, p T , of the parton.
The main observations from the figures are: (i) The central R 1jet AA is quite close to the average R 1jet AA , which is expected since the nuclear modifications at small s are close to the average modifications. (ii) The peripheral R 1jet AA is clearly not unity but there appear almost 10% antishadowing effects at mid-p T at RHIC and even more than 20% shadowing effects at small p T at the LHC, and up to 10% EMC effects at large p T both at RHIC and LHC. [GeV/c] [GeV/c]    suggest that in a precision theory-analysis of the centrality dependence of jet quenching, one needs to account also for the spatial dependence of nPDFs. Finally, regarding the differences between the different nPDF sets applied here, we observe in Figs. 7 and 8 how the stronger shadowing in the EPS09 case (cf. Fig.3) translates into steeper p T slopes of R 1jet AA at small p T than in the EKS98 case.
4.3 Centrality dependence of R π 0 dAu (p T ) at RHIC -comparison with data While the above ratios R 1jet AA and R 1jet CP mainly serve as theoretical pQCD baselines for jet quenching studies, it is important to test our spatially-dependent nPDF framework against some measured centrality-dependent observables. To avoid the complications of hot QCD matter modeling, we turn to the highest-energy d+Au collisions at RHIC and p+Pb at the LHC. For our purposes a promising published data set is the nuclear modification factor R π 0 dAu (p T ) for single inclusive neutral-pion production at mid-rapidity |η| < 0.35, measured by PHENIX [26] at different centrality classes in d+Au collisions at √ s N N = 200 GeV. Since the minimum-bias data precisely from this data set was among the constraints in the global EPS09 fits, it is now very interesting to study, for the first time consistently with EPS09, how well we can reproduce the measured centrality dependence of this ratio.
Analogously with Eq. (4.4), we define the centrality-dependent nuclear modification factors as where the number distribution now involves a further folding over the fragmentation functions, and where dN k dA (b) is obtained from Eq. (2.7) by setting A = d and B = A, and (as we do not assign any nuclear effects for the deuterium PDFs) also r d i ≡ 1. For obtaining a realistic thickness function for deuterium, we use the Hulthen wavefunction formulation [48] given in App. A.1.2. The impact parameter ranges and average numbers of binary collisions at the corresponding centrality classes are obtained again from the optical Glauber model.

Minimum-bias R π 0
dAu (p T ) Setting the spatial integrals in Eq. (4.7) over the whole impact-parameter space gives again the minimum-bias ratios, where dσ π 0 dA,MB = k dσ k dA,MB ⊗ D π 0 /k (z, Q 2 F ) again contains only the spatially averaged nPDFs in dσ k dA,MB which is obtained from Eq. (2.8) by setting A = d, B = A. As noted earlier, in the EKS98 and EPS09 frameworks there are no nuclear modifications to the deuterium PDFs. The p+p baseline dσ π 0 pp is computed correspondingly, but without any nuclear effects. Figure 9 shows the PHENIX data [26] and our NLO (left) and LO (right) results for the nuclear modification factor R π 0 dAu (p T , y = 0) in minimum-bias collisions. For the NLO calculation with EPS09sNLO1 (equivalently one may use EPS09NLO1, since the spatial dependence is here irrelevant) we used the NLO fragmentation functions from KKP [49] 3 , AKK [50] and fDSS [51]. For the free proton PDFs, we use CTEQ6M [47]. Correspondingly, the LO case was computed with EPS09sLO1 and EKS98s, using the KKP and fDSS LO fragmentation functions and CTEQ6.1L PDFs [47]. The renormalization scale µ, factorization scale Q and fragmentation scale µ F are all fixed to p T , the transverse momentum of the produced hadron. For details of the LO calculation, we again refer to Ref. [ Figure 9. The nuclear modification factor R π 0 dAu (p T ) for √ s N N = 200 GeV at y = 0 for minimum bias collisions. Calculations are for the NLO pQCD using EPS09NLO1 with three different fragmentation functions (left), and LO using EPS09LO1 and EKS98 with two different fragmentation functions (right). The blue bands are computed using the EPS09 error sets and fDSS fragmentation functions. The experimental PHENIX data [26] are shown by markers, and their error bars (boxes) stand for the point-to-point statistical (systematic) errors. Notice that the data points and their errors have been multiplied by a factor 1.039 (1.050) for the NLO (LO) case, which is well within the 9.7 % overall normalization error quoted by the experiment (see text for details).
From Fig. 9, we notice the following: (i) The EPS09 uncertainty bands for the NLO results are slightly smaller than in the LO case, reflecting the fact that the EPS09NLO gluons are somewhat better constrained in the antishadowing region than those of EPS09LO.
(ii) In the small p T region there is a difference in the p T slopes between the EKS98 and EPS09LO1 results. This is caused by the weaker shadowing in EKS98. However, also the EKS98 results remain within the EPS09 error bars. (iii) The uncertainty caused by the differences in the fragmentation functions remains conveniently small in all cases.
Regarding the data comparison in Fig. 9, we emphasize the following important point: In addition to the the statistical uncertainties (error bars) and point-to-point systematic errors (boxes), PHENIX quotes a 9.7 % overall uncertainty which originates from the p+p reference and which is not included in the statistical error bars shown. Consequently, allowing for a shift of the data points and their errors by less than 9.7 % and requiring the best possible overall fit to the data (using the fDSS FFs and by minimizing the χ 2 with the point-to-point statistical and systematic errors added in quadrature), we have multiplied the data by a factor 1.039 (NLO) and 1.050 (LO). Such a few-percent shift is well within the uncertainty given by the experiment. As already noticed in the EPS09 analysis [15], the resulting agreement with the data is quite good, both in LO and in NLO. Figure 10 shows the ratio R π 0 dAu (p T , y = 3) at the forward region in minimum-bias collisions. Again, we show both the NLO (left) and LO (right) results with the same set-up as in Fig. 9 above. Now the differences between the fragmentation functions start to be visible, as one is probing their larger-z tails where the uncertainties are larger: for the same pion p T , the differences in the large-z fragmentation functions map to different values of x in the nPDFs. We also notice that the LO calculation gives a stronger small-p T suppression than the NLO case, which is partly due to the different pace of the scale evolution with the NLO and LO nPDFs (cf. Fig. 3) and partly because the NLO computation probes slightly higher values of x than the LO case. Like in Fig. 9 Figure 10. The same as Fig. 9 but for neutral pions at a forward rapidity y = 3.

Centrality dependent R π 0
dAu (p T ) Let us then look at the centrality dependence of the ratio R π 0 dAu . Our NLO and LO results for R π 0 dAu in the centrality bins 0 − 20%, 20 − 40%, 40 − 60% and 60 − 88% are plotted in Figs. 11 and 12, correspondingly, together with the PHENIX data. Table 2 lists the impact parameter ranges and average number of binary collisions for each centrality class, obtained from the optical Glauber model with σ N N inel = 42 mb. Again, it is important to consider the different overall normalization errors in the experimental data. For the centrality-dependent ratios plotted in Figs. 11 and 12 there still is the 9.7 % overall systematic error due to the p+p baseline discussed above. In addition to this, an overall normalization error of 6.6-9.6 % arising from the determination of the average number of binary collisions, is quoted separately for each centrality bin. Following again the same procedure as for Fig. 9, we multiply the data and their point-to-point errors by a factor which minimizes the difference to our calculation. Even the largest upwards shift, 11.3 % for the centralmost bin in the LO case, is well within the acceptable total overall normalization error quoted by the experiment. Note also the systematic decrease of the multiplication factor from central to peripheral collisions, which we believe is due to the difference in the experimental and Glauber-model definitions of the centrality classes.   Figure 11. The nuclear modification factor R π 0 dAu (p T ) for √ s N N = 200 GeV at y = 0 for different centrality classes. Calculations are in NLO pQCD using EPS09sNLO1 and three different fragmentation functions. The blue error bands are computed with the error sets EPS09sNLOx (x=2,..., 31) and fDSS, and the data are from PHENIX [26]. The set-up and labeling are the same as in the left panel of Fig. 9. Notice that the experimental data have been multiplied by a different factor in each panel, which all are well within the total overall normalization uncertainties given by the experiment (6.6, 6.7, 8.5, 9.6 % for the four centrality bins from the Glauberization and 9.7 % from the p+p baseline.) From Figs. 11 and 12 we observe that within the experimental and theoretical uncertainties our calculations are consistent with the measurements. Especially the centrality systematics obtained from our spatially-dependent nPDFs agrees quite well with the data: the nuclear modifications are strongest in the most central collisions and systematically weaken when going to more peripheral collisions. This is especially nicely reflected in the region 1.3 ≤ p T ≤ 4 GeV, where the p T slopes (which are not affected by the overall multiplications) become steeper towards more central collisions. We also see that, like in the minimum-bias case, the EPS09 error bands are slighty smaller for the NLO than for the LO case, and that the uncertainties arising from the fragmentation functions remain small.
In Figs. 13 Figure 12. The same as in Fig. 11 but with LO calculations using EPS09sLO1 and EKS98s and for two different fragmentation functions, and the blue error band is computed with the EPS09sLOx (x=2,...,31) and fDSS. Notice again the overall multiplicative factors for the experimental data.
at a forward rapidity, y = 3, in the different centrality bins. In the forward region the nuclear modifications are larger since we now are probing smaller x values in the nPDFs than in the mid-rapidity region. Like in the minimum-bias case, we notice that the difference between the fragmentation function sets we use, becomes noticeable in the forward region. Again the small-p T suppression is stronger and nPDF-orginating uncertainties are larger for the LO case.

Predictions for p+Pb collisions at LHC
In the heavy-ion program of the LHC at CERN, there are now plans to collide protons with lead nuclei. Such collisions would be very useful for testing the QCD factorization and the universality of nPDFs, as well as for constraining the nuclear PDF modifications further especially at small values of x. Also the centrality dependence of nPDFs could be examined in these collisions via inclusive hadron production, similarly to the RHIC d+Au collisions discussed above but without the theoretical uncertainties arising from modeling the deuterium geometry. Thus, it is interesting to see what are the predictions from our spatially dependent nPDFs for these collisions. In Fig. 15 we plot our EPS09sNLO results for the nuclear modification factor R π 0 pPb (p T ) for neutral pion production in p+Pb collisions at 60-88% Figure 13. The nuclear modification factor R π 0 dAu (p T ) for √ s N N = 200 GeV at y = 3 in different centrality classes. The computation is done in NLO pQCD using EPS09sNLO1 and three different fragmentation functions. The error bands are computed with the EPS09sNLOx (x=2,...,31) and fDSS. centrality classes. 5 We use again the KKP, AKK and fDSS fragmentation functions here. The uncertainty bands arising from EPS09sNLO are computed using fDSS. The inelastic cross section σ N N inel = 70 mb for this √ s N N is obtained from Fig. 5 of Ref. [56]. This leads to the impact parameter values and the average number of binary collisions for each centrality class given in Table 3. For the projectile proton, we have not assumed any spatial size, so that relative to the deuterium case above, in the collision geometry we replace the thickness function T d (s) by δ(s) and the overlap function T dA (b) by the thickness function T Pb (b). As can be seen from Fig. 15, the nuclear modifications are strongest in the small-p T region in all centrality classes. To see the behaviour of R π 0 pPb in this region more clearly, we plot the results also in logarithmic scale in Fig. 16. We again observe the general behavior which follows from the spatial dependence of the nPDFs: the nuclear modifications are stronger in the central collisions and weaker in the peripheral collisions. We also notice that the three fragmentation function sets yield almost identical results. Figure 17 shows the corresponding ratio in minimum bias p+Pb collisions, computed both in NLO (left) and in LO (right). Like in the forward-rapidity case at RHIC, and for 5 Very recently, the LHC moved up to collisions energies √ spp = 8 TeV, hence we take √ sNN = √ spp Z/A ≈ 5.0 TeV. Note also that y is the rapidity in the N N cms frame, i.e. we do not include the rapidity shift due to the antisymmetric collision. 60-88% Figure 14. The same as Fig. 13 but for LO pQCD using EPS09sLO1 and EKS98s with two different fragmentation functions. The error bands are computed with the EPS09sLOx (x=2,...,31) and fDSS. the same reasons, the EPS09NLO leads to a weaker small-p T suppression than EPS09LO and EKS98, and the uncertainty band is clearly smaller for the NLO case. As a probe of nuclear gluons even deeper in the small-x shadowing region, we plot in Fig. 18 our LO results 6 for R π 0 pPb at a forward rapidity, y = 3, for the four centrality classes. Again the KKP and fDSS fragmentation functions are used, and we see that they yield very similar results. We should also point out that the EPS09sLO error band for the peripheral bin in Fig. 18 can be regarded as an underestimate in that it has been computed without 60-80% Figure 15. The nuclear modification factor R π 0 pPb (p T ) for √ s = 5.0 TeV at y = 0 for four different centrality classes, computed in NLO pQCD using EPS09sNLO1 and three different fragmentation functions. The error bands have been obtained with EPS09sNLOx (x=2,...,31) and fDSS. 60-80% Figure 16. The same as Fig. 15 but in a logarithmic scale to emphasize the small-p T region. the error set EPS09sLO7. The reason for this is that the error set EPS09LO7 gives in fact antishadowing at smallest x for the lightest nuclei, and in the EPS09sLO this maps into an antishadowing near the edges of a large nucleus. This unphysical feature can be cured only by redoing the EPS09LO global fit with an improved A-dependence of the fit functions. In the meantime, we suggest that a physically more meaningful upper limit for the LO error band in the small-x region for the peripheral bin can thus be obtained without this LO error set.
Finally, in Fig. 19 we show the minimum-bias R π 0 pPb at y = 3 for the NLO case at p T ≥ 5 GeV and for the LO case starting from p T = 1.3 GeV. Note the linear(logarithmic) p T scale on the left (right). Again we notice the weaker suppression and smaller error bands in the NLO case. Comparing the right panels of Figs. 19 and Fig. 17, we see that the smallest-p T suppressions are of similar magnitude. This is because the ratio R π 0 pPb in the small-p T region at the LHC probes already at y = 0 the flat part of the shadowing assumed as an input in EPS09 (cf. Fig. 4). Hence, a measurement of R π 0 pPb both in the midand forward-rapidities can be expected to serve as a relevant constraint for the smallest-x shadowing region.

Summary and Conclusions
We have developed a framework to determine the spatial dependence of the nuclear modifications of PDFs in such a way that the outcome is consistent with the globally analysed EKS98 and EPS09 nPDFs which in turn are DGLAP-based fits to nuclear hard-process data. Both the LO and NLO cases have been considered, and with EPS09 the spatial dependence has been extracted also for all the 30 error sets. Correspondingly, we call the obtained spatially dependent nPDF sets EPS09s and EKS98s.
The spatial dependence is introduced in terms of powers of the nuclear thickness func- 60-80%

tions T A (s). Regarding the power series r
shown that the 1-parameter approach (n = 1, used e.g. in [19,22]) is not sufficient for reproducing A systematics in the nPDFs, and that we obtain a good overall agreement with the globally analysed averaged nPDFs when we include terms up to [T A ] 4 . The outcome of the performed fits, the sets of coefficients {c i j (x, Q 2 )} for each parton flavor i at each x and Q 2 , are tabulated separately for each of the nPDF sets we considered. These tables along with a routine for interpolation and computing the needed thickness functions are downloadable at [34].
As a concrete application of our framework, we calculated the nuclear modification factor R 1jet AA for LO primary partonic jet production at different centralities in Au+Au collisions at RHIC and in Pb+Pb at the LHC. We observed that while the central R 1jet AA is quite close to the minimum-bias ratio R 1jet AA and the peripheral R 1jet AA differs fairly significantly from unity, the central-to-peripheral ratio R 1jet CP differs clearly from the ratio R 1jet AA . We also compared our NLO and LO calculations of the nuclear modification factor of neutral pion production in d+Au collisions, R π 0 dAu , in different centrality classes at midrapidity with the PHENIX data [26]. Within all the given errors in the experimental data, the nPDF uncertainties, and the possible differences between the experimental and optical Glauber model centrality classes, the EPS09s results are remarkably consistent with the centrality systematics. To our knowledge, this is the first time this has been demonstrated. Especially, our EPS09s results seem to reproduce the low p T slope of the data very well in all centrality classes.
More constraints for the spatial dependence of the nuclear PDFs, and gluons in particular, could be obtained from the scheduled p+Pb collisions at the LHC. We demonstrated this by calculating the NLO and LO predictions from our framework for the ratio R π 0 pPb in different centrality classes both at mid-rapidity y = 0 and forward rapidity y = 3 for √ s N N = 5.0 TeV, which corresponds to the recently achieved p+p cms-energy.
We believe that the nPDF development presented here is an important step forwards, as now a user may for the first time compute the centrality-dependent hard cross-sections more consistently with globally analysed nPDFs. Our spatially dependent nPDFs should also be applicable in Monte Carlo simulations of nuclear collisions, where the analogues of the thickness functions should be straightforwardly obtainable. In addition, our work should also give an idea how the future global analyses of nPDFs could be constructed so that the spatial dependence would be built in right from the start and not afterwards as has been the case here.
eling (see Refs. [59,60]), applied in this study. The calculations of T A (s) and T d (s) are included both in the EKS98s and EPS09s codes.

A.1 Nuclear Thickness Functions
The total amount of nuclear matter in a colliding nucleus A in the beam direction z at a transverse position s is given by the nuclear thickness function where ρ(s, z) is the nucleonic number-density of the nucleus, with a normalization convention In this study we use the standard two parameter Woods-Saxon density profile for ρ A , which is a good approximation for nuclei with A . (A.6)

A.1.2 Deuterium
For the thickness function of a deuterium nucleus, the above Woods-Saxon density profile is obviously not applicable anymore. Instead, one may formulate this with the deuteron wavefunction which describes the probability amplitude for the proton and neutron to be separated by a distance r pn . This can be written in terms of the 3 S 1 -and 3 D 1 -wave components as (see e.g. Ref. [48,61]) where the spin-spherical harmonics Y M JLS (Ω), with S = 1, consist of three components, For the radial parts we use the Hulthen form as in [48,58], where N 2 = 2α 1 − αρ , in which α −1 = 4.316 fm is related to the experimentally measured binding energy, and ρ is fixed by normalization, d 3 r pn |ψ M (r pn )| 2 = 1. For the other parameters, obtained by fitting to experimental data, we use the "set 1" quoted in [58]: The angular-averaged radial probability distribution for the proton-neutron distance r pn in deuteron is given by For computing the thickness function T d (s) as in Eq. (A.1), we need the nucleon density distribution ρ d (r) at a distance r from the center of mass of the deuteron. Assuming identical proton and neutron masses, we have r = r pn /2. In addition, we require the normalization of T d to be in line with Eq. (A.2). We thus have

A.2 Optical Glauber Model
Let us then specify the optical Glauber modeling applied for nuclear collisions in this study. For further discussion, see e.g. Refs. [59,60]. Consider first a nucleon-nucleus (N +A) collision at an impact parameter b. In the eikonal high collision-energy limit the number of binary inelastic collisions is given by where T A (b) is the thickness function defined in Eq. (A.1) and σ N N inel is the inelastic nucleonnucleon cross section. One may interpret N N A bin (b)/A as the probability for an inelastic collision to take place in the A N N collisions that are possible. Consequently, the probability for having no inelastic collisions at all, is 15) and the probabililty for at least one inelastic collision becomes The inelastic cross section for the N +A collision we then obtain as As the probability distribution above is expressable in terms of Poissonian probabilities for n inelastic collisions, |b| [fm] σ N N inel = 64 mb Figure 21. The differential inelastic cross section dσ AB inel /db as a function of impact parameter b for lead-lead (red solid), deuteron-lead (blue long-dashed), and proton-lead (green dashed) collisions with σ N N inel = 64 mb.
The centrality classes in the optical Glauber model can be defined as impact-parameter intervals. The "0-c % central" A+B collisions correspond to the most central collisions, 0 ≤ b ≤ b c which yield c % of the total inelastic cross section, The c 1 -c 2 % centrality class then corresponds to an interval [b 1 , b 2 ] for which For the studies of different hard-process nuclear modification factors in Sec. 4, we also define the average number of binary collisions in a c 1 -c 2 % centrality class: where the denominator is simply (c 2 −c 1 ) % of σ AB inel . For discussing the analogous centrality classes in N +A collisions, we just replace AB by N A in Eqs. (A.25-A.27) above, and also T AB by T A in Eq. (A.27).

B Nuclear modifications r uv and r us
For completeness, we plot here the nuclear modifications from our fits EPS09sNLO1, EPS09sLO1 and EKS98s in a lead nucleus for the u valence quarks in Fig. 22 and u sea quarks in Fig. 23. The corresponding modifications for gluons are shown in Fig. 5.