First Monte Carlo global QCD analysis of pion parton distributions

We perform the first global QCD analysis of parton distribution functions (PDFs) in the pion, combining $\pi A$ Drell-Yan data with leading neutron electroproduction from HERA within a Monte Carlo approach based on nested sampling. Inclusion of the HERA data allows the pion PDFs to be determined down to much lower values of $x$, with relatively weak model dependence from uncertainties in the chiral splitting function. The combined analysis reveals that gluons carry a significantly higher pion momentum fraction, $\sim 30\%$, than that inferred from Drell-Yan data alone, with sea quarks carrying a somewhat smaller fraction, $\sim 15\%$, at the input scale. Within the same effective theory framework, the chiral splitting function and pion PDFs can be used to describe the $\bar d-\bar u$ asymmetry in the proton.

As the lightest QCD bound state, the pion has historically played a central role in the study of the strong nuclear interactions. On one hand, it has been the critical ingredient for understanding the consequences of dynamical chiral symmetry breaking in QCD, and how this dictates the nature of hadronic interactions at low energies. On the other hand, its quark and gluon (or parton) substructure has been revealed through high energy scattering experiments, such as Drell-Yan (DY) lepton-pair creation in inclusive pion-nucleon scattering [1]. In some cases, both aspects are on display, as in the role of the pion cloud of the proton in generating a flavor asymmetry in its light antiquark sea,d =ū [2].
As the simplest qq state, the structure of the pion is relatively more straightforward to compute theoretically than baryons, but the absence of fixed pion targets has made it difficult to determine the pion's structure experimentally. Most information on the partonic structure of pions has come from pion-nucleus scattering with prompt photon or dilepton production at CERN [3,4] and Fermilab [5]. These data have been used in several global QCD analyses [6][7][8][9] to fit the momentum dependence of the pion's parton distribution functions (PDFs) for parton fractions x π 0.1 of the pion's light-cone momentum.
While the DY data constrain reasonably well the pion's valence PDFs, the sea quark and gluon PDFs at small x π values have remained essentially unknown. More recently, leading neutron (LN) electroproduction from HERA [10,11], which at forward angles is expected to be dominated by pion exchange, has been used to study the pion structure function down to very low values of x π ∼ 10 −3 . The interpretability of the LN data in terms of pion PDFs is limited, however, by the model dependence inherent in this process, in which the cross section is given as a product of a proton to neutron + pion "chiral splitting function" and the structure function of the (nearly on-shell) exchanged pion. Consequently the LN data have never been used in global QCD analyses, although recently the first steps toward their inclusion were taken by McKenney et al. [12], who studied the impact of the model dependence on the extracted pion structure function by constraining the p → nπ + splitting function empirically.
In the present work, we combine the strategy of global QCD analysis with an empirical approach to using the DY and LN data in the same fit to determine the pion PDFs in both the high-x π and low-x π regions. We use for the first time a Monte Carlo (MC) approach, based on the nested sampling algorithm [13][14][15], to perform the global analysis at next-to-leading order (NLO) in the strong coupling. In contrast to the previous single-fit analyses based on maximum likelihood methods [6][7][8][9], the MC approach allows a systematic exploration of the parameter space by computing the likelihood function directly, providing a rigorous determination of PDF uncertainties.
An important feature of our analysis is the ability of the MC fit to quantify the uncertainty on the dependence of the extracted pion PDFs on the chiral splitting function model. To test the robustness of the chiral framework, we also perform a simultaneous fit to the pion DY + LN data together with the E866 pd/pp DY data [16], from which thed/ū ratio was extracted, using the same MC methodology. Such an analysis provides the most comprehensive study of pion PDFs and their impact on different observables.
In the pion-induced DY process [1], partons from the pion and target nucleus A annihilate to produce a dimuon pair in the final state, πA → µ + µ − X, with cross section where f π i (f A j ) is the PDF for parton flavor i in the pion arXiv:1804.01965v1 [hep-ph] 5 Apr 2018 (flavor j in the nucleus) as a function of the parton momentum fractionx π (x A ), and C ij is the hard scattering kernel [17,18], with µ the renormalization scale. The cross section is differential in the dilepton invariant mass squared Q 2 and rapidity Y , in terms of which one defines x π,A = √ τ e ±Y , where τ = Q 2 /S and S is the π-target invariant mass squared. At lowest order the parton momentum fractions are given byx π,A = x π,A . Typically, the experimental DY cross sections are analyzed in terms of the Feynman variable x F = x π −x N , where x N = Ax A is the nuclear Bjorken variable scaled per nucleon [19]. The available pion DY data from the CERN NA10 [4] and Fermilab E615 [5] experiments were taken on tungsten nuclei, and for the nuclear PDFs in our analysis we use the parametrization from Eskola et al. [20]. For the LN production process, ep → enX, the pion PDFs enter indirectly, under the assumption that the charge exchange cross section at low values of t and large neutron longitudinal momentum fractions x L is dominated by single pion exchange. The differential LN cross section Here f πN (x L ) is the chiral splitting function for fraction x L ≡ 1−x L = x/x π of the proton's light-cone momentum carried by the pion, and F π 2 (x π , Q 2 ) is the pion structure function, evaluated at NLO. The splitting function is evaluated from chiral effective theory [21][22][23], and for x L > 0 is given by ]/x L , with M and m π the nucleon and pion masses, g A the axial charge, and f π the pion decay constant. The form of the splitting function in Eq. (3) is constrained by chiral symmetry [21,24,25], and its infrared or leading nonanalytic behavior is model independent [26][27][28][29]. The ultraviolet behavior, however, is dependent on the regularization procedure, represented in Eq. (3) by the function F. In the literature various forms have been advocated, including cutoff regularization, Pauli-Villars, and phenomenological πN form factors, and following Ref. [12] we consider several forms, with α π ≈ 1 GeV −2 , and Λ is a cutoff parameter. We also considered a model based on a large-k ⊥ cutoff [23], and the Bishari model [11] (which is analogous to the Regge form but with Λ → ∞), but found poor agreement between the fitted results and data, and do not consider them further.
In addition to the LN structure function data collected by H1 [10], the ZEUS collaboration measured the ratio [11] r(x, of LN to inclusive cross sections, where the latter is expressed in terms of the proton structure function, F p 2 , and ∆x L is the bin size in x L . Consistent with expectations from earlier theoretical calculations [30,31], at large x L ∼ 1 pion exchange is the dominant contribution [12]. Other processes, such as absorption and the exchange of other mesons, play an increasingly important role at smaller x L . Instead of choosing a specific value of x L below which one pion exchange is assumed, we fit the lower value of x L for which the data can be described within this framework. For the data analysis we use a Bayesian Monte Carlo method based on the nested sampling algorithm [13][14][15], which allows a faithful Monte Carlo representation of the probability distribution P(a|data) = L(data|a)π(a)/Z, where a is an n-dimensional array of the pion PDFs shape parameters. Here π(a) is the Bayesian prior distribution for a, which allows the parameter sampling to be restricted to physical regions, L(data|a) = exp[− 1 2 χ 2 (a)] is the likelihood function, and Z = d n a L(data|a) π(a) is the Bayesian evidence parameter, which normalizes the probability distribution. We use a χ 2 function in the likelihood that takes into account correlated systematic shifts, as well as overall normalizations of the data sets [32]. For physical observables O, such as the pion PDFs and functions thereof, from the MC samples {a k } one then obtains expectation values E where {w k } are the MC weights. Similar MC technology based on Bayesian statistics has also been applied recently to study nucleon PDFs [33,34] and fragmentation functions [35], as well as the transverse momentum dependent transversity distribution [36].
For the pion valence PDFs we assume charge symme- , and invoke SU(3) symmetry for the pion sea, q π s ≡ū π + = d π + = s π + =s π + . The valence, sea quark, and gluon PDFs are parameterized at the input scale of the charm quark mass Q 2 0 = m 2 c = (1.3 GeV) 2 by the form where a = {N, α, β} are the fitting parameters and B is the Euler beta function. The valence PDFs are normal-ized such that 1 0 dx π q π v = 1, and the momentum sum rule gives the constraint 1 0 dx π x π (2q π v + 6q π s + g π ) = 1. The fits to the DY and LN data sets are shown in Fig. 1, where for clarity the E615 and HERA points are scaled by 3 i . To avoid the J/Ψ and Υ resonances, the DY data were restricted to the mass region 4.16 < Q < 8.34 GeV, covering the range 0.05 ≤ x F ≤ 0.9. Generally very good agreement is found for the entire set of 250 data points. For the best fit, corresponding to model (i) for the LN cross section with a cutoff Λ = 1.31(4) GeV, the combined χ 2 /N dat is 0.98 (244.8/250). The overall normalizations for the DY data are found to be 0.816, 0.758 and 0.985 for the NA10 (194 GeV), NA10 (252 GeV) and E615 data sets, and 1.17 and 0.964 for the H1 and ZEUS LN data, respectively.
For the LN data good fits were obtained for the cut x L > 0.8; including smaller-x L data deteriorated the fit due to larger non-pionic contributions away from the forward limit [12,30,31]. Fitting the DY data alone yields only marginally smaller χ 2 values, with χ 2 /N dat = 0.97 (55.5/70 for NA10 and 82.6/72 for E615). For the combined DY and LN data sets, the total χ 2 per datum for other models are also close to The resulting pion PDFs are shown in Fig. 2 for the valence, sea quark, and gluon distributions at Q 2 = 10 GeV 2 . Compared with the DY-only fits, which constrain mainly the valence quark PDF, the simultaneous DY+LN analysis yields significantly reduced uncertainties on the pion sea and gluon distributions at low x π . While the addition of the LN data gives a small, ≈ 10% reduction of the valence PDF at intermediate x π , the impact on the sea is more dramatic, with the gluon PDF increasing twofold at x π ∼ 0.001 − 0.1 compared with the DY-only result, but with half of the uncertainty. The sea quark PDF q π s is reduced at x π 0.1, but is slightly larger at x π 0.1 for the full result. Importantly, the model dependence of the combined DY+LN fit (represented in Fig. 2 by the yellow bands) reveals a relatively small uncertainty, especially compared with the scale of the effect induced by the addition of the LN results.
For the valence PDF in the large-x π region, our analysis finds a behavior ∼ (1 − x) at the input scale, consistent with Nambu-Jona-Lasinio [37] and constituent quark models [38] and duality arguments [39], but harder than expectations based on pQCD [40] and Dyson-Schwinger calculations [41], which prefer a (1 − x) 2 fall-off. On the other hand, the present analysis does not include threshold resummation effects, which are known to be important at large x π [9,42], which will be examined in a separate analysis [43].
The inclusion of the LN data into the global analysis allows a more precise breakdown of the pion momentum  into fractions carried by valence quarks, sea quarks, and gluons, shown in Fig. 3. The total valence, sea and gluon momentum fractions for the full analysis at the input scale Q 2 = m 2 c are found to be { x π v , x π s , x π g } = {54(1)%, 16(2)%, 30(2)%} for the best fit with model (i).  fraction carried by gluons is about 3 times larger, but with less than half of the uncertainty. Since the valence fraction remains relatively unchanged, the momentum sum rule forces the sea quarks to carry about 1/2 of the momentum fraction compared with the DY-only fit.
This turns out to be similar to the result from the older SMRS analysis [6] of DY data, which considered several scenarios in which the momentum fraction carried by sea quarks at Q 2 = 5 GeV 2 varied from 10% to 20%. With the valence momentum fraction in [6] constrained to be 46%, the momentum fraction carried by gluons varied from 43% to 34% over this range. At the same scale the analogous fractions from our analysis are {48(1)%, 17(1)%, 35(2)%}, which is closer to the SMRS scenario with maximum sea and minimum glue. In contrast, the GRS analysis [8], which assumes a constituent quark model to constrain the sea quark and gluon PDFs, gives a gluon momentum fraction (44%) that is actually larger than the valence fraction (39%) at this scale.
As an application of our results and a test of the validity of the chiral framework for the LN data, we consider the contribution from the p → n π + dissociation to thē d −ū asymmetry in the proton sea. This asymmetry, predicted from pion loop effects in the nucleon [2], was conclusively established by the pp and pd Drell-Yan data from the E866 experiment at Fermilab [16]. While the effect of these data on PDFs is most rigorously quantified through global QCD fits, an approximate NLO analysis of the E866 data extracted thed −ū difference, shown in Fig. 4, assuming knowledge of the valence quark distributions in the proton. The antiquark asymmetry can be represented within chiral effective theory as a convolution (represented by the symbol "⊗") of the p → n π + and p → ∆ ++ π − splitting functions and the valence pion PDF [44,45], The π∆ splitting function, which does not contribute to the LN cross section but enhances theū PDF overd for inclusive processes, is given in Ref. [23]. Refitting the DY and LN data together with the 15 additional E866 extracted data points, and attributing the entire asymmetry to pion loops, gives a best fit for model (i) with only a marginally larger χ 2 /N dat = 1.03 (272.4/265). The combined fit has additional sensitivity to the valence quark pion PDF, however, the resulting PDFs are relatively stable, and the fitted model (i) cutoff parameter, Λ = 1.35(2) GeV, is consistent with that from the fit without E866 data.
The resultingd−ū asymmetry is shown in Fig. 4, for all models (i)-(v) for the regulator F. Without compromising the description of the DY or LN data, a reasonably good fit to the E866 data points can be achieved in all cases for x 0.2, beyond which all the fits overestimate the data. It is known that the apparent change of sign ind −ū at high x is difficult to accommodate theoretically [46], and the new DY SeaQuest experiment [47] will allow a more precise determination of the asymmetry up to x ≈ 0.5.
In the future, pion PDFs will be further constrained by new πA DY data from COMPASS [48,49], as well as from the Tagged DIS (TDIS) [50] experiment at Jefferson Lab, which will study pion structure through the charge exchange mechanism in leading proton production form a quasi-free neutron in the deuteron, ed → eppX. Theoretically, effects from gluon resummation [9,42] and higher twists [40] will be explored [43] systematically in order to unravel the behavior of pion PDFs at very high x π ∼ 1. Beyond this, an ultimate future goal will be a simultaneous fit of pion, proton and nuclear PDFs within the same MC global QCD analysis.