Parton Distribution Functions of the Charged Pion Within The xFitter Framework

Using xFitter, we extract parton distribution functions of charged pions from currently available Drell-Yan and photon production data. While the valence distribution is well-constrained, we find that the considered data are not sensitive enough to unambiguously determine sea and gluon distributions. In the used approximation, we find the high-x behavior of the valence distribution to be linear in (1-x) at high x. Fractions of momentum carried by the valence, sea and gluon components are discussed.


INTRODUCTION
The pion plays an important role in our understanding of strong interactions. At the same time, it is a mediator of nucleon-nucleon interactions, a pseudo-Goldstone boson of dynamical chiral symmetry breaking and the simplest qq state in the quark-parton model of hadrons. However, from the experimental point of view, the pion structure is currently poorly understood, especially compared to the proton. Parton distribution functions (PDFs) are a primary theoretical construct used to describe hadron structure as it is probed in hard processes. Much progress has been made in mapping out the parton distribution functions of the proton in the last decades [1].
Experimentally, the pion PDF is known mostly from QCD analyses of Drell-Yan (DY) and prompt photon production data [18][19][20][21]. Within a dynamical approach, only the relatively well-known valence distribution is determined from DY data, with the sea and gluon content at a very low initial scale fixed by simplifying assumptions [22] or constraints of the constituent quark model [23,24]. While all modern pion PDF extractions are performed at next-to-leading order (NLO), additional threshold-resummation corrections and their impact on the valence distribution at high x have been studied [25]. In addition to DY data, a recent work by the JAM collaboration [26] included leading neutron electroproduction data obtained from the HERA collider (as suggested in [27]). The latest pion PDF fit by Bourelly and Soffer [28] uses a novel parameterisation at the initial scale Q 0 .
In this work we approach the pion PDFs from a phenomenological perspective, and extract them from the currently available experimental data on pion scattering using xFitter, an open-source PDF fitting framework [29].
The paper is organised as follows: In Section I we briefly discuss the considered data. The adopted PDF parameterisation and decomposition are described in Section II. Calculation of theoretical predictions is discussed in Section III. Section IV is devoted to the statistical treatment used in this work and estimation of the uncer-arXiv:2002.02902v1 [hep-ph] 7 Feb 2020 tainty of the obtained PDFs. Finally, the results of the analysis are presented and compared to results of other studies in Section V.

I. EXPERIMENTAL DATA
This analysis is based on Drell-Yan data from NA10 [30] and E615 [31] experiments, and on photon production data from the WA70 [32] experiment. The NA10 and E615 experiments studied scattering of a π − beam off a tungsten target, with E π = 194 and 286 GeV in the NA10 experiment and E π = 252 GeV in the E615 experiment. The WA70 experiment used π ± beams and a proton target. For the Drell-Yan data, Υ-resonance range, which corresponds to bins with √ τ ∈ [0.415, 0.484], were excluded from the analysis. Here √ τ = m µµ / √ s, m µµ is the invariant mass of the muon pair, and √ s is the center-of-mass energy of pion-nucleon system. Leading order Feynman diagrams for the considered processes are shown in Fig. 1. The Drell-Yan data constrain the valence distribution relatively well, but are not sensitive to sea and gluon distributions. The prompt photon production data complement the DY data by providing some sensitivity to the gluon distribution, but have smaller statistics and large uncertainties in comparison to the DY data. Additionally, the predictions for prompt photon production have significant theoretical uncertainty, as discussed in Section III.

II. PDF PARAMETERISATION
The π − PDF xf (x, Q 2 ) is parameterized at an initial scale Q 2 0 = 1.9 GeV 2 , just below the charm mass threshold m 2 c = 2.04 GeV 2 . Neglecting electroweak corrections and quark masses, charge symmetry is assumed: d =ū, and SU(3)-symmetric sea: u =d = s =s. Under these assumptions, pion PDFs are reduced to three distributions: total valence v, total sea S, and gluon g: which we parameterise using a generic form: where B is the Euler beta function, which ensures that the A S parameter represents the total momentum fraction carried by the sea quarks. The B-parameters determine the low-x behavior, and C-parameters determine the high-x behavior. Quark-counting and momentum sum rules have the following form for π − : The sum rules determine the values of parameters A v and A g , respectively. The constant factors in the definitions of v, S, g were chosen in such a way, that xv , xS , xg are momentum fractions of pion carried by the valence quarks, sea quarks, and gluons, respectively (here xf = 1 0 xf (x)dx). The extension D v x α was introduced in xv(x) to mitigate possible bias due to inflexibility of the chosen parameterisation. This extension was omitted in the initial fits (D v = 0). Afterwards, a parameterisation scan was performed by repeating the fit with free D v and different fixed values of parameter α. The scan showed that only α = 5 2 has noticeably improved the quality of the fit (see Table I and Section V for discussion). The additional free parameter D v changes the shape of the valence distribution only slightly (Fig. 2). Similar attempts to add more parameters of the form (1 + D v x α + E v x β ) did not result in significant improvement of χ 2 . The final presented results use a free D v and α = 5 2 .
The valence distribution when using minimal parameterisation (Dv = 0) and the extended parameterisation with free Dv. The shown uncertainty bands do not include scale variations. The high-x behavior is linear in (1 − x).

III. CROSS-SECTION CALCULATION
PDFs are evolved up from the starting scale Q 2 0 by solving the DGLAP equations numerically using QCDNUM [33]. The evolution is performed using the variable flavornumber scheme with quark mass thresholds at m c = 1.43 GeV, m b = 4.5 GeV. Predictions for the crosssections were calculated as a convolution of the evolved pion PDFs with precomputed grids of NLO coefficients and with PDFs of a proton or tungsten target. The APPLgrid [34] package was used for these calculations. The grids were generated using the MCFM [35] generator. For Drell-Yan, the invariant mass of the lepton pair was used for the renormalization and factorisation scales, namely µ R = µ F = m ll . For prompt photon production, the scale was chosen as the transverse momentum of the prompt photon, namely µ R = µ F = p T (γ).
It was verified that the grid binning was sufficiently fine by comparing the convolution of the grid with the PDFs used for the grid generation and a reference cross-section produced by MCFM. The deviation from the reference cross-section, as well as estimated statistical uncertainty of the predictions, are an order of magnitude smaller than the uncertainty of the data. This check was performed for each data bin.
Both the evolution and cross-section calculations are performed at next-to-leading order (NLO). For the tungsten target, nuclear PDFs from nCTEQ15 [36] determination were used. In the case of a proton target, the PDFs from ref. [37] were employed. These were also used as the baseline in the nCTEQ15 study. The use of another popular nuclear PDF set EPPS16 [38] was omitted because their fit had used the same piontungsten DY data as the present analysis. Considering π − N data, EPPS16 fitted PDFs of tungsten using fixed pion PDFs from an old analysis by GRV [22].
In the case of prompt photon production, the contribution of fragmentation photons cannot be accounted for using the described techniques. The model used in the fit included only the direct photons. We estimate the impact of the missing fragmentation contribution by comparing the total integrated cross-sections computed using MCFM for proton-proton collision at the WA70 energy with and without fragmentation. The relative difference of 32% is treated as the theoretical uncertainty in overall normalization of the WA70 data. In the run without fragmentation, Frixione isolation is used. In the other run the fragmentation function set GdRG LO and cone isolation are used. The isolation cone size parameter is R 0 = 0.4 for both cases.

IV. STATISTICAL TREATMENT AND ESTIMATION OF UNCERTAINTIES
The PDF parameters are found by minimizing the χ 2 function defined as where i is the index of the datapoint and α is the index of the source of correlated error. The measured cross-section is denoted by d i , with δ syst i and δ stat i being respectively the corresponding systematic and statistical uncertainties. The t i 's represent the calculated theory predictions, andt i = t i (1 − α γ iα b α ) are theory predictions corrected for the correlated shifts. γ iα is the relative coefficient of the influence of the correlated error source α on the datapoint i, and b α is the nuisance parameter for the correlated error source α.
The error rescalingδ stat = t i di δ stat is used to correct for Poisson fluctuations of the data. Since statistical uncertainties are typically estimated as a square root of the number of events, a random statistical fluctuation down in the number of observed events leads to a smaller estimated uncertainty, which gives such points a disproportionately large weight in the fit. The error rescaling corrects for this effect. This correction was only used for the Drell-Yan data.
The nuisance parameters b α are used to account for correlated uncertainties. In this analysis the correlated uncertainties consist of the overall normalization uncertainties of the datasets, the correlated shifts in predictions related to uncertainties from nuclear PDFs, and the strong coupling constant α S (M 2 Z ) = 0.118 ± 0.001. The nuisance parameters are included in the minimization along with the PDF parameters. They determine shifts of the theory predictions and contribute to the χ 2 via the penalty term α b 2 α . For overall data normalization, the coefficients γ iα are relative uncertainties as reported by the corresponding experiments, and, in the case of the WA70 data, the abovementioned additional 32% theoretical uncertainty, (listed in Table II). For the uncertainties from nuclear PDFs and α S , the coefficients γ iα are estimated as derivatives of the theory predictions with respect to α S and the uncertainty eigenvectors of the nuclear PDFs as provided by the nCTEQ15 set. This TABLE II. The normalization and partial χ 2 for the considered datasets. The normalization uncertainty is presented as estimated by corresponding experiments. In order to agree with theory predictions, the measurements must be multiplied by the normalization factor. Deviations from 1 in the normalization factor lead to a penalty in χ 2 , as described in Section IV.

Experiment
Normalization uncertainty linear approximation is valid only when the minimisation parameters are close to their optimal values. It was verified that this condition was satisfied for the performed fits.
The uncertainty of the perturbative calculation is estimated by varying the renormalization scale µ R and factorization scale µ F by a factor of two up and down, separately for µ R and µ F . The scales were varied using APPLgrid, and the variations were coherent for all data bins. Renormalization scale variation for DGLAP evolution was not performed. We observe a significant dependence of the predicted cross-sections on µ R and µ F : the change in predictions is ∼ 10%, which is comparable to the normalization uncertainty of the data. This dependence indicates that next-to-next-to-leading order corrections may be significant.
In order to estimate the uncertainty related to the flexibility of chosen parameterisation, the fit is repeated with a varied initial scale Q 2 0 = 1.9 ± 0.4 GeV 2 . This variation leads to only a small change in χ 2 (∆χ 2 < ∼ 1). In order to stay below the charm mass, for variation up to Q 2 0 = 2.3 GeV 2 the mass threshold m 2 c was shifted up by the same amount. The effect of such a change in the charm mass threshold by itself was found to be negligible. Figure 3 shows the obtained pion PDFs in comparison to a recent analysis by JAM [26], and to GRVPI1 [22] -the only set available in the LHAPDF6 [39] library. The new valence distribution presented here is in good agreement with JAM, and both disagree with the early GRV analysis. The relatively difficult to determine sea and gluon distributions are different in all three PDF sets, however, this new PDF and the JAM determination agree within the larger uncertainties of our fit.

V. RESULTS
In the case of valence distribution, the dominant contribution to the uncertainty estimate is the variation of the scales µ R and µ F . For the sea and gluon distributions, the missing fragmentation contribution to prompt photon production is the dominant uncertainty source, and the effect of scale variation is also significant.
A comparison between experimental data and theory predictions obtained with the fitted PDFs is presented in Fig. 5. Reasonable agreement between data and theory is observed, with no systematic trends for any of the kinematic regions.
The extracted value of parameter C v is consistent with unity, meaning that v(x) ∼ (1 − x) as x → 1. However, the data are only sensitive to some region x < 1 and do not constrain the derivative of v(x) at x = 1. Although the addition of the extra parameter D v changes the fitted value of C v , the valence distribution is still linear near x ∼ 1 in the experimentally accessible high-x region (Fig. 2). The change in C v shows that the derivative of v(x) at x = 1 cannot be extracted from the considered data and that the behavior of v(x) can be studied only in some region x < 1. The behavior v(x) ∼ (1 − x) is favored by Nambu-Jona-Lasinio models [40]. At the same time it is in conflict with approaches based on the Dyson-Schwinger equations (DSE) [8,9], which predict v(x) ∼ (1 − x) 2 . This discrepancy between DSE predictions and fits to pion Drell-Yan data is well-known [9,21,25]. It has been shown that soft-gluon threshold resummation -which was not included in this analysis -may be used to account for this disagreement [25]. Alternatively, DSE calculations using inhomogeneous Bethe-Salpeter equations [9] produce PDFs consistent with the linear behavior of the v(x) in the region covered by DY data, pushing the onset of the (1 − x) 2 regime to very high x. Fig.  4 shows the obtained momentum fractions in the pion as a function of Q 2 in comparison to the proton momentum fractions calculated using the NNPDF31 nlo as 0118 [41] set. The points correspond to other studies, which are also listed in Table III. The sum of momentum fractions is forced to one according to Eq.(1) by adjusting the gluon normalization parameter A g . In comparison to the proton, in the pion the valence quarks carry a larger momentum fraction. Above the charm and bottom mass thresholds Q > m c , m b , the c and b quarks and antiquarks are included in the sea distribution.

VI. SUMMARY AND OUTLOOK
PDFs of charged pion have been determined by reanalysing the currently available Drell-Yan and prompt photon production data using modern tools. While the valence distribution is well constrained, the considered data are not sensitive enough to unambiguously determine the sea and gluon distributions. While the data are reasonably well-described by NLO QCD, the sensitivity to µ R and µ F indicates that next-to-nextto-leading order corrections could be significant. The valence distribution behaves as v(x) ∼ (1 − x) as x → 1 in the experimentally accessible region, although the considered data do not constrain the derivative at x = 1. The valence momentum fraction in the pion is found to be large in comparison to the proton. In the future, new FIG. 3. Comparison between the pion PDFs obtained in this work, a recent determination by the JAM collaboration [26], and the GRVPI1 pion PDF set [22].
FIG. 4. Momentum fractions of the pion as a function of Q 2 . The error bands include all uncertainties described in Section IV. Analogous momentum fractions in the proton PDF set NNPDF31 nlo as 0118 are shown for comparison. The labeled green, red, and blue points show respectively valence, sea, and gluon momentum fractions as reported by other studies. The references and numerical values for these points are listed in Table III. data from the COMPASS++/AMBER [43] experiment may allow for more stringent constraints of pion PDFs.
It is planned to publish the extracted pion PDFs in the LHAPDF6 PDFs library and the APPLgrid grid files in the Ploughshare [44] grid library.  [21] 0.434 ± 0.022 27 ChQM-1 [11] 0.428 27 ChQM-2 [13] 0.46 27 this work 0.42 ± 0.04 0.25 ± 0.10 0.32 ± 0.10 27 SMRS [20] 0 Considered experimental data and corresponding theory predictions. The displayed theory predictions include correlated shifts. Bands of different colors correspond to different datasets. Width of the bands shows uncertainty of the theory predictions. The cross-sections are shown in the same format as adopted by corresponding experimental papers. The E615 data is given as d 2 σ/(d √ τ dxF ) in nb/nucleon, averaged over each ( √ τ , xF ) bin. The DY data from the NA10 experiment is d 2 σ/(d √ τ dxF ) in nb/nucleus, integrated over each ( √ τ , xF ) bin. The WA70 data on direct photon production is given as invariant cross-section Ed 3 σ/dp 3 in pb, averaged over each (pT , xF ) bin.