Gluon Parton Distribution of the Pion from Lattice QCD

We present the first determination of the $x$-dependent pion gluon distribution from lattice QCD using the pseudo-PDF approach. We use lattice ensembles with 2+1+1 flavors of highly improved staggered quarks (HISQ), generated by MILC Collaboration, at two lattice spacings $a\approx 0.12$ and 0.15~fm and three pion masses $M_\pi\approx 220$, 310 and 690 MeV. We use clover fermions for the valence action and momentum smearing to achieve pion boost momentum up to 2.29 GeV. We find that the dependence of the pion gluon parton distribution on lattice spacing and pion mass is mild. We compare our results from the lightest pion mass ensemble with the determination by JAM and xFitter global fits.


I. INTRODUCTION
The lightest bound state in quantum chromodynamics (QCD), the pion, plays a fundamental role, since it is the Nambu-Goldstone boson of dynamical chiral symmetry breaking (DCSB). Studies of pion and kaon structure reveal the physics of DCSB, help to reveal the relative impact of DCSB versus the chiral symmetry breaking by the quark masses, and are important to understand nonperturbative QCD. Studying the pion parton distribution functions (PDFs) is important to characterize the structure of the pion and further understand DCSB and nonperturbative QCD. Currently, our knowledge of the pion PDFs is less than the nucleon PDFs, because there are fewer experimental data sets, especially for the seaquark and gluon distributions. The future U.S.-based Electron-Ion Collider (EIC) [1], planned to be built at Brookhaven National Lab, will further our knowledge of pion structure [2,3]. In China, a similar machine, the Electron-Ion Collider in China (EicC) [4], is also planned to make impacts on the pion gluon and sea-quark distributions. In Europe, the Drell-Yan and J/ψ-production experiments from COMPASS++/AMBER [5] will aim at improving our knowledge of both the pion gluon and quark PDFs.
Global analyses of pion PDFs mostly rely on Drell-Yan data. The early studies of pion PDFs were based mostly on pion-induced Drell-Yan data and use J/ψ-production data or direct photon production to constrain the pion gluon PDF [6][7][8][9][10]. There are more recent studies, such as the work by Bourrely and Soffer [11], that extract the pion PDF based on Drell-Yan π + W data. JAM Collaboration [12,13] uses a Monte-Carlo approach to analyze the Drell-Yan πA and leading-neutron electroproduction data from HERA to reach the lower-x region, and revealed that gluons carry a significantly higher momentum fraction (about 30%) in the pion than had been inferred from Drell-Yan data alone. The xFitter group [14] analyzed Drell-Yan πA and photoproduction data using their open-source QCD fit framework for PDF extraction and found that these data can constrain the valence distribution well but are not sensitive enough for the sea and gluon distributions to be precisely determined. The analysis done Ref. [15] suggests that the pion-induced J/ψ-production data has additional constraint on pion PDFs, particularly in the pion gluon PDF in the large-x region. All in all, the pion valence-quark distributions are better constrained than the gluon distribution from the global analysis of experimental data. While waiting for more experimental data sets, the study of the pion gluon distribution from theoretical side can provide useful information for the experiments.
The pion gluon PDF is rarely studied using continuum-QCD phenomenological models or through lattice-QCD (LQCD) simulations. Most model studies only predict the pion valence-quark distribution [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30], but the gluon and sea PDFs are predicted by the Dyson-Schwinger equation (DSE) continuum approach [31]. The prediction of the pion gluon PDF in DSE, based on an implementation of rainbow-ladder truncation of DSE, is consistent with the JAM pion gluon PDF result [12,13] within two sigma. LQCD provides an first-principles calculations to improve our knowledge of nonperturbative pion gluon structure; however, there have been only two efforts to determine the first moment of pion gluon PDF [32,33]. An early calculation in 2000 using quenched QCD predicted x g = 0.37(8)(12), using Wilson fermion action with a lattice spacing a = 0.093 fm, lattice size L 3 × T = 24 4 , a large 890-MeV pion mass and 3,066 configurations at µ 2 = 4 GeV 2 [32]. A more recent study in 2018 using N f = 2 + 1 clover fermion action with a lattice spacing a = 0.1167(16) fm, larger lattice size 32 3 × 96, 450-MeV pion mass, and 572,663 measurements, gave a larger first-moment result, x g = 0.61 (9) at µ 2 = 4 GeV 2 [33]. In principle, a series of moments can be used to reconstruct the PDF. Although there are calculations of the first moment of the pion gluon PDF, there is little chance that sufficient higher moments of the pion gluon PDF can be obtained to perform such a reconstruction. A direct lattice calculation of the xdependence of the pion gluon PDF is needed.
In recent years, there has been an increasing number of calculations of x-dependent hadron structure in lattice QCD, following the proposal of Large-Momentum Effective Theory (LaMET) [34][35][36]. The LaMET method calculates lattice quasi-distribution functions, defined in terms of matrix elements of equal-time and spatially separated operators, and then takes the infinite-momentum limit to extract the lightcone distribution. The quasi-PDF can be related to the P z -independent lightcone PDF through a factorization theorem that factors from it a perturbative matching coefficient with corrections suppressed by the hadron momentum [35]. The factorization can be calculated exactly in perturbation theory [37,38]. Many lattice works have been done on nucleon and meson PDFs, and generalized parton distributions (GPDs) based on the quasi-PDF approach . Alternative approaches to lightcone PDFs in lattice QCD are "good lattice cross sections" [37,[72][73][74][75] and the pseudo-PDF approach [76][77][78][79][80][81][82][83][84][85][86][87][88][89]. However, LQCD calculations of the x dependence of the pion PDFs have only been done for the valence-quark distribution [49,64,74,75,81,90,91].
Only recently have lattice calculations of the gluon PDF become possible, when the necessary one-loop matching relations of the gluon PDF were computed for the pseudo-PDF [92] and quasi-PDF [54,93] approaches. Both approaches make direct calculation of the x dependence of the pion gluon PDF feasible. In this work, we apply the pseudo-PDF method by using the ratio renormalization scheme to avoid the difficulty of calculating the gluon renormalization factors. There is a developed procedure for using the pseudo-PDF method to obtain lightcone PDFs from Ioffe-time distributions (ITDs) by matching through two steps, evolution and scheme conversion [82]. Using the pseudo-PDF method also allows us to use lattice correlators at all boost momenta at small Ioffe-time. There have been a number of successful pseudo-PDF calculations of nucleon isovector PDFs [76,80,85,86] and pion valence-quark PDFs [81]. The earliest calculation was done on a quenched lattice [76], then the pion masses were set closer to the physical pion mass [80,81,85], and the calculation at physical pion mass was done recently [86]. The lattice-calculated PDFs in Refs. [80,81,85,86] show good agreement with the global-analysis PDFs.
In this work, we present the first calculation of the full x-dependent pion gluon distribution using the pseudo-PDF method from two lattice spacings, 0.12 and 0.15 fm, and three pion masses: 690, 310 and 220 MeV. The rest of the paper is organized as follows. In Sec. II, we present the procedure to obtain the lightcone gluon PDF from the reduced ITDs, the numerical setup of lattice simulation, and how we extracted the reduced ITDs from lattice calculated correlators. In Sec. III, the final determination of the pion gluon PDF from our lattice calculations is compared with the NLO xFitter [14] and JAM pion gluon PDFs [12,13]. The systematics induced by different steps are studied, and the lattice-spacing and pion-mass dependence are investigated.

II. GLUON PDF FROM LATTICE CALCULATION USING PSEUDO-PDF METHOD
In this work, we use the unpolarized gluon operator defined in Ref. [92], , z is the Wilson link length, and the field strength F µν is defined as where the a is the lattice spacing, g 0 is the strong coupling constant, and the plaquette P µ, There is an alternative operator i =z,t O(F ti , F zi ; z) corresponding to the same matching kernel in Ref. [92]. We do not choose this operator, because it vanishes at P z = 0 for kinematic reasons, bringing additional difficulties in obtaining the distributions.
Using this operator in Eq. 1, we calculate lattice gluon matrix elements of the ground-state meson |0(P z ) with various boost momenta P z and Wilson-line displacement lengths z. We then study their dependence on Ioffe time ν = zP z , calling M(ν, z 2 ) the Ioffe-time distribution (ITD). To eliminate the ultraviolet divergences in the ITD, we construct the reduced ITD (RITD) by taking the ratio of the ITD to the corresponding z-dependent matrix element at P z = 0, and further normalize the ratio by the matrix element at z 2 = 0 as done in the first quark pseudo-PDF calculation [76], The renormalization of O(z) and kinematic factors are cancelled out in the RITDs. The RITD double ratios used here are automatically normalized to one at z = 0. The RITDs are related to the pion gluon g and quark q S PDFs via the pseudo-PDF matching condition [92] M (ν, where µ is the renormalization scale in MS scheme and x g = 1 0 dx xg(x, µ 2 ) is the gluon momentum fraction of the pion. We can split the gluon-in-gluon R gg contribution in Eq. 5 into two parts, which are introduced in Eqs. 21 and 24 in Ref. [82], where R 1 (y, z 2 µ 2 ) is the term related to evolution, R 2 (y) is the term related to scheme conversion, α s is the strong coupling at scale µ, N c = 3 is the number of colors, and γ E = 0.5772 is the Euler-Mascheroni constant. The z in R 1 (y, z 2 µ 2 ) is chosen to be 2e −γ E −1/2 /µ so that the log term vanishes, suppressing residuals that contain higher orders of the log term, as discussed in the paper on the one-loop evolution of the pseudo-PDF [82]. The gluon-inquark kernel R gq (y, z 2 µ 2 ), along with R B (y), R L (y) and R C (y), are defined in Eqs. 7.21-23 and the paragraph below Eq. 7.23 in Ref. [92].
In this work, we first neglect the pion quark PDF, since the total quark PDF is found to be much smaller than the gluon PDF in global fits [12,14]. We will later estimate the systematic uncertainty introduced by this assumption. The gluon evolved ITD (EITD), G is obtained by using the evolution term R 1 (y, z 2 µ 2 ), The z dependence of the EITDs should be compensated by the ln z 2 term in the evolution formula. In principle, the EITD G is free of z dependence and is connected to the lightcone gluon PDF g(x, µ 2 ) through the schemeconversion term R 2 (y), so the gluon PDF g(x, µ 2 ) can be extracted by inverting this equation.
On the lattice, we use clover valence fermions on three ensembles with N f = 2 + 1 + 1 highly improved staggered quarks (HISQ) [94] generated by the MILC Collaboration [95] with two different lattice spacings (a ≈ 0.12 and 0.15 fm) and three pion masses (220, 310, 690 MeV). The masses of the clover quarks are tuned to reproduce the lightest light and strange sea pseudoscalar meson masses used by PNDME Collaboration [96][97][98][99]. We use five HYP-smearing [100] steps on the gluon loops to reduce the statistical uncertainties, as studied in Ref. [52]. We use Gaussian momentum smearing for the quark fields [101] to reach higher meson boost momenta. Table I gives the lattice spacing a, valence pion mass M val π and η s mass M val ηs , lattice size L 3 × T , number of configurations N cfg , number of total two-point correlator measurements N 2pt meas , and separation time t sep used in the three-point correlator fits for the three ensembles. This allows us to reach the continuum limit and physical pion mass through extrapolation. The total amount of measurements vary in 10 5 -10 6 for different ensembles.  and ηs mass M val ηs , lattice size L 3 × T , number of configurations N cfg , number of total two-point correlator measurements N 2pt meas , and separation times tsep used in the three-point correlator fits of N f = 2 + 1 + 1 clover valence fermions on HISQ ensembles generated by MILC Collaboration and analyzed in this study.
The two-point correlator for a meson Φ is where P z is the meson momentum in the z-direction, χ Φ =q 1 γ 5 q 2 is the pseudoscalar-meson interpolation operator, t is the Euclidean time, and |A Φ,i | 2 and E Φ,i are the amplitude and energy for the ground-state (i = 0) and the first excited state (i = 1), respectively. The three-point gluon correlators are obtained by combining the gluon loop with pion two-point correlators. The matrix elements of the gluon operators can be obtained by fitting the three-point correlators to the energyeigenstate expansion, where t sep is the source-sink time separation, and t is the gluon-operator insertion time. The amplitudes and energies, A Φ,0 , A Φ,1 , E Φ,0 and E Φ,1 , are obtained from the two-state fits of the two-point correlators. 0|O|0 , 0|O|1 ( 1|O|0 ), and 1|O|1 are the ground-state matrix element, the ground-excited-state matrix element, and the excited-state matrix element, respectively. We shown on all plots is the extracted ground-state matrix element from the two-sim fit using tsep ∈ [5,9]. From left to right, the columns are: the ratio of the three-point to two-point correlators with the reconstructed fit bands from the two-sim fit using tsep ∈ [5,9], shown as functions of t − tsep/2, the one-state fit results for the three-point correlators at each tsep ∈ [3,9], the two-sim fit results using tsep ∈ [t min sep , 9] as functions of t min sep , and the two-sim fit results using tsep ∈ [5, t max sep ] as functions of t max sep .
extract the ground-state matrix element 0|O|0 from the two-state fit of the three-point correlators, or a two-state simultaneous "two-sim" fit on multiple separation times with the 0|O|0 , 0|O|1 and 1|O|0 terms.
To verify that our fitted matrix elements are reliably extracted, we compare to ratios of the three-point to the two-point correlator if there were no excited states, the ratio would be the ground-state matrix element. The left-hand side of Fig. 1 shows example ratios for the gluon matrix elements from the lightest pion ensemble, a12m220, at selected momenta P z and Wilson-line length z. We see the ratios increase with increasing source-sink separation going from 0.60 to 1.08 fm. At large separation, the ratios begin to converge, indicating the neglect of excited states becomes less problematic. The gray bands indicate the groundstate matrix elements extracted using the two-sim fit to three-point correlators at five t sep . The convergence of the fits that neglect excited states can also be seen in second column of Fig. 1, where we compare one-state fits from each source-sink separations: the one-state fit results increase as t sep increases, starting to converge at large t sep to the two-sim fit results. The third and fourth columns of Fig. 1 show two-sim fits using t sep ∈ [t min sep , 9] and t sep ∈ [5, t max sep ] to study how the two-sim ground-state matrix elements depend on the source-sink separations input into fit. We observe that the matrix elements are consistent with each other within one standard deviation, showing consistent extraction of the ground-state matrix element, though the statistical errors are larger than those of the one-state fits. We observe larger fluctuations in the matrix element extractions when small t min sep = 3 and 4, or small t max sep = 6 and 7, are used. The ground state matrix element extracted from two-sim fits becomes very stable when t min sep > 4 and t max sep > 7. Figure 2 shows the RITD of the same examples P z = 2 × 2π/L, z = 1 and P z = 4 × 2π/L, z = 4 from two-sim fit results using t sep ∈ [t min sep , 9]. The RITD results, which are constructed to suppress lattice fluctuations, are very stable over the range of different fits considered. For a12m310 and a15m310 ensembles, the t sep dependence of RITDs is milder than those from a12m220 ensemble due to the heavier pion mass. Overall, our ground-state RITDs from the two-sim fit are stable, and we use them to extract the gluon PDF.

III. RESULTS AND DISCUSSIONS
Using the RITDs extracted in the previous section, we examine the pion-mass and lattice-spacing dependence. The top of Fig. 3 shows the η s RITDs at boost momentum around 2 GeV as functions of the Wilson-line length z for the a12m220, a12m310, and a15m310 ensembles. We see no noticeable lattice-spacing dependence. The bottom of Fig. 3 shows the pion RITDs with boost momentum around 1.3 GeV for the same ensembles. Again, there is no visible lattice-spacing or pion-mass dependence.
To extract gluon PDFs, we follow the steps in Sec. II between Eq. 3 and Eq. 9 by first obtaining EITDs and using Eq. 10 to extract g(x). To obtain EITDs, we need the RITD M (ν, z 2 ) to be a continuous function of ν to evaluate the x ∈ [0, 1] integral in Eq. 9. We achieve this by using a "z-expansion" 1 fit [102,103] to the RITD. The following form is used [81]: where τ = √ νcut+ν− √ νcut √ νcut+ν+ √ νcut . Then, we use the fitted M (ν, z 2 ) in the integral in Eq. 9. The z-dependence in the M (uν, z 2 ) term of the evolution function comes from the one-loop matching term, which is a higher-order correction compared to the tree-level term; thus, the zdependence can be neglected in M (ν, z 2 ). We adopt as a12m220, P z =1.96 GeV a12m310, P z =2.14 GeV the best value ν cut = 1, as used in Ref. [81], but we also vary ν cut in the range [0. 5,2], and the results are consistent. We fix λ 0 = 1 to enforce the RITD M (ν, z 2 ) in Eq. 4. The expansion order k max = 3 is used, because we can fit all the data points of P z ∈ [1, 5] × 2π/L (P z ∈ [1, 7] × 2π/L for a12m220 ensemble) and z up to 0.6 fm with χ 2 /dof < 1 using a 4-term z-expansion for each ensemble. The reconstructed bands from "zexpansion" on RITDs are shown in the upper plot in Fig. 4. They describe the RITD data points well for all ensembles.
After we have the continuous-ν fitted RITDs, we obtain the EITDs through Eq. 9. The RITDs M and EITDs G as functions of ν on all ensembles studied in this work are shown in Fig. 4. At some ν values, there are multiple z and P z combinations for a fixed ν value. Therefore, there are points in the same color and symbol overlapping at the same ν from the same lattice spacing and pion mass. To match with the lightcone gluon PDF through Eq. 10, the EITDs G(ν, µ) should be free of z 2 dependence. However, the EITDs obtained from Eq. 9 have z 2 dependence from neglecting the gluon-inquark contribution and higher-order terms in the matching. The EITDs also depend on lattice-spacing a and pion-mass M π . Recall that the RITDs show weak dependence on lattice spacing a and pion mass M π . We see  that the effects of a and M π dependence on the EITDs are also not large; the EITD results from different a, M π are mostly consistent with each other, as shown in the second row of Fig. 4. We also observe a weak dependence on z 2 for the RITDs and EITDs in Fig. 4.
The gluon PDF g(x, µ 2 ) can now be extracted from the EITDs using Eq. 10. We assume a functional form, also used by JAM [12,13], for the lightcone PDF to fit the EITD, for x ∈ [0, 1] and zero elsewhere. The beta function B(A + 1, C + 1) = 1 0 dx x A (1 − x) C is used to normalize the area to unity. Then, we apply the matching formula to obtain the EITD G from the functional form PDF using Eq. 10. We fit the EITDs G(ν, µ) obtained from the parametrization to the EITDs G(ν, z 2 , µ, a, M π ) from the lattice calculation.
We investigate the systematic uncertainty introduced by the different parametrization forms which are commonly used for f g (x, µ) in PDF global analysis and some lattice calculations. The first one is the 2-parameter form in Eq. 15. Second, we consider the 1-parameter form N 1 (1 − x) C used in xFitter's analysis [14] (also used in Ref. [7,8]), which is equivalent to Eq. 15 with A = 0. Third, we consider a 3-parameter form We fit the three different forms to the EITDs of lattice data with z max ≈ 0.6 fm by applying the scheme conversion Eq. 10 to the 1-, 2-and 3-parameter PDF forms. Here, we focus on the result from the lightest pion mass M π ≈ 220 MeV at lattice spacing a ≈ 0.12 fm. The χ 2 /dof of the fits decreases as 1.47(72), 1.08(68), to 1.04(41), shows slightly better fit quality for 2-and 3parameter fits. As shown in Fig. 5, there is a big discrepancy between the f g (x, µ) fit bands from the 1-parameter fit and the 2-parameter fit in the x < 0.4 region, but the discrepancy between the 2-and 3-parameter fits is much smaller. Therefore, we conclude that 1-parameter fit on lattice data here is not quite reliable, and the fit results converge at the 2-and 3-parameter fits. The same conclusions hold for all other ensembles and pion masses. Therefore, using the 2-parameter form defined Eq. 15 (same parametrization as JAM) for our final results is very reasonable.
Another source of systematic uncertainty comes from neglecting the contribution of the quark term in Eq. 5 based on the assumption (motivated by global fits) that the pion q S (x) is smaller than the gluon PDF. Currently, there are no q S (x) results from lattice simulation since only the valence distribution of the pion has been done. Thus, we estimate the systematic due to omitting the q S (x) contribution by using the pion quark PDFs from xFitter [14] at NLO. Using these, we obtain revised RITDs and EITDs including the gluon-in-quark R gq term focusing on example from the a ≈ 0.12 fm, pion mass M π ≈ 220 MeV lattice, repeating the same procedure from Eq. 14 and fitting the EITDs with Eq. 15. On the right-hand side of Fig. 6, we show the mean value of xg(x, µ)/ x g with both gluon-in-gluon (gg) and gluonin-quark (gq) contributions (the blue solid line) compared to the a12m220 results using the gluon-in-gluon contribution only (the green solid line). There are 5 to 10% differences in the mean value including the gluon-in-gluon contribution for x < 0.9, which indicates that the gluonin-quark contribution is relatively small at µ 2 = 4 GeV 2 compared to the current statistical errors in the small-x region. In the x > 0.9 region, the gluon-in-quark contribution becomes more significant, but it remains smaller than the statistical error. Once studies are available with sufficiently reduced statistical uncertainty in the large-x region, the quark contribution will need to be included.
From the above analyses of the choice of fit form and the contribution of the quark term, we conclude that these systematics are negligible relative to the current statistics. Therefore, we adopt the z max ≈ 0.6 fm (z max ≈ 0.75 fm for a15m310 ensembles) fits to the EITDs, neglect the quark contribution term in the matching, and use the Eq. 10 fit form for our final results on all lattice ensembles. The xg(x, µ)/ x g reconstructed fit bands of these ensembles are shown in the left plot in Fig. 6, comparing results from different lattice spacings and pion masses. The reconstructed fit bands with different pion mass M π ≈ {220, 310, 690} MeV are consistent at the same lattice spacing a ≈ 0.12 fm, indicating mild gluon PDF dependence on pion mass. Similarly, when comparing lattice-spacing dependence of pion PDFs using data around pion mass M π ≈ 310 MeV, we find that fitted PDF is slightly smaller in the x > 0.1 region for the 0.12-fm lattice, but still within one sigma, which indicates the lattice-spacing dependence is also mild. We also note that the bands from different ensembles show a differing speed of fall-off as x → 1 in the large-x region. We study this fall-off behavior in more depth below.
The behavior of the gluon PDF fall-off in the large-x region is widely studied in both theory and global analyses. Perturbative QCD studies [104,105] and DSE calculations [29,31] suggest that the gluon distribution g(x, µ 2 ) ∼ (1 − x) C with C ≈ 3 in the limit x → 1. The prediction from perturbative QCD [105] is based on the idea that the gluon PDF should be suppressed at large x relative to the quark PDF, because the quarks are the sources of large-x gluons; that is, g(x, µ 2 )/q v (x, µ 2 ) → 0 as x → 1. Early fits of experimental data gave C ≈ 2 [7,8] or C < 2 [9,10], but the more recent global analysis from JAM collaboration yielded C > 3 [12,13] and xFitter collaboration found C ≈ 3 [14]. Our fitted parameter C is 3.6(1.5), 3.3(2.0), 4.7(2.8) for M π ≈ {690, 310, 220} MeV, respectively, at lattice spacing a ≈ 0.12 fm. These C results are consistent with each other and show a slightly increasing trend as the pion mass approaches the physical pion mass. For lattice spacings a ≈ {0.15, 0.12} fm, C = {2.2(1.5), 3.3(2.0)}, respectively, at M π ≈ 310 MeV, which suggests that C will increase toward the continuum limit. We also investigate the effect of the gluon-in-quark contribution on the C value, and it makes about 0.1 difference, which we neglect. Given that both the pion-mass and lattice-spacing extrapolations seem to show increasing C, it seems rea-sonable to conclude from this lattice-QCD study that C > 3.
We compare our reconstructed gluon PDF to those from global fits on the right-hand side of Fig. 6. It shows the xg(x, µ)/ x g reconstructed fit band of a ≈ 0.12 fm, M π ≈ 220 MeV lattice and the NLO pion gluon PDFs from xFitter [14] and JAM [12,13] at µ 2 = 4 GeV 2 . The JAM band appears somewhat wider than expected, because we reconstruct it by dividing xg(x, µ) by the mean value of x g ; the correlated values needed for a correct error estimation were not available. Note that xFitter uses the fit form of Eq. 15 with A = 0. Our fitted pion gluon PDF is consistent with JAM for x > 0.2 and with xFitter for x > 0.5 within one sigma. We see in this comparison that our results are of similar error size as the global-fit analysis and are useful to provide constraints from theoretical calculation in addition to the experimental data.

IV. SUMMARY
In this work, we presented the first calculation of the pion gluon PDF from lattice QCD and studied its pionmass and lattice-spacing dependence using the pseudo-PDF approach. We employed clover valence fermions on ensembles with N f = 2 + 1 + 1 highly improved staggered quarks (HISQ) at two lattice spacings (a ≈ 0.12 and 0.15 fm) and three pion masses (220, 310 and 690 MeV). These ensembles allowed us to probe the dependence of the pion gluon PDF on pion mass and lattice spacing. In both cases, the dependence appears to be weak compared to the current statistical uncertainty.
We investigated the systematics associated with the functional form used in the reconstruction fits as well as the systematics caused by neglecting the quark contribution in the matching. The effect of the assumed gluon PDF fit form was investigated by using various forms, which are all commonly used or proposed in other PDF works. We observe large effects changing the fit to xg(x, µ)/ x g from 1-to 2-parameter form but convergence at 3 parameters. This implies the 2-parameter fits are sufficient for our calculation, and our finial pion gluon PDF results are presented using the 2-parameter fit results. We used the pion quark PDF from xFitter to make an estimation of the quark contribution to the pion gluon RITD. We found the systematic errors it contributed are smaller than 10% of the statistical errors.
Our pion gluon PDF for the lightest pion mass is consistent with JAM for x > 0.2 and with xFitter for x > 0.5 within uncertainty, as shown in our final comparison plots of the pion gluon PDF. We also studied the asymptotic behavior of the pion gluon PDF in the large-x region in terms of (1−x) C . C > 3 is implied from our study at two lattice spacings and three pion masses. The future study of the pion gluon PDF from the lattice QCD with improved precision and systematic control when combined in global-fit analyses with the results of anticipated ex-  6. The pion gluon PDF xg(x, µ)/ x g as a function of x obtained from the fit to the lattice data on ensembles with lattice spacing a ≈ {0.12, 0.15} fm, pion masses Mπ ≈ {220, 310, 690} MeV (left), and xg(x, µ)/ x g as function of x obtained from lattices of a ≈ 0.12 fm, Mπ ≈ 220 MeV (right), compared with the NLO pion gluon PDFs from xFitter and JAM at µ = 2 GeV in the MS scheme. The JAM error shown is overestimated due to lack of available correlated uncertainties in its constituent components. Our PDF results are consistent with JAM [12,13] for x > 0.2 and xFitter [14] for x > 0.5.
periments [2,4,5] will provide best determination of the gluon content within the pion.