Impactparameter dependent nuclear parton distribution functions: EPS09s and EKS98s and their applications in nuclear hardprocesses
Abstract
We determine the spatial (impact parameter) dependence of nuclear parton distribution functions (nPDFs) using the 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 . To reproduce the dependence over the entire range we need terms up to . 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 hardprocess 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 for primary partonicjet production in different centrality classes in Au+Au collisions at RHIC and Pb+Pb collisions at LHC. Also the corresponding centraltoperipheral ratios are studied. We also calculate the LO and NLO nuclear modification factors for single inclusive neutral pion production, , 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 midrapidity data within the overall normalization uncertainties given by the experiment. Finally, we show our predictions for the corresponding modifications in the forthcoming p+Pb collisions at LHC.
a,b]Ilkka Helenius, a,b]Kari J. Eskola, a,c]Heli Honkanen d,e]Carlos A. Salgado \affiliation[a]Department of Physics, P.O. Box 35, FI40014 University of Jyväskylä, Finland \affiliation[b]Helsinki Institute of Physics, P.O. Box 64, FIN00014 University of Helsinki, Finland \affiliation[c]The Pennsylvania State University, 104 Davey Lab, University Park, PA 16802, USA \affiliation[d]Departamento de Física de Partículas and IGFAE, Universidade de Santiago de Compostela, GaliciaSpain \affiliation[e]Physics Department, Theory Unit, CERN, CH1211 Genève 23, Switzerland \emailAddilkka.helenius@jyu.fi \emailAddkari.eskola@phys.jyu.fi \emailAddhmh17@psu.edu \emailAddcarlos.salgado@usc.es \keywordsNuclear PDFs; hard processes, centrality dependence, nucleus+nucleus collisions, deuterium+nucleus collisions, proton+nucleus collisions \arxivnumberxxxx.yyyy
1 Introduction
In a highenergy hadronic or nuclear collision of particles and the inclusive cross sections for hard processes involving a large interaction scale can be computed using the QCD collinear factorization theorem [1, 2],
(1) 
where represents the perturbatively computable partonic pieces (cross sections in lowest order), and () is the parton distribution function (PDFs) for a given parton flavor in the colliding particle (and correspondingly for the flavor in ). The PDFs are universal, processindependent functions of nonperturbative origin, whose evolution in the scale 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 BNLRHIC and CERNLHC. This holds as well for protonproton collisions as for protonnucleus and nucleusnucleus collisions. To determine the nonperturbative input in the PDFs, one has developed global analyses which exploit a multitude of experimental hardprocess 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 freenucleon PDFs. Analogously to the freeproton case, global DGLAP analyses have been developed also for the nPDFs. The further complication with these is that in addition to the usual and dependences also the nuclear massnumber () dependence of the PDFs needs to be dealt with. The global nPDF fits have so far resulted in leadingorder (LO) nPDF sets EKS98 [10], HKM [11] and HKN04 [12], and nexttoleading 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 hardprocess crosssections.
The global analyses mentioned above have all considered only the spatially averaged nPDFs, probed in minimumbias nuclear collisions with no cuts on the collision centrality (impact parameter). In particular, as the modest amount of available nuclear hardprocess data severely limits the number of possible fit parameters, it has so far been impossible to embed the spatial dependence, or the impactparameter 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 hardprocess crosssections 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 GribovGlauber 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 selfconsistent 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 production [31, 32] in different centrality classes in d+Au collisions at RHIC, as well as by the hardprocess measurements in the forthcoming p+Pb collisions at the LHC. Also the theoretical modeling of the 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 in each nucleus at each and in terms of a power series of the standard nuclear thickness functions . Then, we determine the coefficients of each power of by exploiting the 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 hardprocess crosssections in different centrality bins, we first discuss the centrality dependence of the LO nuclear modification ratios of primary partonicjet production in Au+Au collisions at RHIC and Pb+Pb collisions at the LHC. We also study the nuclear modification factors of inclusive production, , in d+Au collisions at RHIC and in p+Pb collisions at the LHC, both at mid and forward rapidities, and considering both the NLO and LO cases. For 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 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.
2 The Analysis Framework
2.1 Definitions of the Nuclear Modifications
First we need to define how we introduce the spatial dependence to the nPDFs in terms of the hardprocess crosssections. Let us start with the usual spatially averaged nPDFs. The number distribution of an observable produced in a collision of nuclei and at an impact parameter b is given by
(2) 
where is the standard nuclear overlap function normalized to (cf. Eq. (40) in App. A.2, see the nuclear collision geometry in Fig. 20), and is the bindependent inclusive hard crosssection of Eq. (1) containing the nPDFs and perturbative pieces. The spatially averaged nPDFs in a nucleus with protons and neutrons are now given by
(3) 
where the nPDFs of a bound neutron, , may be (approximately) obtained from those of the bound proton, , 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 and the corresponding free proton PDF ,
(4) 
To lighten the notations, we express the nPDFs in Eq. (3) as
(5) 
where the sum runs over all the nucleons .
Decomposing the into the standard nuclear thickness functions (cf. Eq. (40)), and using Eq. (1) we may write
(6) 
From this, we see that a suitable definition of the spatially dependent nuclear modification for the PDF of parton flavor (per nucleon) is
(7) 
where the thickness function is normalized to and where the case of no nuclear effects corresponds to . Using these definitions, we can now generalize Eq. (6) to include the spatially dependent nuclear modifications,
(8) 
As a consistency check, we note that the definition in Eq. (7) guarantees that the minimumbias cross sections, which are obtained by integrating Eq. (8) over the whole b space, become simply times the hard crosssection computed with the spatially averaged nPDFs,
(9) 
The key assumption in the present analysis is that the spatial dependence of is a function of the nuclear thickness . The motivation for this comes mainly from the shadowing region at small , where the partons of sufficiently small values of may interact with partons from any other nucleon near enough in the transverse direction. Also in the GribovGlauber modeling [23] of the initial state nPDFs the nuclear effects become essentially functions of . The functional form we choose to use and test here is a simple power series of the thickness functions,
(10) 
Here we would like to emphasize the following points: First, all the dependence is now in the thickness functions which are fully known, and all the coefficients which will be our fit parameters, depend on and but not on . Second, the power series of the form also fixes by construction that when , 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 systematics, i.e. weighted integrals of Eq. (10) over for each nucleus. Fourth, since the EKS98 and EPS09 global analyses have not been constructed to reproduce any specific theoretically motivated dependence of the nPDFs, we can test the validity of the assumption of Eq. (10), as well as the number of terms needed, only a posteriori.
Using the definitions above, we can see why the simplest 1parameter approach with in Eq. (10) (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 systematics of the nuclear data. In this case, , and from the definition in Eq. (7), one obtains , where is given by the globally analysed nPDFs (i.e. nuclear data). The problem then is that the coefficient may depend in fact quite strongly on , 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 dependence is illustrated in Fig. 2 for gluons in a lead nucleus at 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 systematics in the spatially dependent nPDFs.
2.2 Fitting Procedure
To extract the independent coefficients , we need to introduce a fitting procedure, where we utilize the definition (7) and the dependence of the EKS98 and EPS09 nuclear modifications at different values of and for each parton flavor . To reproduce the systematics in the spatially independent nuclear modifications with the powerseries ansatz of Eq. (10), we minimize the defined as
(11) 
where the spatially averaged modifications from EKS98 and EPS09 now represent the ”experimental” data. The weight factors 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 4thorder polynomials in , i.e. in Eq. (10). Furthermore, best fits were obtained with the weight for EKS98 (this corresponds to fitting the deviations from unity within a constant relative error) and 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 , and we do not expect the fit form of Eq. (10) work for the smallest values of , either. Consequently, we exclude the nuclei from the fit. Thus for EKS98 the sum runs over (i.e. emphasising the large nuclei) and for EPS09 we use all the values for which these sets are currently available.
3 Results
3.1 Quality of the fit
First, we demonstrate that our fit framework manages to reproduce the spatially averaged nuclear modifications and especially their dependence indeed very well. Figure 3 shows the obtained spatially dependent gluon modifications integrated over the transverse plane according to Eq. (7), and the corresponding input modifications at different fixed values of 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 for EPS09NLO, 1 % at for EPS09LO, and 0.2 % at for EKS98. Importantly, the keyfeature here, the 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 centralitydependent nuclear hard crosssections can now be done as before, using the EPS09s sets.
3.2 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 at the initial scale as a function of and , 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 shape of the nuclear modification away from the edge of the nucleus, at , is similar as in the input distribution. This confirms that our fit does not generate any unwanted extra curvature.

The nuclear modification dies out as expected, by construction, when . This feature arises from the vanishing 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.
3.3 Comparison with other approaches
Next, we compare our EPS09s and EKS98s fits with the 1parameter approach described in the end of Sec. 2.1. [20, 22].
In Fig. 6 we plot the nuclear modification for gluons at fixed values of and scale for as a function of from our EPS09sNLO1, EPS09sLO1 and EKS98s fits, from the 1parameter 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 1parameter approach leads to a too steep transverse profile for the modifications in all cases.
4 Applications
Next, we consider some concrete examples of computing the nuclear hardprocess crosssections in different centrality classes using the spatially dependent nPDFs. First, we discuss the centrality dependence of primary partonicjet production in + 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.
4.1 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 impactparameter interval . The generic average number distribution of a hardprocess observable in this centrality class of an + collision is
(12) 
where from Eq. (43), and is obtained from Eq. (8). Using the expansion of in powers of from Eq. (10), 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:
(13) 
where are the free nucleon PDFs, and we have defined and
(14) 
From Eq. (13) we see the most straightforward implementation of the spatially dependent nPDFs. The purely geometric integrals, the in Eq. (14) for each pair of the powers and , can be computed independently of the kinematic variables , , and and also independently of the parton flavors . Thus, in total we have 25 different geometric integrals to do (or 15 if ) but we need to do them only once. In comparison with the spatially averaged case, the fit parameters thus play the role of the nuclear modifications for each of the 25 pairs . To arrive at the final bintegrated result for the number distribution of , we thus need to repeat the computation of the kinematic parts 25 times, each with different sets of the coefficient pairs and a different geometric weight . The EKS98s
and EPS09s
routines which we provide in [34], give in addition to the fit coefficients also the thickness functions (used in the fits here) for the computation of , as well as the combination for other possible implementations.
Note also that for the b integral in Eq. (14) the angular part is trivial, giving just .
4.2 The Nuclear Modification Factors and
Let us now consider the centrality dependence of primary inclusive high parton production in + collisions at RHIC and LHC. Following the generic discussion above, we define the nuclear modification ratio relative to the p+p case for each centrality class as
(15) 
where is the average number of binary collisions in this centrality class given by Eq. (47) and is the inelastic nucleonnucleon 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 minimumbias collisions is obtained from above by setting and , in which case we have
(16) 
where , which contains only the spatially averaged nPDFs, is obtained from Eq. (9) by setting , and the p+p baseline from the same equation by setting .
In addition to the centrality dependence of , we are interested in the centraltoperipheral ratios, defined as
(17) 
where and refer to the central and peripheral bins, correspondingly. The advantage of this ratio (in the experiments) is that the information of the protonproton baseline is not required. In particular, we would like to see exactly how much differs from the modification which is computed with the spatially averaged nPDFs. We perform these examplecalculations 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 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 05% and peripheral 6080% 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.
Central = 05 %  Peripheral = 6080 %  

[GeV]  [mb]  
Au+Au  200  42  0.0  3.355  1083  11.62  13.42  15.10 
Pb+Pb  2760  64  0.0  3.478  1771  12.05  13.91  19.08 
In Fig. 7 we plot the ratio for central, peripheral and minimumbias collisions, as well as in Au+Au collisions at RHIC. Figure 8 shows the same quantities for the LHC Pb+Pb case. The central and peripheral and have been obtained with the spatially dependent nPDFs EPS09sLO1 (left) and EKS98s (right), and the average 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 has been set to be the transverse momentum, , of the parton.
The main observations from the figures are: (i) The central is quite close to the average , which is expected since the nuclear modifications at small are close to the average modifications. (ii) The peripheral is clearly not unity but there appear almost 10% antishadowing effects at mid at RHIC and even more than 20% shadowing effects at small at the LHC, and up to 10% EMC effects at large both at RHIC and LHC. (iii) Consequently, the ratio differs significantly from the average . The results suggest that in a precision theoryanalysis 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 slopes of at small than in the EKS98 case.
4.3 Centrality dependence of at RHIC – comparison with data
While the above ratios and mainly serve as theoretical pQCD baselines for jet quenching studies, it is important to test our spatiallydependent nPDF framework against some measured centralitydependent observables. To avoid the complications of hot QCD matter modeling, we turn to the highestenergy d+Au collisions at RHIC and p+Pb at the LHC. For our purposes a promising published data set is the nuclear modification factor for single inclusive neutralpion production at midrapidity , measured by PHENIX [26] at different centrality classes in d+Au collisions at GeV. Since the minimumbias 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. (15), we define the centralitydependent nuclear modification factors as
(18) 
where the number distribution now involves a further folding over the fragmentation functions,
(19) 
and where is obtained from Eq. (8) by setting and , and (as we do not assign any nuclear effects for the deuterium PDFs) also . 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.
Minimumbias
Setting the spatial integrals in Eq. (18) over the whole impactparameter space gives again the minimumbias ratios,
(20) 
where again contains only the spatially averaged nPDFs in which is obtained from Eq. (9) by setting , . As noted earlier, in the EKS98 and EPS09 frameworks there are no nuclear modifications to the deuterium PDFs. The p+p baseline 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 in minimumbias 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]INCNLO
code [52, 53].
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 region there is a difference in the 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 pointtopoint 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 with the pointtopoint statistical and systematic errors added in quadrature), we have multiplied the data by a factor 1.039 (NLO) and 1.050 (LO). Such a fewpercent 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 at the forward region in minimumbias collisions. Again, we show both the NLO (left) and LO (right) results with the same setup as in Fig. 9 above. Now the differences between the fragmentation functions start to be visible, as one is probing their larger tails where the uncertainties are larger: for the same pion , the differences in the large fragmentation functions map to different values of in the nPDFs.
We also notice that the LO calculation gives a stronger small 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 than the LO case. Like in Fig. 9, the EPS09 error band is smaller for NLO than for LO.
Centrality dependent
Let us then look at the centrality dependence of the ratio . Our NLO and LO results for in the centrality bins , , and 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 .
0.0  3.798  15.57  
3.798  5.371  10.95  
5.371  6.583  6.013  
6.583  8.336  2.353 
Again, it is important to consider the different overall normalization errors in the experimental data. For the centralitydependent 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 pointtopoint 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 Glaubermodel definitions of the centrality classes.
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 spatiallydependent 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 GeV, where the slopes (which are not affected by the overall multiplications) become steeper towards more central collisions. We also see that, like in the minimumbias 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 and 14 we plot also our NLO and LO results for at at a forward rapidity, , in the different centrality bins. In the forward region the nuclear modifications are larger since we now are probing smaller values in the nPDFs than in the midrapidity region. Like in the minimumbias case, we notice that the difference between the fragmentation function sets we use, becomes noticeable in the forward region. Again the small suppression is stronger and nPDForginating uncertainties are larger for the LO case.
4.4 Predictions for p+Pb collisions at LHC
In the heavyion 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 . 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. 16 we plot our EPS09sNLO results for the nuclear modification factor for neutral pion production in p+Pb collisions at at
in four different centrality classes.
0.0  3.471  14.24  
3.471  4.908  11.41  
4.908  6.012  7.663  
6.012  6.986  3.680 
As can be seen from Fig. 16, the nuclear modifications are strongest in the small region in all centrality classes. To see the behaviour of 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 forwardrapidity case at RHIC, and for the same reasons, the EPS09NLO leads to a weaker small 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 shadowing region, we plot in Fig. 18 our LO results
Finally, in Fig. 19 we show the minimumbias at for the NLO case at GeV and for the LO case starting from GeV. Note the linear(logarithmic) 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 suppressions are of similar magnitude. This is because the ratio in the small region at the LHC probes already at the flat part of the shadowing assumed as an input in EPS09 (cf. Fig. 4). Hence, a measurement of both in the mid and forwardrapidities can be expected to serve as a relevant constraint for the smallest shadowing region.
5 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 DGLAPbased fits to nuclear hardprocess 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 functions . Regarding the power series , we have shown that the 1parameter approach (, used e.g. in [19, 22]) is not sufficient for reproducing systematics in the nPDFs, and that we obtain a good overall agreement with the globally analysed averaged nPDFs when we include terms up to . The outcome of the performed fits, the sets of coefficients for each parton flavor at each and , 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 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 is quite close to the minimumbias ratio and the peripheral differs fairly significantly from unity, the centraltoperipheral ratio differs clearly from the ratio .
We also compared our NLO and LO calculations of the nuclear modification factor of neutral pion production in d+Au collisions, , 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 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 in different centrality classes both at midrapidity and forward rapidity for , which corresponds to the recently achieved p+p cmsenergy.
We believe that the nPDF development presented here is an important step forwards, as now a user may for the first time compute the centralitydependent hard crosssections 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.
Acknowledgements
We thank M. Wysocki, J. Rak, T. Lappi, T. Renk and H. Mäntysaari for discussions. We gratefully acknowledge the following financial support: I.H. from the Magnus Ehrnrooth Foundation; K.J.E., I.H. and H.H. from the Academy of Finland, K.J.E.’s Project No. 133005. C.A.S. is supported by the European Research Council grant HotLHC ERC 2001StG279579 and by Ministerio de Ciencia e Innovación of Spain. C.A.S. is a Ramón y Cajal researcher. This work (H.H.) was supported in part by the U.S. Department of Energy under Grant DE FG0293ER40771.
Appendix A Nuclear Collision Geometry
For clarity, we specify here the modeling and parameters of the nuclear collision geometry, i.e. the nuclear thickness functions and (see Refs. [57, 58]) and Glauber modeling (see Refs. [59, 60]), applied in this study. The calculations of and are included both in the EKS98s
and EPS09s
codes.
a.1 Nuclear Thickness Functions
Large nuclei
The total amount of nuclear matter in a colliding nucleus in the beam direction at a transverse position is given by the nuclear thickness function
(21) 
where is the nucleonic numberdensity of the nucleus, with a normalization convention
(22) 
In this study we use the standard two parameter WoodsSaxon density profile for ,
(23) 
which is a good approximation for nuclei with . The parameter values for the WoodsSaxon distribution are
(24)  
(25) 
and for large nuclei the normalization condition (22) fixes the constant as
(26) 
Deuterium
For the thickness function of a deuterium nucleus, the above WoodsSaxon 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 . This can be written in terms of the  and wave components as (see e.g. Ref. [48, 61])
(27) 
where the spinspherical harmonics , with , consist of three components,
(28) 
For the radial parts we use the Hulthen form as in [48, 58],
(29)  
(30)  