{'fraction_non_alphanumeric': 0.03753823436346377, 'fraction_numerical': 0.015040607530851176, 'mean_word_length': 4.455492260774498, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Extraction of Spin-Dependent Parton Densities and Their Uncertainties

We discuss techniques and results for the extraction of the nucleon's spin-dependent parton distributions and their uncertainties from data for polarized deep-inelastic lepton-nucleon and proton-proton scattering by means of a global QCD analysis. Computational methods are described that significantly increase the speed of the required calculations to a level that allows to perform the full analysis consistently at next-to-leading order accuracy. We examine how the various data sets help to constrain different aspects of the quark, anti-quark, and gluon helicity distributions. Uncertainty estimates are performed using both the Lagrange multiplier and the Hessian approaches. We use the extracted parton distribution functions and their estimated uncertainties to predict spin asymmetries for high-transverse momentum pion and jet production in polarized proton-proton collisions at 500 GeV center-of-mass system energy at BNL-RHIC, as well as for W boson production.


I. INTRODUCTION
The last twenty years have witnessed remarkable improvements in the sophistication and precision of methods devoted to the extraction of parton density and fragmentation functions from experimental hard scattering data. These distributions are essential ingredients for any phenomenological study of hard scattering processes involving identified hadrons in the initial and final-state, respectively. Their precise knowledge is not only critical for testing the very successful framework of perturbative Quantum Chromodynamics (pQCD) but, in more general terms, also for quantifying uncertainties for precision studies of the Standard Model and searches of "new physics" at high energy accelerators like the CERN-LHC. At the same time, parton densities and fragmentation functions give fundamental insights into nucleon structure and the hadronization mechanism.
With the gain in precision and refinements of analyses, modern parton distribution functions (PDFs) have often revealed intriguing aspects of hadronic structure, such as the sizable breaking of isospin symmetry in the light sea quark sector, suggestions of differences between the strange quark and anti-quark distributions, the steep rise of the distributions at small momentum fractions, and an interesting pattern of modifications of the distributions in nuclei, to name just a few. Certainly one of the most striking results is the unexpectedly small fraction, about a quarter, of the proton's spin that can be attributed to the intrinsic angular momenta of quarks and anti-quarks. This finding, famously dubbed "proton spin crisis", has triggered a flurry of experimental and theoretical activity aiming at clarifying the contributions of gluons and orbital angular momenta of partons to the spin of the proton [1].
The only way to effectively deconvolute the experimental information on PDFs, which in its raw form is smeared over the light-cone momentum fraction x, summed over many different partonic subprocesses, and taken at different hard scales Q for each data point, is a "global QCD analysis". It treats all available probes simultaneously, in order to extract the set of universal PDFs that yields the optimal theoretical description of the combined data. For the case of polarized PDFs, the available world-data are from polarized deep-inelastic scattering (DIS) [2,3,4,5,6,7,8,9,10,11,12,13], semiinclusive DIS (SIDIS) [10,14,15,16], photo-and electroproduction of hadrons and charm [17,18,19,20,21], and proton-proton (pp) collisions at BNL-RHIC [22,23,24,25,26,27]. The different data sets are complementary in the sense that they probe different aspects of the helicity dependent PDFs. Fully inclusive DIS data from the many different experiments are pivotal in precisely determining the sums of quark and anti-quark distributions, SIDIS data help to tell different quark flavors and quark and anti-quarks apart, and RHIC pp data give a first direct constraint on gluon polarization.
A global QCD analysis of nucleon spin structure at full next-to-leading (NLO) accuracy was completed recently [28]. The present paper gives in large part a more detailed account of the methods and results of [28]. It also addresses the issue of the uncertainties of the PDFs in a more detailed and comprehensive way. As customary for recent unpolarized PDF analyses [29,30], we provide sets of polarized PDFs associated with displacements in the PDF parameter space in the vicinity of the best fit which greatly facilitate the propagation of PDF uncertainties to any observable of interest.
As a new feature over all previous fits based only on DIS [31,32,33,34,35], or combined DIS and SIDIS data [36], the analysis [28] included for the first time also results from polarized pp scattering at RHIC in a NLO framework. It benefitted significantly from an improved knowledge of parton-to-hadron fragmentation functions [37] which are an essential non-perturbative input for the theoretical description of all processes with identified hadrons in the final state, such as SIDIS. For the first time, these fragmentation functions provide a good description of identified hadron yields in the entire kinematic regime relevant for the analysis of polarized SIDIS and pp data [37]. For the time being, hadron and charm production data from fixed-target experiments [17,18,19,20,21], which constrain the gluon polarization at momentum fractions around x ≃ 0.1, were not included in the analysis [28] since NLO calculations of the relevant cross sections are not yet complete [38]. We will provide a comparison to these data in order to demonstrate their consistency with the results of the global fit. RHIC data for charged pion spin asymmetries [27] are also not taken into account as they are still preliminary and statistically not as significant as the neutral pion or jet data [22,23,24,25,26]. With sufficient statistics, however, they can provide an important constraint on the gluon polarization as will be shown below.
The use of parton distributions in predictions for Standard Model benchmark processes, e.g., as "luminosity candles" at the LHC, or in understanding fundamental properties of a nucleon like its spin, not only requires a careful extraction of PDFs from data but also a proper assessment of their uncertainties and how they propagate to other observables of interest. In spite of a great deal of activity and many significant achievements for both parton distribution and fragmentation functions (FFs), this has shown to be a rather formidable task in practice [29,30,37,39,40]. The specific challenge of a global QCD analysis is to incorporate a large body of data from many experiments with diverse characteristics and errors. The complications are compounded by uncertainties inherent to the theoretical framework used to describe the data, which are notoriously difficult to quantify. Examples are the choice of the factorization scale, the functional form used to parameterize the PDFs, or unavoidable approximations and assumptions limiting the parameter space.
Several complementary strategies have been devised and implemented to estimate uncertainties of PDFs and FFs [41,42,43]. In general, one starts with introducing an effective χ 2 function that combines all phenomenological inputs to the analysis as a quantitative measure of the goodness of the global fit. Minimizing this χ 2 function yields the optimal set or "best fit" of parameters in the multi-dimensional space defining the PDFs. The most common method to determine the range of uncertainties is to study the dependence of χ 2 near its global minimum based on a Taylor expansion and keeping only the leading term as characterized by the error matrix or its inverse, the Hessian matrix [43]. This assumes a quadratic form in the displacements of all parameters from their optimum values. The Hessian also determines the uncertainties of any other physical observable O, provided that the dependence of O on the fit parameters is approximately linear around the minimum. Both assumptions are not necessarily adequate in the complex global analysis environment, and their range of applicability needs to be carefully scrutinized.
The more robust method of Lagrange multipliers [42] circumvents all these shortcomings and is free of assumptions concerning the functional dependence of χ 2 on the fit parameters. The idea is to explore directly how the fit to data deteriorates if one enforces certain values for an observable O away from its best fit value. In practice, one performs a series of constrained fits in which χ 2 is minimized for particular values of O, in order to map out the parametric relationship between χ 2 and O. The method is straightforward to implement and can be applied to any combination of physical observables or even to fit parameters themselves. We will pursue and compare both methods, Hessian and Lagrange multiplier, to estimate uncertainties for the shape and truncated first moments of helicity-dependent PDFs using the analysis presented in [28] as the starting point. The Lagrange multiplier technique will provide the necessary benchmarks for testing the accuracy of approximations within the Hessian method. We note that alternative approaches recently proposed in the literature include studies of uncertainties based on neural networks or large samples of PDFs generated with Monte-Carlo methods [40].
In any case, all of these methods require an extensive number of calculations and minimizations of the effective χ 2 function, in order to explore the very complex and entangled sensitivity of the data to variations of the parameters describing the PDFs. This calls for new and more efficient calculational tools to include all observables used in the global fit consistently at NLO accuracy without resorting to potentially unreliable approximations. In particular, numerical computations of NLO cross sections in hadron-hadron scattering are prone to being very time consuming.
In order to deal with this problem, Ref. [28] employed a method based on the Mellin transform technique proposed in Refs. [44,45], which allows to speed up the relevant NLO computations to a level that they can be incorporated in the global analysis. An important improve-ment of this Mellin transform method was the implementation of a novel Monte-Carlo sampling technique. This new computational strategy proved already very useful for the analysis of single-inclusive observables in pp scattering that are presently relevant at RHIC. However, it is completely general and becomes especially powerful when less inclusive observables like two-particle correlations in hadronic collisions need to be incorporated in the global QCD fits, as will soon be the case. In the present paper we will describe the Mellin transform technique and its improvement in detail. We note, that the fast and efficient Mellin technique for incorporating NLO pp processes is, of course, not restricted to analyzing helicitydependent PDFs, but could equally find important applications for QCD processes at the LHC.
The paper is organized as follows: in the next Section we describe all technical details of a global PDF analysis. We first discuss the χ 2 function and the underlying ideas of the Hessian and Lagrange multiplier approaches for estimating PDF uncertainties. We next describe in detail our Mellin moment and Monte-Carlo sampling techniques as implemented for fast evaluations of NLO pp cross sections in our global QCD analysis. In Sec. III we apply all techniques to the global analysis of helicity-dependent PDFs [28]. We discuss the results for the best fit and its uncertainties. We argue that the range of applicability of the Hessian method is limited to estimating uncertainties of helicity-dependent PDFs consistent with only small departures from the best global fit, corresponding to ∆χ 2 ≈ 1. We present sets of polarized PDFs associated with displacements along the eigenvector directions of the Hessian matrix and resulting in ∆χ 2 = 1, which characterize the PDF parameter space in the vicinity of the global minimum in a processindependent way. We also explore the impact of the individual data sets on the results and uncertainties obtained for helicity-dependent PDFs in the global fit. In Sec. IV we study the potential of upcoming measurements at RHIC at √ S = 500 GeV center-of-mass system (c.m.s.) energy for further constraining the polarized PDFs. We focus on predictions for single-inclusive pion and jet production, and on W boson single-spin asymmetries. We conclude in Sec. V.

II. TECHNIQUES FOR NLO GLOBAL PDF ANALYSES
In this Section, we will describe all techniques we use for the global analysis of polarized PDFs. The first two Subsections discuss the χ 2 function and the various methods for the analysis of PDF uncertainties. Much of the discussion here will follow the pioneering work in Refs. [41,42,43]. We then lay out the details of our Mellin moment and Monte-Carlo sampling techniques.
A. The effective χ 2 function Global QCD extractions of PDFs [28,29,30,39,46,47] or FFs [37] are implemented around an effective χ 2 function that quantifies the goodness of the fit to data for a given set of theoretical parameters {a i }, i = 1, . . . , N par that determines the PDFs or FFs at some initial scale µ 0 . The simplest χ 2 function, convenient for the search for optimum PDFs by minimization, is usually taken as where N exp counts the individual experimental data sets and N (n) data the corresponding number of data points in each set. Each data value D j is compared to the corresponding theoretical estimate T j , which depends in general non-linearly on the N par parameters {a i }, weighted with the estimated uncertainties combined in δD j . In Eq. (1) ω j is a special weighting factor for each data point with default value one. It can be set to zero if a certain data point is to be removed from the analysis due to some physics considerations. For instance, such cuts are routinely introduced in a global fit to remove kinematical regions where the framework of perturbative QCD used to compute T j ({a i }) is known to be not adequate for describing the available data. The simple form (1) for χ 2 is appropriate only in the ideal case of data sets with uncorrelated errors, and δD 2 j is then given by statistical and point-to-point systematic errors added in quadrature.
For most experiments, additional information on the fully correlated normalization uncertainty δN n can be found, i.e., on a systematic shift common to the entire data set. Equation (1) is straightforwardly extended to account for such normalization uncertainties: Here, N n are the normalization factors, which can be either fitted along with the {a i }, or even determined analytically in each step of the minimization by solving ∂χ 2 /∂N n = 0 [42]. There are several equivalent methods of extending further the simple χ 2 function in Eq. (1) in the presence of K (n) sources of correlated systematic errors for a data point D j in data set n [42,46]. The numerically most efficient method treats the correlated systematic errors analytically in the optimization procedure like for the global normalization uncertainties discussed in (2). This avoids the construction and inversion of large covariance matrices used in the conventional approach. The result-ing χ 2 function has the form (assuming ω j = 1 for simplicity) [42,46] where Here, (δD (u) j ) 2 is the quadratic sum of the statistical and uncorrelated systematic errors, and β kj specifies the kth correlated systematic error of data point D j .
We note that most global analyses of unpolarized PDFs now start to include correlated systematic errors whenever this information is available from experiment. This is of much importance, as very precise PDF uncertainty studies are a key ingredient for reliably estimating new physics signals and standard model backgrounds at the Tevatron or the LHC. For the time being, and keeping in mind that our analysis is the first fully global one of this kind for polarized PDFs, we stick to an effective χ 2 function based on the simplest expression in Eq. (1). Also, most of the relevant experiments do not publish the full information on correlated systematic errors. Whenever a global normalization uncertainty is provided, we have explored the possibility of normalization shifts to improve the global fit by minimizing χ 2 according to Eq. (2). We have not found any significant improvements of the fit from this. Since data sets are continuously evolving and more and more precise information becomes available, the proper inclusion of correlated systematic errors is certainly needed in future global analyses of helicity PDFs.
As mentioned in the Introduction, there are considerable complications when statistical methods are applied to global QCD analyses based on a large body of diverse data and a theoretical model with many parameters {a i }. In particular, the statistical value of the definitions given in Eqs. (1)- (3) has been under considerable debate [29,30,41,42,43,46,47], since both the theoretical [T j ({a i })] and the experimental inputs [D j , δD j , . . .] are far from being ideally suited for statistical analysis. For instance, uncertainties inherent to the theoretical framework used to describe the data are notoriously difficult to quantify and usually correlated. In addition, it is often the case that even in a "best fit" to data, the values of χ 2 per data point for individual experiments vary considerably around unity, sometimes by much more than the expected amount 2/N (n) data . Therefore, conclusions drawn on a tolerable range of uncertainty from a certain increase ∆χ 2 , must be considered keeping in mind this complex situation.

B. Uncertainty estimates: Hessian and Lagrange multiplier methods
An important objective is to estimate uncertainties of the spin-dependent PDFs obtained from the χ 2 optimization. To this end one must study the behavior of the effective χ 2 ({a i }) function in the neighborhood of the global minimum χ 2 0 ≡ χ 2 ({a 0 i }), using reliable statistical methods rather than some subjective tuning of selected parameters. Basically two complementary tools have been put forward, refined to deal with the complexity of global PDF analyses, and applied to characterize uncertainties in the extraction of unpolarized PDFs. We pursue and compare both, the Hessian [43] and the Lagrange multiplier [42] methods, in our uncertainty estimates. We only briefly recall the main elements of the two approaches and highlight potential problems and limitations. For a detailed account we refer the reader to Refs. [41,42,43].
In the more standard Hessian approach, the exploration of the uncertainties associated with the fit is implemented through a Taylor expansion of Keeping only the leading quadratic terms, the increase ∆χ 2 can be written in terms of the Hessian matrix as where y i ≡ a i − a 0 i and the derivatives in Eq. (6) are taken at the minimum.
Global QCD fits are usually characterized by very disparate uncertainties in different directions of the multiparameter space, so that standard methods to evaluate H ij by finite difference, as implemented in, e.g., the Minuit package [48], tend to be numerically unstable and hence unreliable. To overcome such difficulties, a new iterative algorithm to compute the derivatives in (6) was devised in Ref. [41] and subsequently used in global analyses of unpolarized PDFs. We apply this improved Hessian method also in our studies.
Instead of working in the parameter basis {a i }, the Hessian H ij is expressed in terms of its N par eigenvectors v (k) i and eigenvalues ε k . The displacements y i in Eqs. (6) and (7) are replaced by a new set of parameters {z i } defined by [41,43] The {z i } are appropriately normalized by scale factors s j ∝ 1/ε j such that surfaces of constant χ 2 turn into hyper-spheres in {z i } space, with the distance from the minimum given by Large eigenvalues ε k correspond to steep directions in which χ 2 changes rapidly and the parameters are tightly constrained by the data, while small eigenvalues belong to directions where the parameters are only weakly determined.
Within the eigenvector representation, it is convenient to construct 2N par eigenvector basis sets of PDFs which greatly facilitate the propagation of PDF uncertainties to arbitrary observables O [43]. These basis sets S ± k are defined in {z i } space by and hence correspond to positive and negative displacements along each of the eigenvector directions by the amount T = ∆χ 2 still tolerated for an acceptable global fit. To estimate the error ∆O on a quantity O away from its best fit estimate O(S 0 ) it is only necessary to evaluate O for each of the 2N par sets S ± k [43], i.e., One must keep in mind that the propagation of PDF uncertainties in the Hessian method has been derived under the assumption that a first order, linear approximation is adequate. Of course, due to the complicated nature of a global fit, deviations, also from the simple quadratic behavior in Eq. (7), are inevitable, and error estimates based on the Hessian method are not necessarily always accurate.
A strategy based on Lagrange multipliers avoids all these assumptions and probes the uncertainties on any quantity O({a i }) much more directly. The result is a parametric relationship between one or more observables O j ({a i }) and the effective χ 2 function used to determine the goodness of the global fit. Its application to global QCD analyses was proposed in Ref. [42] not only to estimate uncertainties of observables depending on PDFs but also to test the range of applicability of the Hessian approach described above.
In practice, the method is implemented by minimizing a function with respect to the set of PDF parameters {a i } for fixed values of the Lagrange multipliers {λ j }. Each multiplier is related to one specific observable O j , and the choice λ j = 0 corresponds to the best fit S 0 . By repeating this minimization procedure for many values of λ j one can map out precisely how the fit to data deteriorates as the expectation for the observable O j is forced to change.
Here, contrary to the Hessian method, no assumptions are made regarding the dependence of χ 2 or the observable O j on the parameters {a i } of the fit. The Lagrange multiplier method also generates a large set of sample PDFs along the curve of maximum variation of the observable(s) used in Eq. (12). In principle, the Lagrange multiplier method is superior to the Hessian approach for reliably estimating uncertainties. For a given ∆χ 2 it finds the largest range of O j ({a i }) allowed by the data used in the global fit and the theoretical ansatz, independent of the approximations involved in the Hessian method. The entire parameter space {a i } is explored in minimizing Eq. (12), not necessarily limited to the vicinity of the best fit {a 0 i }. In practice, a large number of global fits is required to map out all χ 2 profiles of interest. Thanks to our novel Monte Carlo sampling technique to be described below, this is no longer a serious limitation computationally. The only drawback to this method is that it requires access to the full machinery of global fitting to estimate uncertainties of a given observable of interest, contrary to the Hessian method for which the eigenvector PDF sets S ± k make it very simple to propagate PDF uncertainties to arbitrary observables.
In this paper, Lagrange multipliers will mainly be the method of our choice for estimating uncertainties and to monitor to what extent the approximations involved in the Hessian approach are justified.

C. Computational technique: The Mellin method
As is evident from the previous Subsections, the extraction of PDFs in a global QCD analysis of a large body of data at NLO accuracy requires an extensive number of time consuming computations of the corresponding observables in each step of the χ 2 minimization procedure. The large number of parameters specifying the functional form of the PDFs in the fit, typically of O (20), and the need for a proper assessment of their uncertainties, add to this. Assuming a representative analysis with about 500 data points, 5000 iterations to reach the optimum set of parameters, and a modest computing time of 1 second per cross section calculation at NLO, one easily realizes that computational improvements are very much in demand.
We stress at this point that approximating NLO corrections by parameterizing them by a K-factor, K ≡ σ NLO /σ LO , which is a possible strategy in the spinaveraged case in order to speed up the analysis [29,30,39,40,46,47], is not a viable option in the polarized case. In the latter, the parton distributions as well as the hard scattering cross sections may have nodes and change sign within the kinematic region of interest. As a result, differ-ent partonic subprocess contributions can have very different NLO corrections that are never well approximated by a single K-factor. This problem is even more prominent at the present time when the spin-dependent parton distributions, in particular the gluon distribution, are still poorly known. Also, from a more theoretical point of view, the NLO corrections are expected to decrease the factorization/renormalization scale dependence of the calculation. In the ideal case that the NLO cross section has relatively little scale dependence, this would imply that the K-factor inherits the full LO scale dependence and thus cannot serve well as a proxy for the NLO corrections. We thus do not really see a useful workaround that would avoid inclusion of the full NLO calculation in the global analysis.
The required numerical calculations also involve the scale evolution of the PDFs from some initial scale µ 0 ∼ O(1 GeV) to each of the scales relevant for the data points. The evolution is governed by a well-known set of integro-differential equations [49,50] that can be solved analytically after an integral transform from Bjorken x space to complex Mellin N -moment space. The Mellin transform of a generic function ϕ depending on the lightcone momentum fraction x is defined as and its inverse reads Here C N denotes a suitable contour in the complex N plane that has an imaginary part ranging from −∞ to +∞ and that intersects the real axis to the right of the rightmost pole of Φ(N ). In practice, it is beneficial to choose the contour to be bent at an angle < π/2 towards the negative real-N axis [51]. The integration in (14) can then be very efficiently performed numerically by choosing the values of N as the supports for a Gaussian integration. The transformation (13) has the welcome property that convolutions factorize into ordinary products, which greatly simplifies calculations based on Mellin moments. It can be carried out analytically not only for the splitting functions governing the scale evolution of the PDFs but even for the partonic hard scattering cross sections for both inclusive and semi-inclusive DIS. The latter case requires straightforward extensions of Eqs. (13) and (14) to double transformations as was discussed in Ref. [44,52]. This reflects the fact that SIDIS depends on two scaling variables, the momentum fractions x and z taken by the struck parton from the parent nucleon and by the finalstate hadron from the scattered parton, respectively. The usefulness of the Mellin technique was demonstrated in practice in a global analysis of helicity-dependent PDFs using all available DIS and SIDIS data [36].
However, the inclusion of observables in hadron-hadron collisions or in less inclusive reactions in lepton-hadron scattering, which are crucial for determining, e.g., the gluon distribution, require a more elaborate framework. They involve multiple convolution integrals between one or more PDFs, partonic cross sections, and, depending on the process, fragmentation functions. More importantly, they typically depend on several kinematic variables, so that there is no direct way of taking Mellin moments of the cross section under which it would become a simple product of Mellin moments of PDFs and partonic cross sections. An example is single-inclusive pion production in proton-proton collisions, pp → πX, at measured transverse momentum p T and rapidity η of the pion. While S is the typical scaling variable of the process, taking Mellin moments of the cross section in x 2 T does not lead to a simple expression involving the moments of the PDFs. An efficient computational scheme that allows to circumvent this complication has been devised in Ref. [44,45]. Its main feature is to use the inverse Mellin transforms of the PDFs in the theoretical expression, which makes it possible to store all numerically time consuming calculations involving the lengthy and complicated expressions for the underlying hard scattering cross sections on large "look-up tables" or "grids" in Mellin moment space. Here, we review the technique of [44] and describe a method to compute these grids very efficiently using Monte-Carlo sampling techniques. The latter allows to generalize the Mellin moment technique beyond the single inclusive processes considered in Ref. [44].
A typical observable of interest, the spin-dependent cross section for pp → πX at RHIC, has the following general factorized structure: Here we have suppressed the dependence on kinematic variables other than momentum fractions, and on the factorization and renormalization scales. In (16), ∆f i and D k denote the spin-dependent parton distribution and spin-averaged fragmentation functions for flavor i and k, respectively, and d∆σ ij→kX the relevant partonic hard scattering cross sections. S represents a "measurement function" that defines details of the observable ∆σ | bin in each bin, such as the kinematical ranges explored and the relevant experimental cuts. The symbols ± in Eq. (15) denote the helicity combinations of the colliding longitudinally polarized protons defining the cross section ∆σ | bin .
Following the strategy developed in Ref. [44], we replace the PDFs in Eq. (16) by their representations as Mellin inverses, Eq. (14), and subsequently interchange integrations to obtain where and Here and in the following, we assume that only the spindependent PDFs are subject to the global analysis and that the fragmentation functions D k are known.
In the next step, one can now calculate the quantities d∆σ ijk (N, M ) in Eq. (19) prior to the actual fit, as they do not depend on the PDFs, and store them in large lookup tables in the Mellin variables M and N along the contours C M,N . In practice, it is convenient to choose M and N on the supporting points of a Gaussian integration, see Ref. [44] for details. As can be seen from (19), this effectively amounts to computing the NLO cross sections for all partonic subprocesses using complex "dummy PDFs" . Even after making use of all symmetry relations, e.g., among different subprocesses ij → kX, setting up all look-up tables for a typical N × M grid size of 64 × 64 is a rather time-consuming step in practice, although it has to be done only once. This is where our new Monte-Carlo sampling technique has considerable advantages over a "brute-force" computation.
First, we recall that a Monte-Carlo algorithm reduces the multiple integrations in (15) to a finite sum over random "events" n, with weights w (n) and scaled by the number of iterations, κ, used to generate large samples of I events. The weight includes the dependence on the parton distribution and fragmentation functions and the hard scattering cross sections for each event. In Eq. (20) we keep the measurement function S separate from w (n) , as it is usually implemented only after an event has been generated. The most efficient strategy to compute all the required d∆σ ijk (N, M ) in Eq. (19) is to choose a suitable set of trial PDFs ∆f i and to perform a single highstatistics Monte-Carlo integration calculation of the cross section in Eq. (16). During the calculation one stores, for each event n, the momentum fractions x ijk for all of the subprocesses ij → kX. Renormalizing each weight by removes its dependence on the assumed set of PDFs. The d∆σ ijk (N, M ) are then straightforwardly obtained as ijk , which corresponds to a profile of the integrand in Eq. (19), allows to simultaneously compute the moments d∆σ ijk (N, M ) for all N, M , without having to do an individual Monte-Carlo integration for each of them as in (19). This greatly reduces the computational burden. To give an example, for our global analysis of polarized PDFs, which currently includes about 50 very timeconsuming NLO calculations of single-inclusive jet and pion production at RHIC in each step of the χ 2 minimization, all the necessary grids in Eq. (23) can be computed within approximately one day of CPU time on a single standard PC with a CoreDuo processor running at 2 GHz. Once these grids are available, a full PDF fit can be finalized in about 15-30 minutes. A detailed PDF uncertainty assessment, which requires a very large number of χ 2 minimizations, can then be performed easily in about 1-2 weeks. We note that in practice one does not even need to physically store the set x ijk , since the summations (23) can in fact be done simultaneously with the generation of the events n.
The formula in Eq. (20) applies to any computation of the theoretical cross section when a Monte-Carlo integration is employed. This is the case for analytical NLO calculations of single-inclusive jet or hadron rates [53,54] where Monte-Carlo integration techniques are just used to perform the convolution with the parton distributions, or for NLO parton-level Monte-Carlo generators evaluating general infrared-safe observables [55,56]. Therefore our method based on Mellin moments and Monte-Carlo sampling techniques outlined in Eqs. (13)-(23) above is completely general and can be straightforwardly applied to any observable for which a perturbative QCD description is available. It can be applied to global analyses of polarized PDFs, pursued in this paper, extractions of fragmentation functions, see Refs. [37], but equally well to analyses of ordinary spin-averaged PDFs incorporating Tevatron and future LHC data consistently at NLO or beyond. We note that for the latter case an alternative method, called "fastNLO" [57], has been developed, which also allows to include NLO corrections to hadronic scattering in a fast and efficient way in a global PDF analysis. Like our method, it amounts to preparing huge look-up tables that contain all the time-consuming NLO calculations prior to the actual fit. In the "fastNLO" case, this is done in x-space, and interpolations between the various support points in x are used. This step appears to be rather time-consuming and the method has so far been tested only for inclusive jet production data. In addition, the evolution of the PDFs needs to be performed as a separate calculation, whereas in our approach it is immediately included in Mellin moment space.
This latter advantage can in fact be used for further improvements of our Mellin technique. The scale dependence of the PDFs, which we have so far suppressed, can be schematically written as Here, E ii ′ (N, µ 0 , µ) denotes the appropriate evolution matrix from the initial scale µ 0 where we parameterize the PDFs to the scale µ relevant for a certain observable O. Inserting (24) into Eq. (17) allows to absorb also the scale evolution of the PDFs into the pre-calculated Mellin grids by extending Eq. (23) to (25) so that the luminosity function in Eq. (18) now only refers to the PDFs at the initial scale µ 0 , L ij (N, M ) = ∆f i (N, µ 0 ) ∆f j (M, µ 0 ). An advantage of this re-shuffling is that now all dependence on the scale µ is contained in the look-up tables (25), eliminating the need to perform the scale evolution later in the fitting code. More importantly, if the experimental observable used in the global fit involves an integration over a bin in, say, the transverse momentum p T of a jet, which also acts as the factorization scale, the µ dependence of the PDFs is correctly taken care of in the integration. While we have not made use of this particular improvement in the present analysis, we expect it to be useful in the future for further optimizing the performance of the global analysis code.

III. GLOBAL ANALYSIS AND UNCERTAINTY ESTIMATES FOR POLARIZED PDFS
In this Section we give a detailed account of the first global analysis of polarized PDFs presented in Ref. [28] which in the following will be referred to as "DSSV". We first discuss the data selection and the determination of the best fit, which we compare to the fitted data. We then focus on the studies of uncertainties, including a comparison of the Lagrange multiplier method used in [28] and the more standard Hessian error matrix approach. For the latter we present a new family of eigenvector PDFs, as described above, which greatly facilitates estimates of the PDF uncertainties of any given observable of interest.
scription of the available hard scattering data. The data sets for the spin asymmetries we use in our analysis are listed in Tab. I, along with the number of data points fitted. We minimize the effective χ 2 function in Eq. (1). Attempts to further improve the global fit by introducing normalization shifts for each experiment and minimizing χ 2 according to Eq. (2) were to no avail. All theoretical spin asymmetries in Eq. (1) are calculated at NLO, using the appropriate factorized leading-twist expressions. We use the MS scheme throughout, and all our results for the polarized PDFs will refer to this scheme.
In case of inclusive DIS, the asymmetries are computed, as in our previous analyses [31,36], as the ratios between the polarized and unpolarized structure functions, with and the corresponding expression for F 1 (x, Q 2 ), both computed at NLO using the appropriate coefficient functions [50,58]. For DIS off a deuteron target, we take into account the ω D = 5.8% D-wave state probability in relating the g 1 structure function of the deuteron to those of proton and neutron: The extension of the expression in (26) to SIDIS is straightforward, using the NLO coefficient functions given in [59].
For the case of pp scattering, the spin asymmetries are computed using the generic expression in Eq. (16) at NLO, and its spin-averaged counterpart. The NLO corrections for high-p T single-jet and hadron production can be found in [53,54], respectively. We have always chosen the renormalization and factorization scales as the transverse momentum p T of the observed final state, a choice that leads to good agreement of NLO calculations [37,53,54] and experimental data from RHIC [22,60] in the spin-averaged case. For the computation of the unpolarized cross sections, we always use the NLO unpolarized PDFs of Ref. [47]. Whenever fragmentation functions are needed, as is the case for SIDIS and the RHIC pp → πX data, we use the DSS set [37] for pions, kaons, and unidentified charged hadrons, which was recently obtained from a global analysis of hadron production data. The use of up-to-date fragmentation functions that are consistent with HERMES and RHIC unpolarized data and have quantified uncertainty estimates [37] is a crucial ingredient of our analysis. It is a major difference with respect to the polarized PDF analysis in Ref. [36], where also SIDIS data were included. In the computation of the χ 2 contribution from SIDIS asymmetries we have taken into account the uncertainty coming from the set of fragmentation functions. In practice, this was done by determining for each data point the maximum variation of the theoretical estimates T j in Eq. (1) due to the FFs within their own uncertainty ranges quoted in [37]. This variation was treated as an additional uncertainty and added in quadrature to the experimental error δD j . Since at least for pions the FFs are fairly well determined, this amounts to additional uncertainties of a few percent at most. We do note, however, that this uncertainty estimate does not reflect possible systematic problems in interpreting the available spin-averaged kaon SIDIS data within a leading-twist pQCD analysis, a point to which we will return later.
We do not systematically include any higher twist (or, more generally, power-suppressed) contributions in the theoretical calculations that enter our analysis, neither for the inclusive and semi-inclusive DIS observables, nor for proton-proton collisions. We employ a cut of Q 2 > 1 GeV 2 for all DIS and SIDIS data, and of p T > 1 GeV for the RHIC high-p T polarized pp data. These cuts have the purpose to exclude regions where contributions beyond the leading twist, factorized framework of pQCD become crucially important. For example, for the SIDIS and RHIC pp data, it is known that the underlying unpolarized cross sections in the same kinematic domain, i.e., for scales above 1 GeV, can be quite successfully described by pQCD [37].
That said, it is known [61] that the relation between A 1 and g 1 /F 1 in Eq. (26) is corrected by a factor (1 + γ 2 ) ≡ (1 + 4M 2 x 2 /Q 2 ) on the right-hand-side (r.h.s.), corresponding to a target mass correction. It has been pointed out [62] that this correction is non-negligible in some kinematic regimes accessed by the lower energy fixedtarget experiments, typically at relatively low Q 2 and high x. We have therefore corrected for this factor where necessary. Specifically, since our set of polarized PDFs is defined and related to the measured asymmetries through the leading-twist relations (26) and (27), our choice is to multiply data sets that are given in terms of measured g 1 /F 1 by the factor (1 + γ 2 ), but to leave data sets for measured A 1 unchanged. The resulting data are confronted with the NLO leading twist calculation. We stress that various other choices have been adopted in the literature [32,33,34], using, for example, parameterizations of experimental data for F 2 (x, Q 2 ) and R(x, Q 2 ) = F 2 /(2xF 1 ) − 1. All choices are equivalent in the strict leading twist sense, but will in general differ in the amount of power corrections needed to describe the data at lower Q 2 and/or higher x. As we shall see below, we find that for our choice without inclusion of any power corrections an excellent description of all data sets within our specified cuts is achieved (cf. also Tab. I), with no visible discrepancies even at rather low scales. For the time being, we thus regard our choice as one that "empirically" alleviates the need for power-suppressed corrections in spin asymmetries. Despite some significant progress on the analysis of higher twist effects in polarized DIS [33,62], further detailed investigations will be needed in this area, including implementation of the full target-mass corrections [63], the effects of yet higher orders [64,65] which generally reduce the need for higher twist contributions [64,66], and so forth.
In order to find the optimal helicity-dependent PDFs from a χ 2 minimization, we parameterize them at an input scale of µ 0 = 1 GeV by the following flexible form [28] x∆f i (x, with independent parameters for ∆u + ∆ū, ∆d + ∆d, ∆ū, ∆d, ∆s ≡ ∆s, and ∆g (note that here and in the following we interchangeably use ∆f u = ∆u, ∆f g = ∆g etc. to denote the polarized PDFs). The minimization is carried out with respect to the set of parameters The PDFs are evolved to the scales µ > µ 0 relevant in experiment. The particular functional form and the value for µ 0 are not too crucial, as long as the parameterization is flexible enough to accommodate all hard scattering data within their ranges of uncertainties. The ansatz (28) deviates considerably from the standard form used in fits to DIS data only [31,32,33,34,36], where γ i = η i = 0, inasmuch as it allows the PDFs to develop nodes and to deviate from an SU(3) flavor symmetric sea. As will be seen from our results presented below, this extra freedom in parameter space {a i } is crucial in a comprehensive analysis of DIS, SIDIS, and RHIC pp data.
In addition to the much more flexible input parameterization proposed in the preceding paragraph, we have repeated the analysis with alternative parameterizations, some of them even more flexible than the one we choose. For example, we have chosen the powers of x in the last two terms different from 1/2 and 1, even allowing the fit to vary them. None of these modifications resulted in any significant improvement in the quality of the fit to data, or changes of the uncertainty bands. This indicates that the present data is not really able to discriminate between various forms of the input distributions, as long as a sufficiently flexible choice is made. Therefore our present choice does not introduce large additional uncertainty in that respect.
Analyses of polarized PDFs routinely use constraints that can be derived from baryonic semi-leptonic β-decays under the assumption of SU(2) and SU(3) flavor symmetries [67]. These relate combinations of the first moments of the PDFs to the constants F and D parameterizing the β-decay rates. We make use of these constraints in our present analysis; however, rather than imposing the exact SU(2) and SU(3) flavor symmetry relations, we allow for deviations in our fit within the uncertainty ranges of the F, D values. Specifically, we set ε SU(2,3) parameterize the departures from exact SU(2) and SU(3) symmetries, and where we use the latest values F + D = 1.269 ± 0.003 and 3F − D = 0.586 ± 0.031 (see, e.g., Ref. [11]). As a practical matter, we trade the input parameters N u+ū and N d+d in Eq. (28) for ε SU(2,3) and fit the latter. Here the relative uncertainties of F + D and 3F − D are assumed to represent the typical ranges of ε SU(2,3) ; we use them to include the deviations from ε SU(2,3) = 0 as additional contributions to χ 2 , similarly to the case of normalization uncertainties shown in the first term of Eq. (2). We note that the relative uncertainties of F +D and 3F −D are rather modest and may not fully reflect the actual breaking of the SU(2) and, in particular, SU(3) symmetries, for which larger breaking effects have been discussed in the literature [68]. This issue may need to be revisited in the future. For now we note that as a result of this the PDFs in our fits will naturally have a tendency to have relatively small ε SU (2,3) .
Rather than determining also the strong coupling α s in the global fit along with the PDFs, we take Λ QCD = 334.2 MeV for n f = 4 flavors from the analysis of unpolarized PDFs in Ref. [47]. The scale dependence of α s is computed by numerically solving its renormalization group equation at NLO accuracy. The charm and bottom quark thresholds are crossed at m c = 1.43 GeV and m b = 4.3 GeV, respectively. As already mentioned, the scale evolution equations for the PDFs are solved analytically in Mellin moment space by explicitly truncating the solutions at NLO. Likewise, all observables used in our fit are computed consistently at NLO accuracy in the MS factorization scheme. All quarks are treated as massless, and charm and bottom PDFs are turned on in the evolution at Q = m c,b . We note that for all presently available spin-dependent observables heavy quarks play a negligible role, and, for the time being, we can safely refrain from introducing more elaborate variable flavor number schemes which take into account quark mass effects near threshold (see, e.g., Refs. [29,30]). The parameters {a 0 i } representing our best global fit of polarized parton densities ∆f i in Eq. (28), henceforth denoted as "set S 0 ", are given in Tab. II. A few additional remarks are in order here. The currently available data do not fully constrain the entire x dependence of ∆f i imposed in Eq. (28), and we are forced to make some restrictions on the parameter space {a i }. For the sea quark and gluon densities we set γ i = 0 in Eq. (28). This only marginally limits the freedom in the functional form and still allows nodes. In addition, we tie the small x behavior, represented by the α i in Eq. (28), of ∆u + ∆ū and ∆d + ∆d to that of ∆ū and ∆d, respectively, which is reasonable as sea quarks likely dominate in this region. The parameter αs always came out close to αd, so we set them equal. Since the SIDIS data are not yet sufficient to distinguish ∆s from ∆s, we assume ∆s = ∆s throughout. In general, the x → 1 behavior of all PDFs is rather unconstrained as is the case even for current sets of unpolarized PDFs. This is because there are no data sensitive to x 0.6. To avoid problems with the fundamental positivity constraint we make sure that all ∆f i vanish at large x at least as fast as the f i of our reference set of unpolarized PDFs [47], which puts constraints on the β i . Choosing any other recent set of unpolarized PDFs, like CTEQ6 [46], does not alter our results. We note that implementing the positivity constraint at the level of PDFs [69], i.e., is strictly valid in LO only, but in the MS scheme is sufficient to guarantee (32) also in NLO. The parameters β u+ū and β d+d come out very close to their corresponding values for the unpolarized case, implying (∆u + ∆ū)/(u +ū) → const. as x → 1, and likewise for (∆d + ∆d)/(d +d). Since the other β i in Eq. (28) are only very weakly determined by the fit, we fix them to their preferred values within the positivity constraints whenever we examine uncertainties of the PDFs, in order to avoid extremely flat directions in parameter space where χ 2 varies only very slowly. Therefore, and as is in general the case in PDF studies, our uncertainty estimates for helicity-dependent PDFs are valid only in the x-region explored by experiment. Notice that the near equality of η d+d and ηd is not imposed but a result of the fit (in fact, the actual values for these parameters before rounding are 98.9384 and 98.9354, respectively). In total this leaves us with 19 free parameters in the fit [or 17, if we fix also β u+ū and β d+d ], which we include later on also in our uncertainty estimates. We tried to relax the imposed constraints discussed above, but found that present data, i.e., the effective χ 2 function, are not really sensitive to them. In Tab. II we have converted the fitted values for ε SU(2,3) , defined in Eqs. (29) and (30), back to N u+ū and N d+d for convenience. For the optimal DSSV fit (set S 0 ) we find which corresponds to only very minor violations of the canonical constraints on the first moments ∆Σ u − ∆Σ d and ∆Σ u + ∆Σ d − 2∆Σ s assumed in most fits so far. As we have discussed above, the smallness of ε SU(2,3) is not really a surprise in view of the relatively small nominal uncertainty of the F + D and 3F − D values in Eqs. (29) and (30). If correct, the small value for ε SU(3) has interesting implications on the behavior of the ∆s(x) = ∆s(x) distribution in the best fit, as we shall see later.

B. Comparison to fitted data
The total χ 2 of the best fit S 0 is 392.6 for 467 data points used in our NLO global analysis. We list in Tab. I also the individual χ 2 values for each experiment. As one can see, there are only very few cases where the χ 2 /N (n) data is significantly larger than one. In each case, this is due to large fluctuations of some of the data points in that particular set which are impossible to accommodate in the fit. Figure 1 shows the comparison of our fit to the fitted DIS data, while the comparison to the SIDIS data is shown in Fig. 2. Notice that in Fig. 1 the plots are generically labeled as asymmetries "A 1 "; however, in the case of the E143, E155, CLAS and Hall A data, they actually show the reported structure function ratios and are compared to the DSSV estimates for the asymmetries, divided by the factor (1+γ 2 ). The overall agreement of the experimental sets in the global analysis is excellent. All data can be very satisfactorily described by a universal set of polarized PDFs as is assumed by the fundamental factorization theorem. The agreement with the RHIC pp data is equally good; we have shown the corresponding comparison in our previous paper [28] and will come back to it in the next Subsection. Figures 1 and 2 also show the results obtained for the set of polarized PDFs of [36] in the following labelled as "DNS". Apart from the fact that not all of the present data sets were available at the time of [36], a main difference between the two analyses resides in the fragmentation functions used when including the SIDIS spin asymmetries in the fit. To illustrate this point, we have used here the new fragmentation functions of [37] also for the calculations with the DNS set [36]. As can be seen from Fig. 2, this leads to significant differences, in particular, in the kaon and, to a lesser extent, the isospin sensitive (proton target) SIDIS asymmetries. This is mainly due to the strange fragmentation functions, which directly affect the strange quark polarization, and also to differences in the light sea quark distributions. Figure 1 shows that there is, however, little difference between the two sets as far as the inclusive DIS asymmetries are concerned. The changes in the strange quark and other sea polarizations are thus compensated here, either mutually or by the other parton distributions. Figure 3 shows the extracted polarized PDFs at Q 2 = 10 GeV 2 , along with estimates of their uncertainties for the Hessian and Lagrange multiplier methods, both for a tolerance of ∆χ 2 = 1. The results for the Lagrange multiplier method were already shown in our previous paper [28] along with a more conservative estimate of the PDF uncertainties based on ∆χ 2 /χ 2 = 2%. For this method, the estimates were obtained by varying the first moments of the distributions, truncated to the region of momentum fractions 0.001 ≤ x ≤ 1 covered by the data included in the fit. The distributions corresponding to the maximum variations of the truncated moments for a given increase ∆χ 2 were then taken as faithful estimates of the range of variation of the PDFs. In the case of the polarized gluon distribution, this procedure was found to be not adequate [28] because of the fact that there is a significant amount of rather precise proton-proton collision data constraining the gluon density in however a relatively narrow region of momentum fraction, 0.05 x 0.2. In this way the variations of the gluon distribution's integral in the full region 0.001 ≤ x ≤ 1 tend to produce distributions that favor the variations outside the pp kinematic region, misrepresenting the uncertainties inside. In order to circumvent this problem, we performed variations of the integral of the gluon distribution in three different x regions, 0.001 ≤ x ≤ 0.05, 0.05 ≤ x ≤ 0.2, and 0.2 ≤ x ≤ 1, allowing them to jointly contribute a change in χ 2 of ∆χ 2 = 1. Clearly, the choice for these regions and the way they share the increase in ∆χ 2 is not unique. In order to specifically focus on the xregion accessed at RHIC, we also performed a dedicated study of the truncated moment of ∆g in this region, allowing variations of ∆χ 2 = 1 from this region alone. The results for the truncated moments of our polarized PDFs,

C. Extracted PDFs, their uncertainties, and their physics
are given in Tab. III. Inspection of Fig. 3 and Tab. III reveals that the Hessian and the Lagrange multiplier methods yield fairly similar ∆χ 2 = 1 uncertainties, except for the spindependent gluon distribution, for which the Lagrange multiplier approach yields a still significant but generally smaller uncertainty than the one predicted by the Hessian method using Eq. (11). The agreement between the two often becomes better when the observable is better constrained by the data, as is the case for the integral at Q 2 = 10 GeV 2 and their uncertainties for ∆χ 2 = 1 obtained with the Lagrange multiplier and the Hessian methods. For future reference, we also recall the results for the Lagrange multiplier method obtained in [28] under the assumption ∆χ 2 /χ 2 = 2%, which are to be considered more realistic estimates of the uncertainties. In the last line, ∆g RHIC  of ∆g over only the x-range probed at RHIC, or for the actual physical observables that determine ∆g. As an example, in Fig. 4 we show the estimated uncertainties for the double-longitudinal spin asymmetry, for pp → π 0 X at RHIC, where the superscripts denote the helicities of the incoming protons, computed with both the Lagrange multiplier and the improved Hessian approaches. As can be seen, the two give very similar re-  FIG. 2: Same as Fig. 1, but for the semi-inclusive DIS asymmetries [10,14,15,16]. In all calculations the fragmentation functions of [37] have been used.
sults. This feature can be traced back to correlations between the parameters, in the sense that some of them can compensate variations forced in the others. We note that such kinds of correlations are fully accounted for in the Lagrange multiplier approach, whereas it is not generally clear how well are they represented by the approximated Hessian matrix. We shall investigate the distinctive features between the two methods later, but will focus first on the physics aspects related to our extracted polarized PDFs. Table IV shows the evolution of the central values for the truncated first moments ∆f 1,[0.001→1] i with Q 2 . ∆Σ denotes the quark singlet combination, i.e., the sum of all quarks and anti-quarks. We also show the evolution of the full first moments ∆f 1 i . These obviously rely on an extrapolation of the PDFs to x-values outside the measured region, and it is difficult to estimate the uncertainty associated with this. Total up and down distributions: ∆u + ∆ū and ∆d + ∆d, which inclusive DIS probes primarily, are the by far best determined distributions. Their uncertainty bands are very narrow, see Fig. 3, and also our results agree very well with the determinations in previous analyses [31,32,33,34,36]. We note that recent lattice QCD results [70] of the full first moments ∆Σ u ≡ ∆u 1 + ∆ū 1 and ∆Σ d ≡ ∆d 1 + ∆d 1 (albeit excluding disconnected diagrams) also agree very well with the values we extract, which may shed light on the validity of assumed extrapolations of the parton distribution functions to small x.
We have mentioned earlier that in our fit R u ≡ (∆u + ∆ū)/(u +ū) and R d ≡ (∆d + ∆d)/(d +d) become constant in the "valence region" as x → 1, where the sea quark contributions become small. Figure 5 shows the ratios R u , R d along with the most relevant experimental data. The information at the highest values of x comes from the Jefferson Laboratory Hall-A experiment [12]. As one can see, our R u goes to unity at high x, which is consistent with expectations in relativistic constituent quark models [71], but also in perturbative QCD, using power counting and hadron helicity conservation [72]. We furthermore find that R d remains negative in the region where it is constrained by data and presently shows no  tendency to turn towards +1 at high x. The latter behavior would be expected for the pQCD based models. We note that it has recently been argued [73] that the upturn of R d in such models could set in only at relatively high x, due to the presence of valence Fock states of the nucleon with nonzero orbital angular momentum that produce double-logarithmic contributions ∼ ln 2 (1 − x) in the limit of x → 1 on top of the nominal power behavior. The corresponding expectation is also shown in the figure. In contrast to this, relativistic constituent quark models predict R d to tend to −1/3 as x → 1, perfectly FIG. 4: Uncertainties of the calculated A π 0 LL at RHIC in our global fit, computed using both the Lagrange multiplier and the Hessian matrix techniques. We also show the corresponding PHENIX data [23]. consistent with the present data. Light sea quark polarizations: The light sea quark and anti-quark distributions turn out to be better constrained now than in previous analyses [36], thanks to the advent of more precise SIDIS data [10,14,15,16] and of the new set of fragmentation functions [37] that describes the observables well in the unpolarized case. Figure 6  and R d = (∆d + ∆d)/(d +d), shown by the solid curves, along with the data from [12,13,15], as shown in [73]. The dashed lines present the predictions based on power counting and perturbative QCD, taking into account nonzero orbital angular momentum Fock states [73]. We note that these results, as well as some of the experimental data, are for the ratios ∆u/u and ∆d/d rather than for Ru and R d , which, however, makes little difference for the large x values we consider here. The arrows show the expectations for ∆u/u and ∆d/d in relativistic constituent quark models [71].   of the fit.
Of particular physics interest is a possible flavor symmetry breaking in the light sea, i.e., ∆ū = ∆d, given the well-established significant difference betweenū andd in the spin-averaged case [29,30]. Figure 3 indeed clearly points to a largely positive ∆ū distribution, but a negative (and larger) ∆d. Figure 7 specifically shows the difference x(∆ū − ∆d), which is positive within uncertainties. Note that we show both the ∆χ 2 = 1 and the more conservative ∆χ 2 /χ 2 = 2% uncertainty bands here.
The pattern of symmetry breaking in the light antiquark sea polarizations shown by Figs. 3 and 7 has been predicted at least qualitatively by a number of models of nucleon structure. A simple intuitive consideration of the Pauli principle roughly gives the observed picture: if valence-u quarks primarily spin along the proton spin direction, uū pairs in the sea will tend to have the u quark polarized opposite to the proton. Hence, if such pairs are in a spin singlet, one expects ∆ū > 0 and, by the same reasoning, ∆d < 0. Expectations based on the Pauli principle have been made quantitative in [74] and the "valence" scenario of [31], and the resulting predictions are shown by the dot-dashed line in Fig. 7. They tend to lie somewhat higher than our extracted ∆ū − ∆d, but are certainly qualitatively consistent, given the still relatively large uncertainties. The same is true for the case of the chiral quark soliton model [75], represented by the dotted line in the figure. Within the large-N c limit of QCD on which this model is based, one in fact expects |∆ū − ∆d| > |ū −d|. As comparison of our extracted x(∆ū − ∆d) with the result of [46] for x(d −ū) in Fig. 7 shows, one can presently not yet decide whether this expectation is fulfilled. Predictions for ∆ū − ∆d have also been obtained within meson cloud models [76]; it has been found in [77] that also here a flavor asymmetry of similar size is possible. Finally, also statistical parton models [35,78] obtain a similar size of ∆ū − ∆d. We note that predictions for the individual ∆ū and ∆d, where available, agree on ∆ū > 0, ∆d < 0, consistent with our results in Fig. 3, but may differ in the relative size of the distributions. For example, the results of [31,74] have |∆d| > ∆ū, as in Fig. 3, while the statistical models find the two distributions to be of more equal absolute size. Strange quark polarization: The polarization of strange quarks has been a focus since the very beginning of the proton spin crisis. The reason is that in the parton model and assuming SU(3) symmetry (see Sec. III A) one has where the ∆Σ f are as defined in Eq. (31) but now for arbitrary scale Q, and ∆Σ is the total quark and antiquark spin contribution to the proton spin. If the latter is found to be small experimentally, ∆Σ ∼ 0.25, the implication is that strange quarks make a significant negative contribution to the proton spin. Indeed, most fits to only inclusive DIS data have preferred a large and negative strange quark polarization. The same was found in Ref. [36], even though here the SU (3)  The difference between x∆ū and x∆d at Q 2 = 10 GeV 2 , along with the uncertainty bands for ∆χ 2 = 1 and ∆χ 2 /χ 2 = 2%. The dot-dashed and dotted lines show the predictions of the valence scenario of [31] and the chiral quark soliton model of [75], respectively. We also show the result obtained in an earlier global analysis [36] of DIS and SIDIS data (light dotted), for which the fragmentation functions of [37] were not yet available. The dashed line displays for comparison the flavor asymmetry x(d −ū) in the spin-averaged case, using the PDFs of [46].
was not enforced.
At variance with these results, the best fit in our present analysis has a polarized strange distribution ∆s that is positive at large x, but negative at small momentum fractions. Before we discuss the origin and significance of this result, we note that a prerequisite for it is that we have adopted a more flexible parameterization for the strange quark distribution in this work, which permits a node. This is again in contrast with the previous fits in which the initial ∆s always had the same sign for all x. We have assumed however ∆s = ∆s, since the fit is unable to discriminate strange quarks from anti-quarks. This is really an assumption: unlike the spin-averaged case where the distributions s ands will be rather similar (the integral of s −s has to vanish), there is a priori no need for ∆s and ∆s to have the same size or even the same sign.
Qualitatively, the main features of our extracted strange sea distribution arise in the following way: the (kaon) SIDIS data, within the leading-twist framework we employ, turn out to prefer a small and likely positive ∆s at medium x, while inclusive DIS and the constraints from β-decays demand a negative integral of ∆s and so force ∆s to turn negative at low x. Given the importance of ∆s, we address these constraints and their significance and implications in more detail in the following.
We start by we analyzing the behavior of the trun- cated first moment, ∆s 1,[0.001→1] , around the minimum defining the best fit. Figure 8 shows the increase of χ 2 of the fit against variations of ∆s 1,[0.001→1] , along with the partial contributions of the various data sets. Evidently, the best fit has a truncated moment close to zero and only slightly negative, as we also saw in Tab. III. The shape of ∆χ 2 around the minimum is dominated by the SIDIS data, and here primarily by the data for kaon production. All other data sets, pion SIDIS, inclusive DIS, and RHIC pp data, play less important roles, as expected (here one has of course to keep in mind that the impact of individual data sets seen in the Lagrange multiplier scans is always estimated in the "presence" of the other data sets, and therefore should not be construed as an independent fit result). As can be seen from Fig. 8 and Tab. III, the truncated moment of ∆s remains close to zero if changes of ∆χ 2 = 1 are permitted. For the more realistic choice ∆χ 2 /χ 2 = 2%, one finds that a much larger range of ∆s 1,[0.001→1] is allowed, extending from significantly negative to positive values. The size and even the sign of the considered truncated moment are, therefore, presently not well constrained. Nonetheless, there is a trend for ∆s(x) to be positive at medium x ∼ 0.1, even for the choice ∆χ 2 /χ 2 = 2% (see Fig. 2 of [28]). We note that the COMPASS experiment has recently presented a LO extraction of the polarized strange distribution from their kaon SIDIS data [79], which are not yet included in this work. These are consistent with small strangeness polarization down to below x ∼ 0.01, but also allow a significantly negative ∆s at x ∼ 0.005. Furthermore, an extraction of the integral of ∆s over the range 0.02 ≤ x ≤ 0.6 by the HERMES collaboration [80] yields 0.037 ± 0.019 ± 0.027, consistent with our result. We stress that beyond the "data-driven" uncertainties that we find for the polarized strangeness distribution, there could well be effects that are outside the leadingtwist framework we are using here and that may have a significant impact on the extracted ∆s. Given the generally low hadron multiplicities in kaon events in the present SIDIS measurements, it is not ruled out that the kaon SIDIS data are affected by higher twist contributions and not suited for an extraction of leading twist strangeness distributions. We note that the information on the parton-to-kaon fragmentation functions in [37] also primarily comes from unpolarized kaon SIDIS data and would not be reliable in the latter case either. A recent determination of the unpolarized strange distribution in the nucleon by HERMES from their SIDIS multiplicities shows an unexpected shape of the distribution [80]. SIDIS measurements at smaller x, as well as at presently available x, but higher Q 2 , will likely be vital for clarifying these issues. These would become available at an electron-ion collider [81].
As can be seen from Fig. 8, the effects due to SU(2) and SU(3) flavor symmetry breaking in usage of the baryon semi-leptonic decay data, see Eqs. (29), (30), have only a very limited impact on the truncated moment of ∆s. This, however, changes dramatically when the full first moment of ∆s is considered, i.e., the contribution to its integral from x < 0.001. This region is presently not constrained by any DIS or SIDIS data, but we remind the reader that the breaking parameters ε SU(2,3) come out very small, see Eq. (34), in our analysis, as a result of the relatively small nominal uncertainty in the F, D values, as we discussed in Sec. III A. This implies that the strange sea distribution will have a large and negative total first moment, ∆Σ s = ∆s 1 + ∆s 1 = −0.114 as seen from Tab. IV, which in turn can only occur if the distribution shows a sign change to negative values at small x, visible in Fig. 3.
It will clearly be important in the future to better understand the strange contribution to nucleon spin structure. If the full first moment ∆Σ s is small, SU(3) symmetry in relating hyperon β-decays to nucleon spin structure would have to be broken at the 40% level or so, which is not ruled out [67,68]. If, on the other hand, SU(3) symmetry is not broken significantly, the implication is that either ∆s turns large and negative at small x, as in our fit, or that the present kaon SIDIS data do not allow a reliable extraction of ∆s(x). On the theoretical side, there have been very recent lattice determinations of the integral ∆Σ s [82], which point to small values. Models of nucleon structure, on the other hand, have led to quite varied predictions for the integral of ∆s, ranging from small to large negative values [83]. We note that the "valence scenario" of [31] has a first moment of the polarized strange distribution very close to zero, which is consequently at the expense of significant violation of the SU(3) flavor symmetry relation in Eq. (30). We finally stress that the size of ∆s is not a topic of interest just for nucleon spin structure enthusiasts: as was pointed out recently [84], the uncertainty in ∆Σ s provides the single largest uncertainty in predictions of the spin-dependent elastic scattering cross sections of supersymmetric dark matter particles on protons and neutrons. Total quark and anti-quark spin contribution ∆Σ: In Fig. 9 we show the χ 2 profile associated with variations of the truncated moment of the quark singlet distribution, ∆Σ  come from the DIS and SIDIS data. The value for the truncated first moment obtained in the best fit is significantly higher than that for the full first moment given in Tab. IV, which is a manifestation of the large negative contribution from strange quarks and anti-quarks that arises in our fit at small x. Thus, keeping in mind the discussion about strangeness above, we conclude that if SU(3) flavor symmetry in relating hyperon β-decays to nucleon spin structure is strongly broken, ∆Σ would be as large as ∼ 0.36 or so, whereas it will be about 30% smaller if SU(3) holds well and the first strange moment ∆Σ s turns out to be large and negative. We note that such lower values of ∆Σ ∼ 0.24 or so have usually been obtained in previous analyses relying on the use of SU (3) symmetry [31,32,33,34,36]. In any case, ∆Σ is certainly much smaller than the typical expectation of ∆Σ 0.6 in quark models. Spin-dependent gluon distribution ∆g: We have already noted in our DSSV paper [28] that the polarized gluon distribution ∆g(x, Q 2 ) comes out rather small in the presently accessed range of momentum fraction x, and prefers to have a node. At variance with the findings of Ref. [34], we do not find any non-overlapping bestfit solutions with gluon polarizations of opposite signs. This duplicity is readily excluded by the RHIC pp data. The RHIC data in fact turn out to play a crucial role in constraining ∆g [28]. The result is shown again in Fig. 3. We do not repeat the plot of the χ 2 profile as a function of the truncated first moment of ∆g here, which may be found in [28]. As can be seen from Tab. III, the integral of ∆g over the RHIC x-region 0.05 to 0.2, ∆g RHIC , is found to be almost zero, while Tab. IV shows that extrapolation to all x results in the gluon spin contribution ∆g 1 = −0.084 at Q 2 = 10 GeV 2 . We stress, FIG. 10: Comparison of ∆g/g for our best fit, at two representative Q 2 , to the extracted ∆g/g from photon-gluon fusion processes investigated by SMC [19], HERMES [18], and COMPASS [20,21]. These data were not included in our global analysis since a consistent NLO framework is not available at present.
however, that this result is not yet reliable due to the large uncertainty in extrapolation to x → 0. In any case, there are presently no indications of a sizable contribution of gluon spins to the proton spin. This is in line with recent theoretical expectations obtained within an effective low-energy theory of broken scale invariance of QCD [85]. Recent bag model estimates also point to relatively modest (but positive) values [86]. Very large values of the integral of the spin-dependent gluon distribution, ∆g 1 ∼ 1.5 or so at Q 2 = 1 GeV 2 , as predicted based on considerations of the QCD axial anomaly [87], become increasingly disfavored, unless ∆g would show a steep rise at small x. Future data from RHIC for spin asymmetries in forward production of correlated hadron or jet pairs, and from running at 500 GeV c.m.s. energy, are expected to shed light on ∆g at lower momentum fractions [88]. Again, also a polarized electron-ion collider [81] would be ideally suited to address this important question and to quantify the amount of gluon polarization at small x from measurements of scaling violations of the structure function g 1 . Other promising channels are, for instance, the polarized photoproduction of single-inclusive hadrons [89] or jets [90].
We have shown the comparison to some of the RHIC data in Fig. 4 (see also Ref. [28]). A way to access ∆g in lepton-nucleon scattering is to measure final states that dominantly select the photon-gluon fusion process, heavy-flavor production, ℓp → hX, and ℓp → h + h − X, where the hadrons have large transverse momentum. Figure 10 shows the corresponding results for the ex-tracted ∆g/g from SMC, HERMES, and COMPASS [18,19,20,21], which have not yet been included in our global analysis. We also show in the figure our result, for two representative Q 2 scales. It should be noted that this comparison is not quite consistent, as the extraction of ∆g/g by the experiments was performed at LO level based on Monte-Carlo generators. Nonetheless, a small ∆g at x ≃ 0.08 − 0.2 as found in our analysis is also well consistent with the data from lepton-nucleon scattering. We expect that the data for the measured spin asymmetries will be included in our global analysis in the future, after the NLO framework for them has been fully developed and been compared to data for the corresponding spin-averaged cross sections.

D. Exploring the fit parameter space
In this Section we briefly present a few more details of the behavior of our total χ 2 near its minimum, which has ramifications, in particular, for the use of the Hessian matrix method for estimating uncertainties. As we noted before, an advantage of the Hessian technique is that it allows to produce sets of "eigenvector PDFs" [43], which in turn can be straightforwardly used in computations of other observables, in order to estimate their PDF uncertainty based on Eq. (11). For this, however, it is very important to know the range of validity of the method, i.e., to which degree χ 2 is parabolic around its minimum. As described in Sec. II B, the first step in the Hessian method is to transform the fit parameters {a i } to a new set {z i } such that surfaces of constant χ 2 turn into hyperspheres in {z i } space, see Eqs. (8), (9) [43]. Figure 11 shows the overlap of each of the original fit parameters {a i } with the eigenvector directions {z i }; the larger the box size the larger the contribution of a certain eigenvector direction to a fit parameter a i . The z i are ordered in terms of the eigenvalues of the Hessian matrix: z 1 corresponds to the largest eigenvalue, i.e., a direction in parameter space where χ 2 changes rapidly, whereas z 19 is only very weakly constrained by data. One can see that in many cases there is a fairly strong correlation between a given original fit parameter a i and a single eigenvector direction z i . The parameters which appear to be constrained best by current data are the normalizations of the sea quarks, Nū, Nd, and Ns, ε SU (2,3) controlling the breaking of SU(2,3) symmetry, and α u+ū , α d+d related to the small x behavior of ∆u+∆ū, ∆d+∆d. Parameters determining the gluon distribution, N g , α g , and η g are less well constrained and mainly correlated with eigenvector directions z 7 to z 12 . ηd, γ d+d , η d+d , and ηs receive contributions from eigenvector directions which are only weakly constrained by data. As we shall see below, this general picture agrees rather well with results for the χ 2 profiles for each fit parameter obtained with the Lagrange multiplier method.
In Figure 12 we investigate the behavior of χ 2 around its minimum, making use of the transformed parameters {z i }. We vary one of the parameters z i at a time, keeping all others fixed. Of course, since each z i has in principle overlap with all fit parameters {a i }, the latter all vary in this procedure. The variation is done in such a way that a given change of ∆χ 2 = T is produced. For truly quadratic behavior near the minimum, as is the underlying assumption in the Hessian approach, the quantity T 2 − ∆χ 2 i , where ∆χ 2 i is the change in χ 2 contributed by the parameter z i that is varied, is trivially zero. This can be compared to the actual dependence of χ 2 on the varied parameter, making no use of the quadratic expansion in (6). Any deviation of T 2 − ∆χ 2 i from zero will signal a departure from the quadratic behavior near the minimum. One can see from the figure that a choice ∆χ 2 = 1 works reasonably well overall, in the sense that overall only fairly small deviations from zero occur. This implies that the Hessian matrix method is reliable for ∆χ 2 = 1 and our eigenvector sets S ± k will produce faithful uncertainty estimates. Some eigenvector directions starting from z 12 and higher do show a certain departure from the ideal behavior even for ∆χ 2 < 1. This is most pronounced for z 17 to z 19 which are the least constrained parameters. In general we have found that the Hessian method breaks down rapidly once one goes beyond ∆χ 2 = 1. Therefore we cannot provide eigenvector sets S ± k corresponding to the more conservative error estimate ∆χ 2 /χ 2 = 2% preferred in [28]. Figure 13 shows the χ 2 profiles including the individual contributions from the DIS, SIDIS, and RHIC pp data sets and from the F , D values for the fit parameters {a i }, obtained with the Lagrange multiplier approach. Clearly, while for some of the parameters the profiles are smooth and parabolic as expected in the simplest approach, for others they are not, showing not only non-parabolic be- havior but variously asymmetric shapes, multiple minima or almost flat regions. It is worth pointing out that these behaviors are not related to a lack of flexibility of the input parameterizations, but to features of the data itself. For example, the double minima observed for N d and η d are associated with two possible "best-fit solutions" to the pion sidis asymmetries, which show strong fluctuations.
In most cases, the behavior is still reasonably quadratic within ∆χ 2 < 1, however, which further justifies the applicability of the Hessian method for ∆χ 2 = 1. Beyond that, simple extrapolation based on an assumed quadratic behavior may give misleading results. We recall that the central values for the parameters can be found in Tab. II.
Reducing the number of parameters of the fit would improve the constraints on the remaining ones, however at the expense of reducing the quality of the fit. The resulting constraint in that case would be strongly dependent on the functional form assumed for the PDFs. In this sense, the robust error analysis based on Lagrange multipliers allows to use more flexible functional forms.

IV. APPLICATIONS OF OUR UNCERTAINTY ESTIMATES
The 38 eigenvector PDF sets S ± k that we have constructed for the Hessian matrix method are ideally suited for estimating the PDF uncertainty of other observables. In this way, one can for example gauge the accuracy that future additional measurements will need to have in order to have a significant impact on our knowledge of the PDFs. In this Section, we will present a few examples for the case of RHIC. In view of the fact that RHIC has just completed its first physics run at √ S = 500 GeV, we will focus on predictions for this c.m.s. energy. Figure 14 shows the NLO double-spin asymmetries A LL defined in Eq. (36) for pp → hX (h = π 0 , π ± ) and pp → jetX, for our central DSSV fit (solid lines), including the Hessian uncertainty bands for ∆χ 2 = 1 using Eq. (11). One can see that the asymmetry for π 0 remains very small until about p T ∼ 20 GeV, as could be expected from a simple scaling of the asymmetry A LL at √ S = 200 GeV shown in Fig. (4) with x T ≡ 2p T / √ S. It then rapidly increases. The asymmetry for negatively charged pions remains small for all transverse momenta and in fact turns slightly negative at high p T . In contrast, A π + LL is higher than A π 0 LL , reaching about 3% at the highest p T shown. The behavior of the various pion asymmetries is closely tied to that of the polarized gluon distribution: at high p T , relatively large values of x are relevant, where our ∆g is positive, and quark-gluon scattering dominates. An important contribution to the spindependent cross section thus involves the combination (∆u ⊗ D π u + ∆ū ⊗ D π u + ∆d ⊗ D π d + ∆d ⊗ D π d ) ⊗ ∆g of parton distributions and fragmentation functions. For π + production, the u quark andd anti-quark contributions are expected to dominate, as these are valence quarks in a π + . The combination ∆u + ∆d is positive as Fig. 3 shows. The large negative contribution associated with ∆d is suppressed here. For π 0 production, the participating fragmentation functions are all equal, and one probes the sum of up and down quark and anti-quark distributions, which is positive but smaller than ∆u + ∆d. Finally, for π − production, the main contribution involves ∆d+∆ū, which explains the downturn of A π − LL to negative values at high p T . Clearly, the three pion asymmetries are also sensitive to the sign of ∆g. We note that preliminary results for A π ± LL at √ S = 200 GeV have recently been reported from RHIC [27].
Similar features as for the π 0 asymmetry are observed for jets. Very roughly, one finds that A jet LL (p T ) ≈ A π 0 LL (kp T ), where k ≈ 0.5 or so, corresponding to the fact that on average only the fraction k of the total jet momentum is taken by an observed π 0 . This implies that, at a given p T , the jet spin asymmetry is smaller than that for π 0 .
We have seen in the previous Section that the SIDIS data have given some first insights into the flavor structure of the polarized sea distributions of the nucleon. On the other hand, the uncertainties in SIDIS are still quite large, and it is in particular difficult to quantify the systematic uncertainty of the results related to the fragmentation mechanism at the relatively modest energies available so far. Complementary and clean information on ∆u, ∆ū, ∆d, and ∆d will come from pp → W ± X at RHIC, where one will exploit the maximally parityviolating couplings of produced W bosons to left-handed quarks and right-handed anti-quarks [88,91]. The high scale set by the W boson mass makes it possible to extract quark and anti-quark polarizations from inclusive lepton single-spin asymmetries in W boson production with minimal theoretical uncertainties, as higher order and sub-leading terms in the perturbative QCD expansion are suppressed [92,93,94,95].
As a further application of our Hessian uncertainty PDFs, we show in Fig. 15 the single-longitudinal spin asymmetries, for the processes pp → ℓ ± X, where the arrow denotes a longitudinally polarized proton and ℓ = e or µ. The FIG. 14: Double-spin asymmetries ALL for jet and pion production at RHIC at √ S = 500 GeV as functions of the transverse momentum pT of the produced final state. We show the results for the best-fit parton distributions from our global analysis, along with the uncertainties estimated using the Hessian method, allowing changes of one unit in χ 2 . We also show the results for the "standard" scenario of [31] (dashed lines; in the lower plot the result for π + (π − ) is given by the top (bottom) curve). We have used the CTEQ6M unpolarized parton distributions [46] for the calculation of the denominator of the asymmetry. For the pions, we have assumed pseudo-rapidity coverage of |η| < 0.35, and for the jets of −1 ≤ η ≤ 2. charged lepton is assumed to have been produced by a leptonic decay of the W ± boson. The asymmetries are shown as functions of the charged lepton's rapidity η lept , with η lept counted positive in the forward direction of the polarized proton. We have integrated over p T > 20 GeV, where p T is the lepton transverse momentum. The results shown in the figure are based on a simple LO calculation of the processes qq ′ → W ± → ℓ ± ν; the NLO corrections which we should in principle include for consistency are negligible for this observable [92,93,94,95].
For W − production, neglecting all partonic processes but the dominantūd → W − one, the spin-dependent cross section in the numerator of the asymmetry is found to be proportional to the combination [93] ∆ū(x 1 )d(x 2 )(1−cos θ) 2 −∆d(x 1 )ū(x 2 )(1+cos θ) 2 , (39) where θ is the polar angle of the electron in the par- Single-longitudinal spin asymmetries for chargedlepton production at RHIC through production and decay of W bosons. The bands correspond to our uncertainty estimates based on the Hessian ∆χ 2 = 1 eigenvector PDFs. We also show in the figure the spin asymmetries obtained for the "standard" and "valence" scenarios of [31]. For the case of W − , we also show the results of two fits for which the ratio ∆d(x)/d(x) is forced to turn to +1 as x → 1, see text.
tonic c.m.s., with θ = 0 in the forward direction of the polarized parton. At large negative η lept , one has x 2 ≫ x 1 and θ ≫ π/2. In this case, the first term in Eq. (39) strongly dominates, since the combination of parton distributions, ∆ū(x 1 )d(x 2 ), and the angular factor, (1 − cos θ) 2 , each dominate over their counterpart in the second term. Since the denominator of A L is proportional toū(x 1 )d(x 2 )(1 − cos θ) 2 + d(x 1 )ū(x 2 )(1 + cos θ) 2 , the asymmetry provides a clean probe of ∆ū(x 1 )/ū(x 1 ) at medium values of x 1 . Indeed, the Hessian uncertainty band for ∆ū shown in Fig. (3) is directly reflected in the band we show in Fig. 15. We also show in the figure the spin asymmetries obtained for the "standard" and "valence" scenarios of [31]. The latter has a large and positive ∆ū distribution at the relevant x ∼ 0.1, which clearly shows in the asymmetry. By similar reasoning, at forward rapidity η lept ≫ 0 the second term in Eq. (39) dominates, giving access to −∆d(x 1 )/d(x 1 ) at relatively high x 1 . We have discussed in Subsec. III C that there is interest in the question if the polarized down-quark distribution turns positive in the large-x region, for which there are currently no indications, see Fig. (5). As Fig. 15 shows, the asymmetry for W − production becomes large and positive at high η lept , which precisely reflects the fact that ∆d(x) remains negative at high x in our DSSV fit. It is interesting to investigate how the asymmetry might look if ∆d/d were to turn to +1 as x → 1. In order to do this, we have produced two fits where ∆d/d is forced to have this behavior. The two fits are characterized by the value x 0 where ∆d(x, M 2 W ) changes sign from negative to positive values. We have chosen x 0 = 0.67 and x 0 = 0.55. The χ 2 values for these two fits are of course significantly worse than for our DSSV best fit, by about four units for x 0 = 0.67 and about 25 units for x 0 = 0.55. The results for the two fits are shown by the dotted lines in Fig. 15. It should be well possible at RHIC to measure the asymmetry for values η lept out to 2 [88]. We note that the behavior of ∆d/d at high x will also be further addressed by experiments at Jefferson Lab after the 12 GeV upgrade [96].
For W + production, one has the following structure of the spin-dependent cross section [93]: ∆d(x 1 )u(x 2 )(1+cos θ) 2 −∆u(x 1 )d(x 2 )(1−cos θ) 2 . (40) Here the distinction of the two contributions by considering large negative or positive lepton rapidities is less clear-cut than in the case of W − . For example, at negative η lept the partonic combinationd(x 1 )u(x 2 ) will dominate, but at the same time θ ≫ π/2 so that the angular factor (1 + cos θ) 2 is small. Likewise, at positive η lept the dominant partonic combination ∆u(x 1 )d(x 2 ) is suppressed by the angular factor. So both terms in Eq. (40) will compete essentially for all η lept of interest. This is reflected in the behavior of the calculated spin asymmetry A W + L shown in Fig. 15, which does not show as clear features as the one for W − bosons. Nonetheless, the W + measurements at RHIC will of course still be of great value. In fact, our global analysis technique is precisely suited for extracting information on the polarized PDFs even if there is no single dominant partonic subprocess.

V. CONCLUSIONS
We have presented details of a recent study of the helicity parton distribution functions of the nucleon, which used experimental information available from inclusive and semi-inclusive polarized deep-inelastic leptonnucleon scattering and from polarized proton-proton scattering at RHIC. The data sets were used jointly in a next-to-leading order global QCD analysis, which allows to extract the set of parton distributions that provides the optimal overall description of the data, along with estimates of its uncertainties. We have presented techniques and computational methods that speed up the next-to-leading order calculations for pp scattering to the level required in practice for a global analysis. Our technique is formulated in Mellin moment space. A key feature is that the computationally most challenging parts are done only once, prior to the fit. Use of a Monte-Carlo sampling method allows us to perform this one-time calculation very efficiently.
Our extracted parton distributions show particularly interesting features in the sea quark and gluon sector. We find evidence for a mostly positive ∆ū and a negative ∆d distribution, so that ∆ū−∆d is positive. This behavior has been predicted by a number of models of nucleon structure. The polarized strange quark distribution ∆s comes out slightly positive at medium x, which is driven by the semi-inclusive kaon DIS data and could be subject to rather large systematic uncertainties. ∆s turns negative at x 0.02 as a result of constraints from SU(3) symmetry, which have a relatively small nominal error. If true, this means that ∆s acquires its large negative integral essentially completely from the small-x region. As a further consequence, quark and anti-quark spins combined contribute about a fourth to a third of the proton's spin, with the lower value arising if strange quarks and anti-quarks are indeed strongly negatively polarized at low x. Finally, we have found that the gluon helicity distribution ∆g(x, Q 2 ) is small in the region of momentum fraction accessed directly so far by RHIC, with likely a node and an almost vanishing integral over that region. Reliable statements about the full gluon spin contribution to the proton spin are presently not yet possible.
We have performed uncertainty estimates for our polarized parton distributions, using both the Lagrange multiplier technique and the improved Hessian approach. To obtain these, a large number of additional fits are necessary, for which the computational techniques we have developed are particularly important. We find that both approaches yield consistent results for moderate departures from the best fit, typically ∆χ 2 = 1. For larger ∆χ 2 , significant differences develop as a result of departures from parabolic behavior of χ 2 around its minimum. This implies that the Hessian matrix method becomes unreliable. We have produced a set of 38 "eigenvector" parton distributions for the Hessian method with ∆χ 2 = 1 [97], which may be used to estimate the uncertainty of any observable that depends on the distri-butions. We stress, however, that we presently prefer a more conservative choice of ∆χ 2 /χ 2 = 2% as a tolerance criterion for acceptable parton distributions. Unfortunately, the behavior of χ 2 around its minimum does not warrant use of the Hessian method for producing eigenvector parton distributions in this case.
We have used the ∆χ 2 = 1 eigenvector distributions to obtain predictions for spin asymmetries for high transverse momentum pion and jet production in polarized proton-proton collisions at 500 GeV center-of-mass energy at BNL-RHIC, as well as for W boson production. The former would give information on ∆g at lower x, while the latter would provide a clean new probe of the polarized quark and anti-quark distributions, which is important in view of the uncertainties inherent in semiinclusive DIS. Our results indicate that there is significant potential for RHIC to provide further important insights into nucleon helicity structure. It will be straightforward to include all the forthcoming data in the global analysis.