Towards the three-dimensional parton structure of the pion: Integrating transverse momentum data into global QCD analysis

We perform a new Monte Carlo QCD analysis of pion parton distribution functions, including, for the first time, transverse momentum dependent pion-nucleus Drell-Yan cross sections together with $p_{\rm T}$-integrated Drell-Yan and leading neutron electroproduction data from HERA. We assess the sensitivity of the Monte Carlo fits to kinematic cuts, factorization scale, and parametrization choice, and we discuss the impact of the various datasets on the pion's quark and gluon distributions. This study provides the necessary step towards the simultaneous analysis of collinear and transverse momentum dependent pion distributions and the determination of the pion's three-dimensional structure.


I. INTRODUCTION
The pion is one of the most enigmatic particles in nature, with unresolved questions about its most fundamental properties and behavior. Unlike all other hadrons, pions play a crucial role in QCD as the nearly massless Goldstone bosons associated with the breaking of chiral symmetry. At the same time, as strongly interacting bound states of quarks, antiquarks, and gluons (or partons), they reveal a familiar spectrum of momentum distributions at high energies, as do other hadrons and nuclei. Understanding how these contrasting manifestations of the same qq bound state arise dynamically at different energy scales from first principles remains a major challenge in QCD, and progress towards the ultimate resolution of the partonic-hadronic duality of the pion at present requires phenomenological input [1,2].

Recently the Jefferson Lab Angular Momentum (JAM) collaboration performed the first
Monte Carlo global QCD analysis of pion parton distribution functions (PDFs) [3] constrained by the traditional Drell-Yan dilepton production data in pion-nucleus scattering and data on leading neutron production in inclusive lepton-proton collisions. While the Drell-Yan data from CERN [4] and Fermilab [5] provide constraints on the valence quark PDFs at intermediate and large values of the parton momentum fraction x in the pion, the leading neutron data [6,7] were important in establishing the role played by gluons in the pion at low x. Global QCD studies of Drell-Yan and prompt photon data were previously performed in Refs. [8][9][10][11][12][13][14] in the context of single-fit analyses. In addition to the phenomenological studies, there has also been considerable interest recently in computing the pion's valence quark distribution in lattice QCD [15][16][17][18][19].
From another direction, recent progress in the theory of transverse momentum dependent (TMD) distributions has inspired exploratory studies of the transverse momentum structure of partons in the pion [20], which have attempted to describe the p T dependence of the differential pion-induced Drell-Yan cross section. The theoretical formalism developed in the 1980s by Collins, Soper and Sterman (CSS) [21] allows one to describe p T -dependent observables within the so-called "W+Y" framework (where p T is the transverse momentum of the virtual photon relative to the beam axis of the colliding hadrons). In this construction, the "W" term is computed in terms of TMD parton distributions in the region where p T Q, where Q is the invariant mass of the lepton pair. At such kinematics, the transverse momentum of the virtual photon is sensitive to the intrinsic transverse momentum of partons in the hadron. Crucially, the TMDs here are sensitive to the vacuum expectation value of Wilson lines, which is expected to be independent of the type of hadron involved. Confirmation of such universality for all hadrons would be one of the most important validations of QCD factorization theorems and their predictive power.
On the other hand, in regions of large transverse momenta p T ∼ Q, the spectrum is dominated by hard QCD radiation in which the photon recoils from hard radiated partons.
In such situations, the cross sections can be described by collinear factorization involving collinear PDFs. In practice, the CSS framework has been very successful in applications in collider environments such as W and Z production in large-p T pp and pp collisions at the Tevatron and LHC, respectively. One of the main reasons for this success is that collider setups offer sufficient phase space for the final states to allow perturbative QCD calculations to provide good descriptions of the data, in contrast to low-energy reactions, such as those in fixed-target experiments, where the applicability of the CSS framework is more limited.
In Ref. [22], it was found, for example, that even with the inclusion of O(α 2 s ) corrections, the theoretical predictions in the large-p T region significantly underestimate the measured cross section for low-energy experiments. Similar trends have been observed in semi-inclusive deep-inelastic scattering [23][24][25].
The situation is somewhat different in the small-p T region, where several analyses have attempted to describe the data in terms of TMDs [26][27][28]. These efforts have been partially successful, in that the analyses have been able to describe hadron multiplicities at small transverse momenta, making cuts on data with p T above some maximum value. A globally successful analysis, on the other hand, would seek to include the entire p T spectrum, describing the data in terms of the full W + Y treatment.
Such an analysis would be important for several reasons, including the fact that, on formal grounds, the applicability of the W-term is limited for regions where p T /Q 1. Despite that, data beyond the region of applicability, up to p T /Q ∼ 1, are still often included in lowenergy TMD analyses, thereby potentially biasing the extraction of the TMDs. Additionally, the large transverse momenta give access to PDFs at large parton momentum fractions, making it more feasible, for instance, to probe gluon distributions at large x, through gluon initiated subchannels that arise at the same order in α s as quark channels.
Given that TMDs are now considered to be a key tool in describing the structure of hadrons, the question arises of whether it can be confirmed that collinear factorization can serve as an equally reliable tool for describing inclusive particle production at high p T .
Motivated by the state of the art of TMD phenomenology, as well as the broad interest in exploring the inner structure of pions and the recent success in extracting pion TMDs, in this paper we analyze the large transverse momentum region of Drell-Yan lepton-pair production in fixed target pion-nucleus scattering. We will demonstrate that, contrary to proton-induced Drell-Yan reactions, in the case of the pion, it is indeed possible to describe the large transverse momentum region and include the available p T -dependent data consistently within a global analysis of pion PDFs constrained by p T -integrated data, as in Ref. [3].
This analysis is a first step towards a full combined analysis of all available pion-induced data at low and high p T to simultaneously determine collinear and TMD pion PDFs within the same theoretical framework.
In Sec. II of this paper, we outline the theoretical formalism used in our analysis to describe p T -dependent and p T -integrated Drell-Yan lepton-pair production cross sections in pion-nucleus scattering, along with the inclusive electroproduction of forward neutrons in electron-proton collisions. We perform the analysis within the JAM Monte Carlo fitting framework, used in previous successful analyses of pion PDFs [3], unpolarized nucleon PDFs and fragmentation functions [29,30], and spin-dependent nucleon PDFs [31], as well as the nucleon's transversity distributions [32] and other TMD functions [33]. The datasets used in this analysis are described in Sec. III, which include Drell-Yan lepton pair production cross sections in pion-nucleus scattering from CERN [4] and Fermilab [5], including both p T -differential and p T -integrated data, as well as leading neutron electroproduction cross sections from HERA [6,7]. In Sec. IV we describe the Monte Carlo methodology used to extract the PDFs through Bayesian inference. The results for the pion's valence and sea quark and gluon PDFs are presented in Sec. V, where we examine the dependence of the analysis on kinematic cuts, factorization scale, and choice of PDF parametrization. In Sec. VI we compare the results of the analysis with PDFs in the proton, and finally we summarize our findings and discuss future extensions of this work in Sec. VII.

II. THEORETICAL FRAMEWORK
In this section we present the theoretical framework used in the present analysis of PDFs in the pion, including the relevant formulas for p T -differential and p T -integrated cross sec-tions for the Drell-Yan process in pion-nucleus scattering [34], along with the semi-inclusive forward neutron production cross sections in electron-proton collisions. We provide only the relevant formulas needed for our immediate calculations; more details about the formalism for both types of reactions can be found in Refs. [3,[35][36][37].
A. Drell-Yan lepton-pair production Within the CSS framework the p T -differential cross section for the inclusive production of a lepton pair + − in the high-energy collision of hadrons A and B, A B → + − X, can be expressed schematically in the form where "W", "FO" and "ASY" refer to the W-term, the fixed-order term, and the asymptotic contribution, respectively [38], and the kinematic variables are the lepton pair's invariant mass squared Q 2 , rapidity y, and transverse momentum p T . The corrections to the formula (1) appear in the form of ratios of masses m to the large scale Q.
The W-term in Eq. (1) is computable in terms of TMD PDFs, and is a valid approximation to the cross section in the low-p T region, p T /Q 1, with errors O(p T /Q). In contrast, the FO term is applicable in the region where p T ∼ Q, with corresponding errors O(m/p T ). Finally, the asymptotic term ASY in Eq. (1) is, by construction, the large transverse momentum approximation to the W-term, as well as the small transverse momentum approximation to the FO-term. Defining the "Y-term" by Y ≡ FO − ASY, in the small-p T limit one has lim p T /Q→0 Y = 0 and the cross section is dominated by the W-term, while at large p T by construction lim p T ∼Q (W − ASY) = 0 and the cross section is dominated by the FO-term. In the present analysis, we will focus on the latter region and examine in detail the FO-term. Further discussion about the W and ASY contributions can be found at Ref. [21].
To produce a virtual photon with a large transverse momentum, the photon must recoil against a hard radiation, examples of which are shown in Fig. 1. Since the scale of such radiation is large, it cannot be associated with an intrinsic transverse momentum of quarks and gluons in the nucleon, but instead involves hard radiation, which is calculable perturbatively in QCD. It will be convenient to parametrize the virtual photon's four-momentum q in terms of its invariant mass Q, transverse momentum p T and rapidity y, where Q T ≡ Q 2 + p 2 T is the transverse mass of the photon. Schematically, the FO-term can then be written, to all orders in the strong coupling α s , as where f A a (f B b ) is the collinear PDF for flavor a (b) in hadron A (B), evaluated at parton momentum fraction x a (x b ) at a renormalization scale µ, and H FO a,b is the perturbatively calculable short-distance cross section. Note that the integration limits in Eq. (3) depend on all three variables y, Q 2 and p T .
The production of a large transverse momentum photon can proceed via two types of partonic hard scattering subprocesses, qq → + − g ( Fig. 1(a)) and qg → + − q ( Fig. 1(b)), which start at O(α s ) in the QCD coupling constant. At this order, the final state phase space constrains one of the parton momentum fraction integrals so that x a and x b are not independent, but are related by where and the minimum value of x a is given by As p T increases, the lower limit of x min a increases, so that the cross section becomes increasingly more sensitive to PDFs with large parton momentum fractions. While the perturbative QCD calculations for H FO a,b are known up to O(α 2 s ) [39], in the present analysis we will restrict ourselves to O(α s ) in order to test the applicability of the large transverse momentum region using the lowest order in pion induced Drell-Yan reactions.
It is instructive to contrast the theoretical calculation for the p T -differential cross section (1) with the corresponding cross section integrated over all p T . In this case the p Tintegrated cross section can be written in the more familiar form of collinear factorization as where in this case the limits of integration no longer depend on p T . The hard cross section H DY a,b starts at O(α 0 s ), and in our analysis we compute corrections up to O(α s ). Our study is the first attempt to include both p T -differential and p T -integrated pion-nucleus Drell-Yan data [4,5] on the same footing, taking advantage of the fact that the p T -dependent cross sections provide access to a larger region of parton momentum fractions relative to the p T -integrated case.

B. Leading neutron electroproduction
In addition to the inclusive lepton-pair production cross sections in pion-nucleus scattering in Eqs. (1) and (7), following Ref. [3], we supplement the analysis with data from leading neutron (LN) electroproduction in ep collisions at HERA [6,7] in the very forward region, ep → enX, in order to better constrain the pion PDFs at small values of parton momentum fraction x in the pion. By detecting a forward neutron and lepton in the final state, the hard scattering is expected to occur between the exchanged virtual photon (fourmomentum q = − , where and are the incident and scattering lepton momenta) and the pion in the limit where the momentum transfer squared t = (p − p ) 2 from the proton (p) to the final neutron (p ) is small [40,41].
The triply differential LN cross section is typically parametrized by the LN structure function, F LN 2 , which is a function of the Bjorken scaling variable x B = Q 2 /2p · q, the invariant mass squared of the virtual photon, Q 2 = −q 2 , and the longitudinal momentum fraction carried by the detected neutron, x L = p · q/p · q, where y e = q · p/ · p ≈ Q 2 /xs is the lepton inelasticity, and s = ( + p) 2 is the total center of mass energy ( √ s ≈ 300 GeV at HERA kinematics).
In the limit where |t| → 0 andx L ≡ 1 − x L → 0, one expects the charge exchange process γ * p → nX to be dominated by the emission and absorption of pions [40][41][42][43][44]. In this region, the chiral symmetry properties of QCD suggest that the peripheral structure of the nucleon can be described through interactions between the probe and the "cloud" of pseudoscalar mesons associated with the long-range structure of the nucleon. Formally, the effects of this structure on the total LN cross section can be computed within chiral effective field theory, by matching twist-two partonic operators with corresponding hadronic operators having the same quantum numbers [45][46][47].
This leads to a representation for the matrix elements of the twist-two operators between nucleon states in terms of matrix elements of the hadronic operators and matrix elements of the twist-two operators between the hadronic states, which can in turn be expressed through moments of the corresponding light-cone correlation functions. The validity of this construction for all moments implies its veracity also in momentum space, which allows the LN structure function to be written as where f πN is the light-cone momentum distribution of pions in the nucleon (or nucleon → nucleon + pion splitting function) [48][49][50][51][52], and F π 2 is the inclusive structure function of the pion evaluated at x π = x B /x L . Note the factor 2 in Eq. (9) is an isospin factor for the application of the p → nπ + fluctuation, whose distribution is related to the p → pπ 0 fluctuation by f π + n = 2f π 0 p = 2f πN . Within the chiral effective theory framework, at the one-pion loop level the splitting function f πN is given by [40,41,[48][49][50][51][52][53] where D πN is the pion propagator, and g A = 1.267 is the axial charge, f π = 93 MeV is the pion decay constant, and M and m π are the nucleon and pion masses, respectively. The function F regulates the ultraviolet divergences as k 2 ⊥ → ∞, and it is determined phenomenologically in terms of a cutoff mass parameter, Λ [3,36,54]. The pion structure function F π 2 can be written in the standard collinear factorized form, in terms of the PDFs of flavor i in the pion, f π i , and the inclusive DIS hard scattering function, H DIS i . Having summarized the theoretical formulas relevant for our global PDF study, in the next sections we focus on the details of the datasets that will be fitted, and the methodology employed for the Monte Carlo analysis.

III. DATASETS
Following the previous JAM pion PDF analysis [3], we fit pion-tungsten Drell-Yan leptonpair production cross sections from the NA10 experiment at CERN [4] and the E615 experiment at Fermilab [5], together with LN electroproduction data from the ZEUS [7] and H1 [6] experiments at the HERA ep collider. In addition, we include, for the first time, data on the p T -differential Drell-Yan cross sections from E615, which have not been utilized in global PDF analyses to date. The kinematic coverage for the Drell-Yan and LN datasets is shown in Fig. 2.
The Drell-Yan cross sections are typically presented in terms of √ τ = Q 2 /s, where Q 2 is the invariant mass squared of the virtual photon and s is the total incoming center of mass energy, and the Feynman scaling variable x π = 1 2 (x F + x 2 F + 4τ ), and the H1 [6] and ZEUS [7] leading neutron electroproduction data, for which x π = x B /x L . The NA10 datasets include two different pion beam energies, 194 and 286 GeV (see Table I).
values for the tungsten PDFs from the nuclear PDF analysis by Eskola et al. [55]. Using other nuclear PDF sets, such as from the nCTEQ group [56], leads to indistinguishable results on pion PDFs or goodness of fit.
To ensure that the Drell-Yan lepton-pair production process is dominated by leadingtwist partonic subprocesses, we avoid contributions from the J/ψ and Υ resonances by only selecting events with m J/ψ Q m Υ . Furthermore, since with increasing Q the available phase space for the hard factorization becomes smaller, to avoid edges of phase space we restrict the kinematical reach at a fixed s to Q values that are not too large. In practice, we find these conditions to be optimally satisfied for the Q range given by 4.16 < Q < 7.68 GeV.
The upper limit is slightly more conservative than that which was adopted in Ref. [3] in order to avoid increasing deviations between data and theory in the highest-Q bin.
For the range of x F considered, to avoid contamination from exclusive channels and remain within the region where factorization is valid, we implement the same restrictions as in Ref. [3], namely, x F < 0.9. To test if a smaller x F would be more appropriate, we also performed fits to the Drell-Yan data with x max F cuts of 0.6, 0.7, and 0.8, but we found that the χ 2 and resulting PDFs were essentially unchanged from the results with x max F = 0.9. To retain the largest number of data points that can be simultaneously described in the global analysis, we present results for data in the range 0 < x F < 0.9. This leaves 61 data points from the E615 experiment and 56 from NA10, for a total of 117 data points representing the p T -integrated Drell-Yan cross section.
The p T -dependent lepton-pair production cross sections are from the same E615 Drell-Yan pion-tungsten experiment [5], which measured the p T spectra both integrated over x F , d 2 σ DY /dQdp T , and integrated over Q 2 , d 2 σ DY /dx F dp T . The p T -dependent E615 data extend up to p max T = 4.875 GeV, which is ≈ 25% of the maximum p T allowed for the Q 2 ranges used, and the data are still far from the kinematic edge. The lower limit on p T is taken to ensure the applicability of the collinear approximation for the FO term in Eq. (3), and the value p min T = 2.7 GeV was determined phenomenologically, as discussed in Sec. V A below.
For the x F -integrated cross section, the fitted data span the range 4.3 < Q < 7.2 GeV, with a total of 34 data points, while for the Q-integrated cross section (which is integrated over 4.05 < Q < 8.55 GeV) data for x F 0.9 are included, with a total of 49 data points.
For the LN electroproduction data, the H1 experiment at HERA [6] measured the F LN 2 structure function, related to the triply differential LN cross section in Eq. (8), while the ZEUS Collaboration [7] measured the ratio of LN to inclusive cross sections, for x L bin sizes ∆x L . For the denominator of Eq. (13) we use the JAM19 PDFs [29] to compute the inclusive proton cross sections. As seen in Fig. 2, the LN data extend to much lower values of x π than the Drell-Yan data, x π 10 −3 , allowing more direct determination of the pion's sea quark and gluon distributions. Indirect constraints on the gluon PDFs are obtained from Q 2 evolution over the large range of Q 2 values covered by the HERA data, from Q 2 ∼ a few GeV 2 to 10 3 GeV 2 . For the x L dependence, we restrict the analysis to the region of large x L , where most of the incident proton's longitudinal momentum is carried away by the detected neutron, and one can approximate the LN cross section by the exchange of a soft pion [43,44]. We follow the systematic studies of χ 2 performed as a function of the x L in Refs. [3,36] and apply the cut x L > 0.8 found there.

IV. MONTE CARLO METHODOLOGY
In this section we describe the analysis procedure used to extract the pion PDFs through Bayesian inference. The basic strategy follows previous JAM analyses of pion [3] and proton [29,31] PDFs, with some improvements specific to the present application.

A. Parameter inference
As in earlier global PDF analyses, we use the following standard template for the functional form of the nonperturbative distributions, which are parametrized at the scale µ 2 0 , where B is the Euler beta function, and the normalization N i is defined to correspond to the average momentum fraction of the pion carried by parton i. The input scale is defined For the pion PDFs we parametrize 3 degrees of freedom: the valence quark distribution, q π v , the sea quark distribution, q π s , and the gluon distribution, g π . All of the light-quark PDFs can be expressed in terms of the distributions q π v and q π s , using isospin and SU(3) flavor symmetry. In particular, for the valence quark distribution while for the sea quarks we take q π s ≡ū π + = d π + = s π + =s π + . The heavy quark flavors in our analysis are generated perturbatively via QCD evolution by solving the DGLAP equations using the zero-mass variable flavor number scheme evolved up to next-to-leading-logarithmic accuracy.
In our numerical analysis we have found that fitting the available data with the 5- The free parameters are inferred by sampling the Bayesian posterior distribution where a = {N i , a i , b i , . . .} represents the vector of parameters, and the likelihood function L is given by The argument of the exponential in (16) is given by the χ 2 function, where e labels different experimental datasets with data points d e i and corresponding theory values t e i (a), α e i are point-to-point uncorrelated uncertainties added in quadrature, and β e k,i are the corresponding point-to-point correlated uncertainties. We allow the theory to be shifted additively by an amount k r e k β e k,i , with nuisance parameters that are optimized via dχ 2 /dr e k = 0. In addition to fitting the shape parameters in Eq. (14), we also fit the overall multiplicative normalization factors n e in Eq. (17) for all 7 datasets used in this analysis, with a Gaussian penalty controlled by the quoted experimental normalization uncertainties δn e . The total number of free parameters is therefore 15, comprising the 7 PDF shape parameters and 1 pion splitting function parameter for the LN data, and the 7 data normalizations (3 for the p T -integrated E615 and NA10 194 GeV and 286 GeV data, 1 each for the x F and Q-integrated E615 data, and 1 each for the H1 and ZEUS LN datasets).
where {a k } are the ensemble of Monte Carlo parameters of size N drawn from the posterior distribution (15). This definition avoids assuming Gaussianity for the parameter distributions, as in the traditional Hessian method [57][58][59].
The Monte Carlo samples are constructed from Eq. (15) using data resampling and performing multiple χ 2 minimizations. The starting values of the parameters for each minimization are selected randomly from flat priors within the nonvanishing region of the prior, and each data point d e i is shifted Gaussianly within the quoted uncorrelated uncertainties α e i added in quadrature. In practice, the ranges of the a i parameters that control the x → 0 behavior are restricted to guarantee integrability of the pion PDFs for the valence and momentum sum rules.

B. Multiple solutions
As is the nature of Monte Carlo analyses, typically multiple solutions are obtained for many of the PDF parameters. In order to discriminate between the various solutions, we employ a k-means clustering algorithm [60,61], similar to that used in the recent JAM19 analysis of proton PDFs [29]. Discriminating between the various clusters of solutions is then achieved on the basis of the total χ 2 for each of the clusters, or eliminating clusters that build up at edges of certain parameters with physically motivated boundaries.
Initially, fits were performed with the ranges of the b i parameters for the gluon and sea quark PDFs restricted so that these distributions do not exceed the valence quark PDF at large x. After identifying the cluster with the smallest χ 2 , optimal priors were generated from the replicas in that cluster. Each of the parameters in the priors was generated in a random normal distribution with the mean and standard deviation of the best cluster's parameters. Solutions were removed in which the b s parameter for the sea quark PDF was smaller than the valence quark b v parameter, to ensure the natural dominance of the valence distribution at large x. Cuts were also made to avoid edge effects in parameter space from the artificial buildup of replicas at the boundaries of the parameter ranges associated with the gradient descent algorithm.
The spurious buildup at the boundaries is illustrated in Fig. 3(a) and (b), where the distribution of the gluon a g and b g parameters is shown versus the gluon normalization parameter N g in scatterplot form. The k-means algorithm was then used to separate clusters of solutions based on these parameters. A buildup of solutions is seen at the lower boundaries of both the a g and b g parameters, with the small-x parameter a g > −2 constrained to ensure a finite integrated momentum fraction carried by gluons, and b g > 0 to avoid a diverging gluon PDF as x → 1. The solutions corresponding to parameters displaying obvious edge effects should not be regarded as true solutions, and merely reflect of the imperfections of the Monte Carlo algorithm adopted.  x F -dependent (blue circles) p T -differential Drell-Yan E615 data [5] used in this analysis.

V. QCD BAYESIAN ANALYSIS
Having outlined the theoretical framework for our analysis, along with the datasets used and methodology employed to fit them, in this section we present the results of the fits, including data-to-theory comparisons and the resulting pion PDFs together with their uncertainties. We begin the presentation of the results with a discussion of the dataset selection and justification of the cuts employed to ensure the integrity of our theoretical treatment of the data.

A. Data selection
As discussed in Sec. III, for the p T -integrated Drell-Yan lepton-pair production data, we implement mostly the same kinematic cuts as in the previous JAM analysis [3], namely, 0 < x F < 0.9 and a slightly modified range 4.16 < Q < 7.68 GeV, in order to avoid edges of phase space and ensure dominance of leading power factorization. Similarly, for the HERA LN data, to restrict the analysis to the region where the pion-exchange process dominates the neutron production, we apply the cut x L > 0.8, as in the earlier studies [3,36].

B. Data versus theory comparison
With these cuts on the data, the results of the global fit are summarized in Table I, where for each experiment we list the number of data points fitted, the data normalization parameter, and the corresponding χ 2 dat value. The overall χ 2 is 262 for a total of 308 data points, with the total χ 2 dat = 0.85 per datum. Moreover, we find that fairly good fits are obtained to all of the datasets, with χ 2 dat ≈ 1 for each experiment. An exception is the ZEUS LN dataset, with a χ 2 dat ≈ 1.5, which partly reflects the relatively small uncertainties on these data compared with the H1 results.
The comparison between theory and experiment is illustrated in Figs. 5 -7, which show data over theory ratios for the p T -integrated Drell-Yan data from E615 [5] and NA10 [4], the leading neutron electroproduction data from H1 [6] and ZEUS [7], and the p T -differential   for two bins of momentum fractionx L carried by the exchanged charged particle (pion), restricted tox L < 0.1 and 0.1 <x L < 0.2 to ensure pion exchange dominance [3,36].
Within the quoted uncertainties, the H1 data can be well described by our fit across all the kinematics shown, with a χ 2 dat value of ≈ 0.4 per datum. The systematic uncertainties on the H1 data are largest for the lowest-x L bin, and since the magnitude of the F LN 2 structure function increases withx L , the relative uncertainties for the data are thus largest forx L < 0. For the ZEUS data, since these are presented as a ratio of semi-inclusive to inclusive structure functions, the uncertainties are smaller than for the absolute structure functions from the H1 experiment, with a number of the correlated systematic uncertainties canceling.  Data-to-theory ratios for the p T dependence of the Drell-Yan cross sections d 2 σ DY /dp T dQ (top) and d 2 σ DY /dp T dx F (bottom) from the E615 experiment [5]. The yellow bands represent the theoretical uncertainties, and the scale is set to µ = p T /2 for the PDFs.
Consequently, these data provide the strongest constraints, and give the highest overall χ 2 of all the datasets fitted in this analysis, with χ 2 dat ≈ 1.5. There is some tension with the data at the smallest x π values, for the lowest Q 2 bin, but generally the agreement between theory and data is quite good.
The p T -differential Drell-Yan cross sections from E615 [5], which includes the new data considered in our analysis, are displayed in Fig. 7 as data-to-theory ratios versus the p T of the lepton pair. Choosing the QCD factorization to be µ = p T /2 (see Sec. V E), fairly good fits are obtained for both the x F -integrated data (χ 2 dat = 1.08) and the Q-integrated data (χ 2 dat = 0.85) shown in Fig. 7, although a somewhat small value of the overall data normalization factor is necessary for the Q-integrated (n e = 0.50) compared with the x F -integrated (n e = 0.83) cross sections. The difference between the normalization factors for the various datasets reflects possible tensions among the data, which can affect deviations of n e outside of the normalization uncertainties reported by the experiments. The uncertainties on the p T -dependent data are generally much larger than on the p T -integrated cross sections in Fig. 5 (note the vertical scale on the data/theory ratios in Fig. 7), and grow with p T , increasing markedly as x F → 1.
Future, higher-precision data on the p T dependence of Drell-Yan lepton-pair production cross sections would provide stronger constraints on the fits than the currently available data.
A program of pion-induced and kaon-induced Drell-Yan experiments is being discussed in connection with the proposed COMPASS++/AMBER facility at CERN [62]. Nevertheless, given the present difficulty in reconciling the p T dependence of proton-induced Drell-Yan data in the large-p T region with the latest theoretical tools [22], it is noteworthy that a consistent description of p T -differential and p T -integrated pion-nucleus cross sections can be achieved within the same collinear factorization framework.

C. Parton distributions
The resulting spectrum of PDFs from our Monte Carlo analysis of Drell-Yan (p Tintegrated and p T -differential) and LN data, is illustrated in Fig. 8, where we present the individual valence quark q π v , sea quark q π s , and gluon g π distributions (where g π is scaled by a factor 1/10) at a scale µ 2 = 10 GeV 2 . For clarity, we show a random sample of 150 replicas out of the total ≈ 800 replicas from our analysis. Since some of the samples are faulty because of the imperfections of the Monte Carlo algorithm rather than viable physical solutions, the solutions displaying edge effects in Fig. 3 have been removed from the final sample. The larger spread of solutions for sea quark and gluon PDFs compared with the valence quark distribution reflects the weaker constraints on the pion sea at small x values from existing (mostly LN) data compared with those from the Drell-Yan data. More data at low x values would clearly be useful for constraining the sea quark and gluon PDFs.
A comparison of the JAM pion PDFs and their 1σ uncertainty bands with results from previous global analyses is shown in Fig. 9 at a scale µ 2 = 10 GeV 2 . Generally good agreement, especially for the valence quark distribution, is found with the recent xFitter parametrization [14], which fitted the p T -integrated Drell-Yan data along with data on prompt photon production from the WA70 experiment [63]. At low x values the xFitter sea quark PDF has somewhat larger uncertainty, which reflects the fact no LN data were used in this fit, hence the diminished constraining power for the PDFs in this region. For reference, the older GRV parametrization [11], which preceded the HERA data and does not have PDF uncertainties, is also shown. Compared to the JAM result, the GRV fit has a slightly softer valence PDF, ∼ 25% smaller at intermediate  9. Comparison of the pion valence quark q π v , sea quark q π s , and gluon g π (scaled by 1/10) PDFs from the current JAM analysis (red bands) at µ 2 = 10 GeV 2 with the xFitter results [14] (yellow bands) and the GRV parametrization [11] (blue lines). The uncertainty bands represent 1σ CL. large uncertainties. Clearly the biggest overall impact on the PDFs uncertainties is scenario (ii), in which the addition of the HERA LN data constrains significantly the small-x region for the gluon and the sea distributions, with modest effect on the valence distribution. This is consistent with what was previously observed in Ref. [3].
The novel addition of the p T -dependent Drell-Yan data in scenario (iii), has a modest impact on the shapes of the pion PDFs and their uncertainties. The strongest impact is on the gluon distribution at large values of x, x 0.3. This may be expected, given the sensitivity of the p T -differential cross section on the pion's gluon PDF at lowest order in α s . However, since the cross section at large x is still mostly dominated by contributions from valence quarks, the overall impact on the glue is not overwhelming. In other kinematic regions, the reduction in the PDF uncertainties after inclusion of the p T -dependent Drell-Yan data is also relatively small, which reflects the larger errors of these data in Fig. 7 than for the p T -integrated Drell-Yan and LN data in Figs. 5 and 6, respectively.
Interestingly, the behavior of the valence PDF at large x is consistent with a ∼ (1 − x) shape, as was found in Ref. [3]. The p T -dependent Drell-Yan does not alter this conclusion.
Further discussion on the large-x behavior of the pion PDF in the presence of threshold resummation will be discussed elsewhere [37].
The impacts of the different datasets can be further explored by considering the momentum fractions carried by the individual flavors i (= v, s, g) of the pion, defined as at a scale µ 2 . The Monte Carlo samples for the pion momentum fractions are shown in Fig. 11 at the input scale µ 2 = m 2 c as histograms, for each flavor, for the scenarios (ii) and (iii) described above, with the "normalized yield" defined by the area under the histogram being unity. The momentum fractions with PDFs extracted from only the Drell-Yan data (scenario (i)) are not shown in Fig. 11 to avoid overcluttering, but their momentum fractions are similar to what was found previously in Ref. [3]. Namely, in this scenario the total quark valence momentum fraction is relatively well constrained to be 0.59 (1), while the sea quark and gluon fractions, 0.28(10) and 0.13 (11), respectively, have significantly larger uncertainties.
For the fit with the p T -integrated Drell-Yan and LN data, labeled as "no DYp T " in approximately switch in magnitude with the addition of the LN data, and the uncertainties on both decrease by a factor from ∼ 2 to ∼ 4. With the better determined sea quark and gluon fractions, the valence quark fraction decreases by ≈ 10% in order to satisfy the momentum sum rule.
The addition of the p T -dependent Drell-Yan data, for the fit labeled as "full" in Fig. 11, does not affect the momentum fraction appreciably, resulting in a slight reshuffling of strength between the sea quark and gluon sectors, but within the 1σ uncertainties, with    Table II summarizes the momentum fractions for all three scenarios at both the input scale and at the scale µ 2 = 10 GeV 2 . Notably, when the PDFs are evolved to the higher scale, the momentum fraction carried by valence quarks decreases by ∼ 0.1, whereas the gluon momentum fraction increases considerably, and a smaller increase occurs for sea quarks. The uncertainties on the moments decrease as the moments are evolved to the higher scale.

D. Flavor decomposition of observables
A deeper understanding of the impact on the pion PDFs from the various observables considered in our analysis can be obtained by examining the relative flavor contributions to the observables, especially from the less well constrained sea quark and gluon distributions.
In particular, we wish to study how the different flavors build up the LN structure function, which is sensitive to PDFs at small x, and the p T -differential Drell-Yan cross section, in which the gluon PDF enters at the lowest order through the qg channel at O(α s ). In Fig. 12 we show the valence quark, sea quark and gluon contributions to the p T -integrated and p Tdifferential Drell-Yan cross sections and to the LN structure function, relative to the total cross sections, at selected kinematics.
Firstly, for the p T -integrated Drell-Yan cross section, d 2 σ DY /d √ τ dx F , shown in Fig. 12 versus x F at fixed √ τ = 0.288, the most striking observation is the near dominance of the valence quark contribution to the total cross section across the entire x F range. The sea quark contribution grows at small x F , but does not exceed ≈ 15% of the total, while the gluon component is almost negligible over all x F . Analyses that include only the p Tintegrated Drell-Yan data yield valence quark PDFs that are reasonably well determined, but leave the sea quark and gluon distributions almost totally unconstrained.
As observed in Fig. 10, the greatest impact on the sea quark and gluon distributions is from the HERA F LN 2 LN structure function data, shown in Fig. 12 versus x π at a representative value of Q 2 = 24 GeV 2 . Interestingly, the dominant contribution to F LN 2 at low x π (x π 0.01) is from sea quarks, while the valence quark contributions become larger only at x π 0.1. The gluon contribution is strongly suppressed across all x π , largely because it  F LN 2 versus x π at Q 2 = 24 GeV 2 , (right) triply-differential p T -dependent Drell-Yan cross section d 3 σ DY /dp T dx F dQ versus Q 2 at p T = 2.7 GeV and x F = 0.05. The gray horizontal band represents the range of Q 2 for current and upcoming pion-induced DY measurements.
appears only at NLO, and at intermediate x π is actually negative. However, since with the Drell-Yan and LN data both the valence quark and sea quark PDFs become more strongly constrained, as Fig. 10 indicates, the uncertainties on the gluon PDF also decrease because of the momentum sum rule. Thus the LN and Drell-Yan data play rather complementary roles in constraining valence and sea sectors of the pion.
More direct constraint on the gluon PDF is expected from the p T -differential Drell-Yan cross section, which, because of the large p T of the exchanged virtual photon involved in the process, requires a hard gluon radiation (see Fig. 1). In Fig. 12, we show the ratio of flavor contributions to the triply differential Drell-Yan cross section, d 3 σ DY /dp T dx F dQ, versus Q 2 at fixed x F = 0.05 and at the minimal transverse momentum of the E615 data used in this analysis, namely p T = 2.7 GeV. Note that triply differential cross section data for the pion-induced Drell-Yan process currently do not exist, but are shown here simply to illustrate the kinematics. based on past experiments, such as E615 at Fermilab [5], and projections for the future COMPASS++/AMBER experiment at CERN [62]. Accessing and interpreting lower-Q 2 data in terms of collinear factorization is problematic, however, since here meson resonances appear, which do not lend themselves to simple partonic interpretations. Nevertheless, the p T -dependent Drell-Yan data are the most challenging to describe, as Fig. 7 illustrates, and more precise data on the p T dependence, especially for the triply differential cross section, data would be valuable [62].

E. Scale dependence
Before concluding the discussion of the phenomenological results of our analysis, we review the problem of factorization scale setting for the p T -dependent Drell-Yan data. Inclusive processes, such as deep-inelastic scattering and (p T -integrated) Drell-Yan lepton-pair production, are based on a single hard scale that characterizes the reaction, usually taken to be the invariant mass of the exchanged virtual photon. In those cases, the scale dependence, introduced in QCD factorization to evaluate the short-distance cross section and the evolution of the PDFs, is chosen to be the hard scale in order to optimize the perturbative convergence. In contrast, a second scale enters when considering cross sections differential in p T , and the scale choice for the factorization at large transverse momentum can be chosen to be either set by p T or by Q.
Formally, the choice of scale should not affect the results appreciably, since the p Tdependent cross sections at large p T are factorized under the condition that p T ∼ Q. In practice, however, using a lowest order perturbative calculation for the large transverse momentum region does not guarantee that the resulting cross section will be insensitive to the choice of scale. It is important, therefore, when performing an analysis with lowest order accuracy, to vary the scale around the two hard scales in the problem in order to estimate the sensitivity of the extracted PDFs to the scale choice.
We performed fits to the p T -dependent Drell-Yan data by selecting several different scale choices, ranging from µ = p T /2 to µ = 2p T and including µ = Q, as illustrated in Fig. 13.
The scale dependence was found to correlate strongly with the normalization of the datasets, but the overall effect on the extracted PDFs was mild. The best agreement with the p Tdependent Drell-Yan data was found for the smallest value of the scale, µ = p T /2, for values of the scale, with the fit for µ = p T yielding χ 2 dat = 2.11, while for µ = 2p T giving χ 2 dat = 3.63 for the p T -dependent data. When using the scale µ = Q, the quality of the fit is also worse, with χ 2 dat = 3.18. The PDFs with the scale µ = 2p T are almost identical to those for µ = p T , and are not shown in Fig. 13. The fitted normalization factors, n e , for all the scale settings other than µ = p T /2 were found to be away from unity and closer to 0.5, which was the minimum value that we have allowed in our analysis. The correlation between the scale dependence and the fitted normalization factor indicates potentially the need to include higher-order corrections that can restore the normalization factor to be closer to unity, but we leave such studies for future work. the gluon, are different. Attempts at comparing the inclusive proton structure function and the pion structure function extracted from neutron electroproduction were made in the HERA analyses of LN cross sections [6,7], using data obtained under the same experimental conditions. Here we revisit the proton versus pion structure function comparison in the context of a global QCD analysis, contrasting the results from the current pion analysis with the recent JAM19 proton PDFs obtained using similar Monte Carlo methodology [29].

A. Parton distributions
While the valence quark PDFs in the proton and pion are normalized differently, the momentum sum rule should hold for both hadrons, so that the momentum fractions carried by valence quarks, sea quarks and gluons can be compared directly. In Fig. 14, the momentum fractions carried by each of these species in the proton and in the pion are shown in the form of histograms for the normalized yield, evaluated at the input scale, µ 2 = m 2 c . The pion fractions correspond to the "full" results in Fig. 11, while the proton fractions are computed from the JAM19 proton PDFs [29], with mean values and uncertainties given by x p s = 0.147(4), for the valence, sea and gluon, respectively at the input scale. The most striking difference between the pion and proton results is the significantly narrower distributions for the latter, indicating much smaller uncertainties, about an order of magnitude, on the proton PDFs.
Interestingly, the valence quark momentum fraction in the proton is ∼ 10% smaller than that in the pion (which is ≈ 3σ, with the small proton uncertainties), while the central value of the gluon momentum fraction in the proton is ∼ 25% larger, although compatible at the ≈ 1.5σ level with the much larger gluon pion uncertainties. The total momentum carried by sea quarks is ≈ 15% for both hadrons, with ≈ 7 times smaller uncertainty for the proton.
For the case when the scale is evolved to µ 2 = 10 GeV 2 , the proton fractions are given by and as was the case in the pion, the valence quark momentum fraction decreases, while the gluon fraction increases more so than the sea quark with evolution.
The larger uncertainties on the x dependence of the pion PDFs compared to those for the proton are illustrated in Fig. 15, where we show the pion and proton total sea quark and gluon densities 1σ uncertainty bands relative to the mean values of the proton PDFs, at a scale µ 2 = 10 GeV 2 . Although the momentum fractions carried by sea quarks in the pion and proton are similar, as witnessed in Fig. 14, the shapes of the PDFs themselves are rather different. In particular, the pion sea is significantly harder than the proton sea, with the latter ∼ 30% − 40% larger for x 0.01. The pion and proton sea quark PDFs are comparable at x ∼ 0.1, above which the pion's sea quark distributions become larger. The larger sea quark PDFs at low x in the proton has implications for the structure function comparison, as discussed below.
For the gluon distribution, the comparison in Fig. 14  and ZEUS experiments reported that the ratio r was approximately a function of x L only, which suggested that F LN 2 (x B , Q 2 , x L ) ∝ F p 2 (x B , Q 2 ) at a given x L , and the proportionality was shown to hold for intermediate to high values of x L .
In the factorized approximation (9) in terms of pion exchange, the F LN 2 structure function is expressed in terms of the pion-nucleon splitting function f πN (x L ) and the pion structure function F π 2 (x π , Q 2 ), with x π = x B /x L . For the ratio r to be a function of x L only, the x B (and Q 2 ) dependence in the numerator and denominator need to cancel, and since the splitting function depends only on x L , the F π 2 and F p 2 structure functions would be proportional at a given x L , with a proportionality factor α, F π 2 (x π , Q 2 ) ≈ αF p 2 (x B , Q 2 ). Using the results from our analysis, we are able to confront directly whether the equality holds for a scaled proton structure function by computing F p 2 (x B , Q 2 ) and F π 2 (x π , Q 2 ) explicitly. An alternative model that gives reasonable agreement with the LN data is the color-dipole 0.001 0.01 0.1 1 16. Comparison of the pion structure function from the full analysis, F π 2 (red lines and 1σ uncertainty bands) and with Drell-Yan data only, F π,DY 2 (gray lines and 1σ yellow uncertainty bands) with models of F π 2 based on the proton structure function, evaluated at Q 2 = 10 GeV 2 (top) and at the input scale Q 2 = µ 2 0 = m 2 c (bottom). The models include F p 2 scaled by a factor α = 0.65 (blue dashed lines) and scaled as 2 Here, the 2/3 normalization factor comes from the ratio of the pion to proton eigenfunctions of the color-dipole BFKL equation, whereas the x π dependence comes from an inverse proportionality of the effective dipole-dipole scattering energy, a quantity which is 3/2 times larger in the pion than in the proton because of the number of available valence quarks.
In Fig. 16 we compare the x π dependence of the pion structure function from our global QCD analysis with the two models in which F π 2 is computed from the proton F p 2 structure function, both at the input scale Q 2 = µ 2 0 and evolved to Q 2 = 10 GeV 2 . For the simple scaling model we fix the proportionality constant at the value α = 0.65 fitted to give the best agreement at Q 2 = 10 GeV 2 , which happens to be similar to the 2/3 factor of the NSZ model. We also show the pion structure function from the analysis including only Drell-Yan p T -integrated data, F π,DY 2 . As the uncertainties on F π,DY 2 are very large at low x π , the inclusion of the LN data in our analysis is vital to obtain better constraint on the pion's structure function. We also note that since no data currently exist for pion structure functions at Q 2 = m 2 c , the F π 2 in Fig. 16 at the input scale involves an element of extrapolation.
Interestingly, at Q 2 = 10 GeV 2 the scaled proton structure function αF p 2 (x B , Q 2 ) follows closely the pion structure function at small values of x π 0.2 and lies within its uncertainties.
At large x π and x B , on the other hand, the pion F π 2 has a relatively slow falloff as x π → 1, while the proton F p 2 is much softer and falls faster as x B → 1. This effect simply reflects the different behaviors of the valence quark PDFs in the pion and proton. In contrast, for the NSZ model the large-x B behavior of the rescaled F p 2 is more consistent with the pion structure function, due to the 2/3 rescaling factor in the argument. However, at lower x B the agreement with the pion F π 2 is not as good in this model. These observations suggest that a hybrid model, in which the low-x B behavior (x B 0.2) is given by αF p 2 (x B , Q 2 ) and the high-x B behavior (x B 0.2) is given by 2 3 F p 2 ( 2 3 x B , Q 2 ), may be a better representation of F π 2 . Note that whatever agreement exists between the proton and pion structure functions at Q 2 = 10 GeV 2 does not translate to lower Q 2 values. At the input scale, Q 2 = µ 2 0 , for example, the shapes of F π 2 and the two models based on F p 2 differ substantially, with opposite slopes. In particular, at low values of x B 0.1 the latter are significantly larger than the F π 2 from the global analysis, with neither model able to accurately capture the low-x π behavior of F π 2 . The agreement between F π 2 and F p 2 for low x π and x B at Q 2 = 10 GeV 2 , and the disagreement at Q 2 = µ 2 0 , can be understood in terms of the flavor decomposition of the structure functions into valence, sea, and gluon components and the different roles played by each.
Firstly, since the valence quark PDFs in the pion and proton have rather different x dependence, the structure functions at large x π or x B are expected to differ. At low x values, the sea quark PDFs in the pion and proton are also quite different, with a nontrivial x dependence, as illustrated in Fig. 15. This is directly responsible for the difference between F π 2 and F p 2 at low x π and x B at the input scale. As Q 2 increases, the role of evolution becomes more important, and the contributions to the structure functions from gluons become more prominent, even though formally these enter at higher order in α s . Since gluon PDFs in the pion and proton are fairly similar at low values of x, as observed in Fig. 15, their contributions to the structure functions at low x π and x B must also be similar, which is reflected in the better agreement between the pion and proton F 2 functions at higher Q 2 values. Overall, however, our results would suggest that the structure functions at low x π and x B are strongly dependent on the nonperturbative sea quark distribution and on the Q 2 evolution of the PDFs, and it is difficult to draw general, scale-independent conclusions about their characteristics.

VII. OUTLOOK
One of the most important achievements of the present work has been the demonstration that it is indeed possible to describe, within the same O(α s ) collinear factorization framework, the p T -dependent spectrum of Drell-Yan lepton pairs in pion-nucleon collisions, along with the more traditional p T -integrated Drell-Yan data and LN electroproduction cross sections. To our knowledge, this is the first demonstration of its kind for the pion, and paves the way to the quantitative study of the three-dimensional parton structure of the pion.
The inclusion of the p T -dependent data contributes to the reduction of uncertainties on the gluon distribution at large x, although the impact of the existing data is relatively modest compared to the p T -integrated Drell-Yan and LN data, and calls for future, high-precision measurements of large-p T pion-induced lepton-pair production cross sections.
Future improvements to this work will implement threshold resummation in the perturbative calculations of p T -integrated Drell-Yan cross sections [37], which can produce nonnegligible enhancement of the cross section from the emission of soft gluons near threshold.
Previous work [65] indicated potentially significant impact on the large-x behavior of the valence pion PDF, and it will be important to explore the systematic uncertainties of the calculation [66][67][68]. Beyond this, future theoretical improvements in the analysis of p Tdependent data should incorporate p T smearing as well as higher-order gluon emissions, both of which affect the energy and kinematics of incoming partons.
On the experimental side, the planned program of pion-and kaon-induced Drell-Yan measurements at the future COMPASS++/AMBER [62] facility at CERN should provide precision data to better determine the gluon and sea quark PDFs in the pion, and allow first glimpses of the partonic structure of kaons. Additional p T -dependent Drell-Yan data at high p T may help to constrain the valence quark and gluon distributions at large values of x. Other planned experiments, such as the Tagged Deep-Inelastic Scattering (TDIS) experiment at Jefferson Lab, will provide complementary information to the HERA LN data by measuring a leading proton produced in DIS from a neutron in a deuteron. Utilizing the Sullivan process, this reaction will probe the structure of the exchanged pion in a kinematic region between the (low-x) HERA and (high-x) Drell-Yan domains. Experiments at the future Electron-Ion Collider (EIC) may also provide better statistics for LN measurements at low x π values in tagged electron-proton collisions.