Global QCD Analysis of Pion Parton Distributions with Threshold Resummation

We perform the first global QCD analysis of pion valence, sea quark, and gluon distributions within a Bayesian Monte Carlo framework with threshold resummation on Drell-Yan cross sections at next-to-leading log accuracy. Exploring various treatments of resummation, we find that the large-$x$ asymptotics of the valence quark distribution $\sim (1-x)^{\beta_v}$ can differ significantly, with $\beta_v$ ranging from $\approx 1$ to $>2.5$ at the input scale. Regardless of the specific implementation, however, the resummation induced redistribution of the momentum between valence quarks and gluons boosts the total momentum carried by gluons to $\approx 40\%$, increasing the gluon contribution to the pion mass to $\approx 40$ MeV.

Introduction.-As the lightest known hadron, the pion presents itself as a dichotomy of nature. On the one hand, as the pseudo-Goldstone boson associated with chiral symmetry breaking it is fundamental for understanding low-energy hadronic interactions [1,2]; on the other, as a QCD bound state composed of quarks and gluons, its partonic structure reveals itself in high-energy scattering experiments much like other hadrons. Since pions do not exist as free stationary targets from which one could scatter, such as protons or stable nuclei, details of their partonic structure have consequently been more difficult to establish empirically. While experiments with secondary pion beams on nuclear targets [3,4] have provided intriguing glimpses into the pion's valence quark structure, many open questions remain.
In parallel developments, the Jefferson Lab Angular Momentum (JAM) Collaboration recently explored [49] the inclusion of leading neutron (LN) electroproduction data from HERA [50,51], in addition to the DY data, to constrain the valence, sea quark and gluon distributions at low and high x values, using Bayesian Monte Carlo methods. A subsequent fixed order analysis also included high-p T DY data [52], together with p T -integrated data. Other phenomenological analyses have utilized DY and prompt photon data to constrain pion PDFs [48,[53][54][55][56][57][58], and the growing number of recent lattice calculations [59][60][61][62][63][64][65], some including threshold resummation [66], is a testament to the importance of better understanding the pion's PDFs.
In this Letter we bring these strands together to perform a global QCD analysis of DY and LN data within the JAM framework that includes for the first time a systematic study of soft gluon resummation effects on pion PDFs at next-to-leading logarithmic (NLL) accuracy. In particular, we critically examine the universality of the resummation impact on the effective β v parameter, by considering several viable resummation prescriptions, including the traditional Mellin-Fourier (MF) method [45,46,67,68] and the more recently developed double Mellin (DM) method [47].
Threshold resummation of Drell-Yan.-In the inclusive pion-induced DY process [69], a pion beam incident on a nuclear target, with total center of mass energy √ S, produces an inclusive µ + µ − pair with invariant mass Q. The cross section differential with respect to τ = Q 2 /S and rapidity, Y , can be factorized into a convolution of perturbatively calculable hard scattering coefficients C ij arXiv:2108.05822v2 [hep-ph] 30 Nov 2021 and the collinear pion and nuclear PDFs [70], where the PDFs f π(A) i(j) for parton flavor i(j) in the pion (nucleus) are functions of the parton's light-front momentum fraction x π(A) with respect to the parent hadron, and µ is the factorization scale. The integration variables are y = (u−z)/[(1−z)(1+u)], where u = e −2Y (x π /x A ), and z = Q 2 /ŝ, withŝ = x π x A S the partonic invariant mass squared.
At leading order (LO), only the qq channel contributes, and the hard coefficient, C (0) qq = δ(1−z), allows the PDFs to be evaluated at momentum fractions x 0 π,A = √ τ e ±Y . At next-to-leading order (NLO), the tree-level qg channel enters, while the gg and qq channels appear at next-tonext-to-leading order (NNLO). Focusing on the qq channel, the hard coefficients can be schematically organized as a sum of a Born term that describes the LO qq annihilation, a virtual term that includes interference of the one-loop and Born diagrams, and a real emission term describing the radiation in the final state. The qq channel plays the dominant role, and is the focus of this work.
The phase space for real gluon radiation vanishes in the limit z → 1. In the infrared region where the real gluon radiation becomes soft, mass singularities from propagators in the real and virtual diagrams cancel up to terms of the form α k s log 2k−1 (1 − z) /(1 − z), known as threshold logarithms, which appear order by order in perturbation theory. At NLO, for example, in the rapidityintegrated DY cross section the hard coefficient includes terms C Because the PDFs are steeply falling at large x, the partonic cross section in the threshold region plays a substantial role in the overall hadronic cross section [45,47,71,72]. Changes in the perturbative C ij coefficients are compensated by changes in the PDFs to produce the same physical cross section [45]. Near the partonic threshold, the logarithms become increasingly important, and must be resummed in order to maintain the integrity of the perturbative expansion. The threshold resummation framework of Sterman [73] and Catani and Trentadue [74] allows one to systematically include the large logarithmic contributions to the DY process from soft gluon emissions to all orders of α s .
The calculation of threshold resummation is not straightforward in momentum space. In Mellin space, however, the phase space integrals of multiple soft gluon emissions decouple, allowing convenient organization of the threshold logarithms to be resummed. For the rapidity-dependent DY cross section, an additional transformation is needed beyond the single Mellin, which we take as either an additional Fourier transform in Y [67,75] or a double Mellin transform in x 0 π and x 0 A [47], The convolution form of the cross section (1) decouples in conjugate space into a simple product of hard coefficients and PDFs. In the MF case, defining a partonic rapidity Y = Y − 1 2 log(x π /x A ) = − 1 2 log u allows the partonic cross section to be written as a Mellin transform in z and Fourier transform in Y . Note that the threshold region z ≈ 1 maps into large N in Mellin space.
The formal definition [72,74] of the Mellin transform for the resummed hard coefficients implies an evaluation at the Landau pole. Although the Landau pole can be found explicitly in the resummed coefficients as a function of N , there is ambiguity in how it should be treated [46,47,72,76]. The large logarithms being resummed appear at z ≤ 1−1/(N e γ E ), where γ E is the Euler constant, truncating the full Mellin transform and giving rise to analytic terms [45,47]. Since the large Mellin logarithms are far from the Landau pole, N L , the minimal prescription (MP) [72] makes use of special contours C MP that avoid N L , while enclosing poles from the PDFs in the Mellin inversion. To use the Mellin expressions for the threshold resummation we solve DGLAP evolution in Mellin space, which avoids numerical instability of the hard coefficients in an x-space computation [46,47].
To compare with experiment, we invert Eqs.
(2) to momentum space to obtain the differential cross section, Since the resummed expansion is performed at NLL, the threshold logarithms at NLO also appear in the resummed calculation. Expanding the NLL terms in orders of α s and subtracting up to the O(α s ) terms ensures that the NLO pieces are not double counted. In addition to the difference between the MF and DM methods, another ambiguity appears in the MF case. The real emission terms proportional to the threshold logarithms include a factor δ Y ± log(1/ √ z) , which after the Fourier transform produces a cos(M log(1/ √ z)) term. Expanding the cosine reveals O((1 − z) 2 ) corrections, which are subleading near threshold, and the cosine may be approximated as 1 [68]. In this work, we refer to this as the "expansion" method.
Alternatively, the cosine may be kept without expanding [45,67], which we refer to as the "cosine" method. Similar subleading O((1 − z) n ) terms appear in non-qq channels, so one may argue that a consistent treatment requires inclusion of the other channels. On the other hand, because of its dependence on M , the cosine method provides additional information about the rapidity distribution. In our analysis, we systematically explore all three methods, and quantify the dependence of the fitted PDFs on the resummation method choices.
Bayesian inference.-Our numerical analysis uses the Bayesian Monte Carlo methodology developed by the JAM Collaboration [77][78][79][80][81][82][83][84], whereby the probability of the set a of best fit parameters conditioned on the data is proportional to the likelihood of the data given the parameter set. We sample the posterior probability distribution P(a|data) ∝ L(data|a)π(a), where L(data|a) = exp (− 1 2 χ 2 (a, data)) is the likelihood function and π(a) is the prior distribution. Data resampling is used with Gaussian noise added to the data within the quoted uncorrelated uncertainties [77,78], and the posterior distribution is sampled via multiple likelihood regressions. The expectation value E and variance V of an "observable" O (such as a PDF or cross section) are computed from E where the sum is over N rep number of replicas. The pion PDF for a given flavor i is parametrized at the input scale, µ 0 = m c = 1.27 GeV, by the form for the set of parameters a i = {N i , α i , β i , γ i }. We assume charge symmetry for the valence quark PDF, q v ≡ū π − v =ū π − − u π − = d π − v , and a flavor symmetric sea, q s ≡ u π − =d π − = s π =s π . Valence quark number conservation, 1 0 dx q v (x, µ) = 1, and the momentum sum rule, 1 0 dx x 2q v (x, µ) + 6q s (x, µ) + g(x, µ) = 1, constrain the normalizations N v and N s , respectively. We further set γ s,g = 0, as these could not be constrained by existing data, and choose N g > 0 and γ v > −1 to avoid PDFs becoming negative.
The pion valence, sea quark, and gluon PDFs are fitted to the available pion-nucleus Drell-Yan data from E615 [3] (61 data points) and NA10 [4] (56 points), pre- The analysis is limited to the range 4.16 < Q < 7.68 GeV and 0 < x F < 0.9 to avoid the J/ψ and Υ resonances and edges of phase space. For the nuclear PDFs we use the Eskola et al. parametrization [85], although using the nCTEQ fit [86] showed no differences in final results. As in previous analyses [49,52], we include leading neutron electroproduction data at small x from HERA (58 points for H1 [50] and 50 points for ZEUS [51] data), which at very forward angles are expected to be dominated by pion exchange [87,88].
Monte Carlo analysis.-Fairly good fits to the data can be obtained for all resummation prescriptions con- sidered, including no resummation. For the fixed order NLO analysis, a total χ 2 per datum was found of χ 2 dat = 0.81 for the combined DY and LN datasets, consistent with the recent analysis [52], while including resummation with the cosine, expansion, and double Mellin methods gives a total χ 2 dat of 1.29, 0.95, and 0.80, respectively. The χ 2 dat variation comes almost exclusively from the DY datasets, with the low-x LN data mostly unaffected by the perturbative QCD treatment of the DY cross sections.
The ratio of the experimental DY data to the calculated cross sections is shown in Fig. 1 versus x F in bins of √ τ , for a representative subset of the E615 data [3]. The comparison of the fixed order NLO fit with those with resummation shows most variance of the data descriptions in the large-x F region, while at small x F no discernable differences are seen among the prescriptions. The NLO and double Mellin methods give similar descriptions of the data, while the cosine and expansion methods have larger contributions to the theory at high x F , rendering a ratio below unity. Based on the χ 2 dat criteria, we conclude that the NLO, NLL expansion and double Mellin prescriptions give equally good fits, while the cosine method is slightly disfavored.
The fitted valence quark, sea quark, and gluon (scaled by a factor 1/10) PDFs are shown in Fig. 2 for the various resummation prescriptions at the input scale µ = µ 0 . The extracted valence distributions with NLL resummation all clearly show a softer falloff as x → 1, as magnified in the inset. The NLO valence PDF has the hardest distribution, followed by the double Mellin method, while the expansion method is softer, and the cosine method yields the softest falloff. As seen in Fig. 1, the cosine method tends to overpredict the data at large x F , indicating that the contribution to the hard coefficients becomes too large and the PDFs cannot sufficiently adjust because of restrictions from the sum rules and LN data. To compensate for the large hard coefficients, the PDF becomes suppressed to bring the calculated cross section down to match the data. Interestingly, when positivity of the PDFs is not enforced, the cosine method admits slightly negative valence solutions at x 0.85. For the sea quark and gluon PDFs, there is obviously greater spread, with a difference that the double Mellin method admits a sea quark shape that is somewhat larger at high x than the other methods. Additionally, the central values exhibit larger gluon distributions for the cosine and expansion methods compared to the NLO for 0.01 x 0.1, whereas the double Mellin resummation favors a larger gluon at higher x.
To quantify the behavior of the valence PDF at large x, we compute the effective β v parameter reflecting the exponent of the (1−x) βv term in Eq. (4), defined by [89][90][91] The results for the various resummation scenarios are shown in Fig. 3 as a function of x at the input scale, µ 0 , and compared with the ASV [45] analysis that fit the valence PDF to DY data using the cosine method of threshold resummation. In contrast to this work, ASV set x v between 0.55 and 0.7 and fixed the sea quark and gluon distributions from the earlier GRS analysis [48]. Consistent with previous studies [49,52,55,58], our NLO analysis shows a linear falloff of the valence PDF with β eff v ≈ 1 for x → 1. Inclusion of threshold resummation results in a wide variety of β eff v values, with the cosine and expansion methods yielding β eff v > 2, consistent with ASV [45], and as large as ≈ 2.6.
On the other hand, with the DM method the effective exponent is much closer to the NLO case, with β eff v ≈ 1.2 as x → 1. In fact, a recent study [41] pointed out problems with the MF methods neglecting some leading power effects, which are, in contrast, able to be accommodated with the DM method. This would suggest that β eff v values ∼ 1 are preferred in the more consistent DM approach, although it would be desirable to have empirical confirmation of this with additional observables sensitive to pion PDFs at large x.
Momentum fractions and pion mass decomposition.-A consequence of applying the NLL corrections to the DY cross section is that the large-x momentum of the valence quarks is redistributed to gluons at small x. The values of the total momentum fractions x i ≡ 1 0 dxxf i (x) for the different flavors are shown in Table I. Interestingly, while the shapes of the PDFs for the various resummed fits differ, the momentum fractions are rather stable, with ≈ 5%-6% of the momentum moving from the valence quark to the gluon sectors.
This has important implications for the decomposition  (7) 0.40 (5) of the pion mass into the quark and gluon energy and momentum, and trace anomaly contributions [92]. In particular, the gluon contribution to the mass is given by 3/4 of its momentum fraction, which amounts to 40(6) MeV, or ≈ 30% of the pion mass. This represents an increase of ≈ 14% on the gluonic fraction of the mass from the NLO analysis without resummation.
Outlook.-In the future, theoretical improvements will extend the treatment of resummation to NNLO corrections, allowing the analysis to be generalized to the qg, gg and qq channels and resummation effects on sea quark and gluon PDFs [38,41]. Concurrently, planned high luminosity tagged deep-inelastic scattering experiments at Jefferson Lab [93] and the future Electron-Ion Collider [94] on leading proton and neutron production at kinematics complementary to HERA will help isolate pion exchange contributions to valence and sea quarks, and test the Sullivan mechanism [87] in the low-mass region of the pion structure function.
The proposed COMPASS++/AMBER [95] experiment at CERN to measure pion-nucleus DY cross sections, including its p T dependence, would allow the large-x region to be further probed, providing a vital check on the Fermilab DY data and sensitivity to different nuclear targets. The availability of kaon beams would also allow global QCD analysis to elucidate for the first time the kaon's quark and gluon structure.