A Nuclear Collision Geometry

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

## 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 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 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 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 mid-rapidity 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, FI-40014 University of Jyväskylä, Finland \affiliation[b]Helsinki Institute of Physics, P.O. Box 64, FIN-00014 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, Galicia-Spain \affiliation[e]Physics Department, Theory Unit, CERN, CH-1211 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 high-energy 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],

 dσAB→k+X=∑i,j,X′fAi(Q2)⊗fBj(Q2)⊗d^σij→k+X′+O(1/Q2), (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, process-independent 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 hard-process 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 and dependences also the nuclear mass-number () 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 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 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 hard-process cross-sections in different centrality bins, we first discuss the centrality dependence of the LO nuclear modification ratios 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 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 hard-process cross-sections. 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

 dNAB→k+X(b)=TAB(b)dσAB→k+X, (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 b-independent inclusive hard cross-section of Eq. (1) containing the nPDFs and perturbative pieces. The spatially averaged nPDFs in a nucleus with protons and neutrons are now given by

 fAi(x,Q2)=ZAfp/Ai(x,Q2)+A−ZAfn/Ai(x,Q2), (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 ,

 fp/Ai(x,Q2)≡RAi(x,Q2)fpi(x,Q2). (4)

To lighten the notations, we express the nPDFs in Eq. (3) as

 fAi(x,Q2)=1A∑NRN/Ai(x,Q2)fNi(x,Q2), (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

 dNAB→k+X(b)=∑i,j,X′1AB∑NA,NB∫d2s1TA(s1)RNA/Ai(x1,Q2)fNAi(x1,Q2)⊗∫d2s2TB(s2)RNB/Bj(x2,Q2)fNBj(x2,Q2)⊗d^σij→k+X′δ(s2−s1−b). (6)

From this, we see that a suitable definition of the spatially dependent nuclear modification for the PDF of parton flavor (per nucleon) is

 Extra open brace or missing close brace (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,

 dNAB→k+X(b)=∑i,j,X′1AB∑NA,NB∫d2s1TA(s1)rAi(x1,Q2,s1)fNAi(x1,Q2)⊗∫d2s2TB(s2)rBj(x2,Q2,s2)fNBj(x2,Q2)⊗d^σij→k+X′δ(s2−s1−b). (8)

As a consistency check, we note that the definition in Eq. (7) guarantees that the minimum-bias cross sections, which are obtained by integrating Eq. (8) over the whole b space, become simply times the hard cross-section computed with the spatially averaged nPDFs,

 dσAB→k+XMB=∫d2bdNAB→k+X(b)=AB∑i,j,X′fAi(x,Q2)⊗fBj(x,Q2)⊗d^σij→k+X′. (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 Gribov-Glauber 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,

 rAi(x,Q2,s)=1+n∑j=1cij(x,Q2)[TA(s)]j. (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 1-parameter 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 power-series ansatz of Eq. (10), we minimize the defined as

 χ2i(x,Q2)≡∑A⎡⎣RAi(x,Q2)−1A∫d2sTA(s)rAi(x,Q2,s)WAi(x,Q2)⎤⎦2, (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 4th-order 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 key-feature 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 centrality-dependent nuclear hard cross-sections 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.

• In the center of the nucleus, , 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 . 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 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 suppression e.g. in Refs. [42, 31, 43, 33]2 and inclusive hadron production in d+Au collisions at RHIC in Ref. [22]. We also compare with the leading-twist formulation [44, 23] 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- 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 [45, 23], and choose the value of 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 and scale for as a function of 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.

## 4 Applications

Next, we consider some concrete examples of computing the nuclear hard-process cross-sections in different centrality classes using the spatially dependent nPDFs. First, we discuss the centrality dependence of primary partonic-jet 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 impact-parameter interval . The generic average number distribution of a hard-process observable in this centrality class of an + collision is

 ⟨dNkAB⟩b1,b2=∫b2b1d2bdNkAB(b)∫b2b1d2bpinelAB(b), (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:

 ∫b2b1d2bdNkAB(b)= 4∑n,m=0TnmAB(b1,b2)∑i,j,X′1AB∑NA,NBcin(x1,Q2)fNAi(x1,Q2)⊗ cjm(x2,Q2)fNBj(x2,Q2)⊗d^σij→k+X′ (13)

where are the free nucleon PDFs, and we have defined and

 TnmAB(b1,b2)≡∫b2b1d2b∫d2s[TA(s−b/2)]n+1[TB(s+b/2)]m+1. (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 b-integrated 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 R1jetAA and R1jetCP

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

where is the average number of binary collisions in this centrality class given by Eq. (47) and 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 and , in which case we have

 Extra open brace or missing close brace (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 central-to-peripheral ratios, defined as

where and 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 differs from the modification 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- 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 for central, peripheral and minimum-bias 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 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 slopes of at small than in the EKS98 case.

### 4.3 Centrality dependence of Rπ0dAu(pT) 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 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 for single inclusive neutral-pion production at mid-rapidity , measured by PHENIX [26] at different centrality classes in d+Au collisions at  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. (15), we define the centrality-dependent nuclear modification factors as

where the number distribution now involves a further folding over the fragmentation functions,

 dNπ0dA(b)=∑kdNkdA(b)⊗Dπ0/k(z,Q2F) (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.

#### Minimum-bias Rπ0dAu(pT)

Setting the spatial integrals in Eq. (18) over the whole impact-parameter space gives again the minimum-bias ratios,

 Extra open brace or missing close brace (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 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 and fragmentation scale are all fixed to , the transverse momentum of the produced hadron. For details of the LO calculation, we again refer to Ref. [46], while the NLO computation was performed by using the 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 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 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 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- 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.4

#### Centrality dependent Rπ0dAu(pT)

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 .

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.

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  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 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 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 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- suppression is stronger and nPDF-orginating uncertainties are larger for the LO case.

### 4.4 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 . 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.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 for this 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 by and the overlap function by the thickness function .

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 forward-rapidity 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 results6 for at a forward rapidity, , 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 the error set EPS09sLO7. The reason for this is that the error set EPS09LO7 gives in fact antishadowing at smallest 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 -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- region for the peripheral bin can thus be obtained without this LO error set.

Finally, in Fig. 19 we show the minimum-bias 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 forward-rapidities 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 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 functions . Regarding the power series , we have shown that the 1-parameter 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 minimum-bias ratio and the peripheral differs fairly significantly from unity, the central-to-peripheral 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 mid-rapidity 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 mid-rapidity and forward rapidity for , 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.

## 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- 2001-StG-279579 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- FG02-93ER40771.

## 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

 TA(s)=∫∞−∞dzρA(s,z), (21)

where is the nucleonic number-density of the nucleus, with a normalization convention

 A=∫d2sTA(s). (22)

In this study we use the standard two parameter Woods-Saxon density profile for ,

which is a good approximation for nuclei with . The parameter values for the Woods-Saxon distribution are

 d = 0.54 fm (24) RA = 1.12A1/3−0.86A−1/3 fm, (25)

and for large nuclei the normalization condition (22) fixes the constant as

 n0=34AπR3A1(1+(πdRA)2). (26)

#### 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 . This can be written in terms of the - and -wave components as (see e.g. Ref. [48, 61])

 ψM(rpn)=u(rpn)rpnYM101(Ω)+w(rpn)rpnYM121(Ω), (27)

where the spin-spherical harmonics , with , consist of three components,

 [YM101(Ω)]mS=±1,0=⟨Ω,mS|LSJM⟩=∑ML,MS⟨LSMLMS|LSJM⟩YLML(Ω)δmSMS. (28)

For the radial parts we use the Hulthen form as in [48, 58],

 u(rpn)= N√1−ϵ2[1−e−β(αrpn−xc)]eαrpnθ(αrpn−xc) (29) w(rpn)= Missing or unrecognized delimiter for \right (30) [1+3(1−e−γαrpn)αrpn+