Gluon distributions and their applications to Ioffe-time distribution

We investigate unpolarized and polarized gluon distributions and their applications to the Ioffe-time distributions, which are related to lattice QCD calculations of parton distribution functions. Guided by the counting rules based on the perturbative QCD at large momentum fraction $x$ and the color coherence of gluon couplings at small $x$, we parametrize gluon distributions in the helicity basis. By fitting the unpolarized gluon distribution, the inferred polarized gluon distribution from our parametrization agrees with the one from global analysis. A simultaneous fit to both unpolarized and polarized gluon distributions is also performed to explore the model uncertainty. The agreement with the global analysis supports the $(1-x)$ power suppression of the helicity-antialigned distribution relative to the helicity-aligned distribution. The corresponding Ioffe-time distributions and their asymptotic expansions are calculated from the gluon distributions. Our results of the Ioffe-time distributions can provide guidance to the extrapolation of lattice QCD data to the region lacking precise gluonic matrix elements. Therefore, they can help regulate the ill-posed inverse problem associated with extracting the gluon distributions from discrete data from first-principle calculations, which are available in a limited range of the nucleon momentum and the spatial separation between the gluonic currents. The possibility of investigating higher-twist effects and other systematic uncertainties in the contemporary first-principle calculations of parton distributions from phenomenologically well-determined Ioffe-time distributions in the large Ioffe-time region is also discussed.


I. INTRODUCTION
One of the outstanding problems in nuclear and particle physics is to understand the structure of hadrons in terms of quarks and gluons, the fundamental degrees of freedom in QCD. Gluons, which serve as mediator bosons of the strong interaction while also carrying the color charge, play a key role in the nucleon's mass and spin. With great progress in the extraction of nucleon parton distribution functions (PDFs) in the past decades, especially the quark distributions, the understanding of the gluon distribution and its role in hadron structures remains one of the most challenging but fundamental issues in nuclear and particle physics.
Given the fact that no free quarks or gluons have been observed due to the confinement, most analyses of hadron involved high energy scatterings rely on the QCD factorization, where PDFs play an important role. Compared to quark distributions, the gluon distribution is less accurately extracted, which may affect the calculation of the cross section of a process dominated by the gluonfusion channel, e.g. the Higgs boson production at the LHC [1]. While the precision of the extracted g(x) has been improved a lot during the last decade, there are still some issues like the suppression in the momentum fraction region 0.1 < x < 0.4 when ATLAS and CMS jet data are included [2]. Obtaining a more precise determination of g(x) is subject to ongoing efforts in global analyses of PDFs. In contrast to the unpolarized PDFs, the polarized PDFs, especially the polarized gluon dis-tribution ∆g(x) as well as sea quark distributions, are poorly determined, even the sign is not fully determined. One of the main physics goals of the upcoming Electron-Ion-Collider (EIC) [3] is to have precise measurements of the nucleon spin structure, particularly the gluon and sea quark distributions.
It has been 30 years since the EMC experiment [4], which found that only a small fraction of the nucleon spin is carried by the quark spin and triggered the socalled proton spin puzzle. The remaining part of the proton spin is usually assigned to the orbital angular momenta and the gluon spin. After significant efforts in the last decades, the quark spin part was found to contribute only about 30% to the proton spin. Recent experimental data from RHIC and lattice QCD calculation suggest that gluon may contribute a large amount to the proton spin. For a recent review, see [5]. There have been several global analyses using different experimental data sets and different types of parametrizations [6][7][8][9] to impose constraints on ∆G. A recent extraction with updated data sets and PHENIX measurement [10] of double helicity asymmetry in inclusive π 0 production in polarized p − p collision obtained ∆G = 0.2 with a constraint of −0.7 < ∆G < 0.5 for the sampled gluon momentum fraction x ∈ [0.02, 0.3]. Excluding the x < 0.05 region, the value of ∆G = 0.2 0.05 dx ∆g(x) = 0.23 (6) [11] and ∆G = 1 0.05 dx ∆g(x) = 0.19 (6) [12] were obtained. Future experimental measurement of ∆g(x) in the x < 0.02 is required to reduce the uncertainty in ∆G. We note that ∆g(x) extracted mostly from the double longitudi-2 nal spin asymmetry is always limited to some x min no matter how high the energy of the experimental set up is and some theoretical calculations are required to constrain ∆g(x) at low x [13]. Fortunately, one of the major goals of the upcoming EIC [3] is to precisely explore ∆g(x) at low x and provide stringent constraints on the gluon helicity distribution.
Since ∆G is not related to the local matrix element of a gauge invariant operator, it could not be directly calculated in the lattice QCD (LQCD) calculations. However, following the formalism proposed in [14], it has become possible to calculate ∆G in terms of a local and gauge invariant operator in LQCD. Since then, there has only been one direct LQCD calculation [15,16] of the gluon spin content in the nucleon. With leading-order matching using Ji's large-momentum effective theory (LaMET) [17,18], it was determined that ∆G(µ 2 = 10GeV 2 ) = 0.251 (47) (16), i.e. about 50% of the proton spin comes from the gluons. However, a refined study with the investigation of the convergence of matching beyond 1-loop, the estimate of power corrections, and other sources of systematic uncertainties are warranted to obtain an unambiguous determination of ∆G in the future LQCD calculations. Now, the PDFs of the quarks and gluons contain the nonperturbative structure of hadrons, especially at the low-resolution scale Q 2 . However, the possibility that even at low Q 2 one can obtain a reasonable shape and distribution of the hadron structure functions by transcribing our knowledge of the perturbative QCD (pQCD) based counting rules [19] as x → 1 and color coherence of gluon couplings as x → 0 has been shown to provide promising outcomes in many theoretical calculations. For example, calculations of the unpolarized and polarized quark and gluon PDFs in [20,21] showed practical application of these limiting behaviors of PDFs by obtaining PDFs in agreement with the analysis in [22]. The momentum fraction carried by the gluon in the nucleon x g ≈ 0.42 determined in [21] is in remarkable agreement with recent global analyses [23][24][25]. To emphasize further, most of the earlier and present-day global analyses [23,24,26,27] use some functional forms similar to x α (1 − x) β F(x) where the asymptotic behavior of the PDFs at small x (x α behavior) is adopted from the observed Regge behavior [28] in particle colliders and the large-x behavior ((1 − x) β fall-off) based on the power counting rules for hard scattering [19] with some interpolating function F(a i , b i , · · · , x) with unknown parameters a i , b i , · · · between these two limits that varies in different parametrizations of PDFs. Once these PDFs are determined at some initial scale, their Q 2 -evolution is well predicted in pQCD through the DGLAP equation [29][30][31]. These PDFs determined at the low initial scale have been shown to be universal between different reactions with their scale-dependent modifications governed by pQCD evolution, and therefore indicated that these nonperturbative universal PDFs can indeed be well approximated even at low Q 2 by incorporating pQCD constraints at large x and Regge behavior at small x. Very good agreement and consistent behavior of the nucleon unpolarized PDFs and precise prediction of nucleon polarized distributions from the unpolarized PDFs have also been possible in recent calculations [32,33] where the PDFs are governed by these limiting behaviors. Similarly, recent synergies between LQCD and phenomenological calculations have provided useful constraints in the study sea-quark asymmetry in the nucleon with higher precision than either theory or experiment alone could attain [34,35].
In light of the above discussions, we revisit the calculations in [20,21] which incorporated pQCD constraints at large x and coherent correlations of partons at low x to determine unpolarized and polarized gluon distributions. These calculations [20,21] demonstrated that many properties of the exclusive reactions can be calculated by incorporating the knowledge of asymptotic freedom, power-law scaling, and helicity conservation rules of pQCD without explicit knowledge of the nonperturbative light-front wave function.
The main goal of this article is to transcribe these insights from the small and large x physics and compare how adequate and compatible they are with the recent determinations of gluon distributions. We emphasize this calculation does not aim to provide the precise determination of the gluon PDFs, rather our main focus is to determine the shapes of the gluon PDFs based on these simple constraints in the small and large x-regions. Here we point out that, we do not focus here on the important aspects of gluon distributions at extremely small x-values which has been discussed in the literature [36][37][38][39][40][41][42]. Another important goal of the upcoming EIC is to explore the very low-x region where saturation of gluon densities sets in [43,44] and has not yet been conclusively observed. Parton distributions at extremely smallx is an active field of research and we avoid the discussion of the related complication here. We first determine the unknown coefficients in the parametrization of helicity aligned g + (x) and anti-aligned g − (x) gluon distributions using the global fits of unpolarized gluon distribution and use those to calculate the polarized gluon distribution and gluon asymmetry distribution ∆g(x)/g(x). We calculate the corresponding Ioffe-time distributions (ITDs) [45][46][47] of the unpolarized and polarized gluon distributions and demonstrate how these can provide valuable information and important constraints in the determination of full x-dependence of PDFs and also their higher moments in the future LQCD calculations. In particular, we determine the asymptotic behavior of the unpolarized and polarized gluon ITDs which are not accessible within the current reach of LQCD calculations and can provide complementary information to reconstruct the full x-dependence of the unpolarized and polarized gluon distributions.

II. GLUON DISTRIBUTIONS FROM HELICITY-BASIS PARAMETRIZATION
To construct the parametrization of the helicity-basis gluon distributions, g + (x) for helicity-aligned distributions and g − (x) for helicity-antialigned distribution, we consider the counting rules based on perturbative QCD analysis [21]. Compared to valence quark distributions, which fall off as (1 − x) 3 as x → 1, g + (x) is suggested to fall-off faster as (1 − x) 4 and g − (x) is expected to be further suppressed by (1 − x) 2 . Instead of strictly imposing the power counting as x → 1 and the Pomeron intercept as x → 0, we only take them as guidance and phenomenologically introduce two parameters α and β to allow the variation of the power behavior at small and large x regions as usually adopted in global analyses. For a good description of the gluon distribution in the fullx region, we also include a polynomial (1 + γ √ x + δx) with parameters γ and δ to be fitted. As a modification of the functional form utilized in [21] by including the polynomial, we parametrize the helicity-aligned and the helicity-antialigned gluon distributions as where A and B are normalization parameters to be determined. The inclusion of the subleading term in the power of (1 − x) is to account for the contribution from higher Fock state. For each term, the power of (1 − x) differs by 2 as suggested by the pQCD analysis [21]. We refer to the parametrization form Eq. (1) as the ansatz-1.
As a phenomenological exploration, we consider another option of g − (x) being suppressed by one power of (1 − x) in comparison with the g + (x). This results in the parametrization, which we refer to as the ansatz-2. As we will discuss later, fixing the (1 − x) power difference between g + (x) and g − (x) introduces a model bias, which leads to an underestimation of the uncertainties. To investigate the model uncertainty, we consider a more flexible parametrization, where the (1 − x) exponents and the polynomial coefficients in g + (x) and g − (x) are independent parameters.
We refer to this parametrization as ansatz-3. We note that all these ansatzes have the ∆g(x) approaching to 0 as x → 0. This indicates that the helicity correlation between the gluon and its parent nucleon disappears when x → 0, where the relative rapidity becomes infinity. The saturation effect may suppress the evolution of helicity distributions at small-x and consequently leave a small amount of spin contribution in the small-x region [13,48].
Since the goal of this paper is not the small-x distribution, we limit to the assumption above in this study and restrict the subsequent analyses in the x ≥ 10 −3 region. With the parametrization of the helicity-aligned and the helicity-antialigned gluon distributions, one can directly obtain the unpolarized and polarized gluon distributions from the sum and the difference of them, To determine the parameters in ansatz-1 and ansatz-2, we fit the unpolarized gluon distribution from the NNPDF global analysis [25] at the factorization scale µ = 2 GeV. Our procedure described here can be applied to any other gluon distribution given by global analyses [23,24,26,27] or model calculation. To fit the distribution in the full-x range, we select 200 points in x values. 100 of them are equally separated from 10 −4 to 10 −1 in the logarithmic scale and the other 100 x-values are equally separated from 10 −1 to 1 in the linear scale. Each point is weighted by the inverse square of its uncertainty given by the global analysis. We take the 100 replicas of the gluon distribution from NNPDF3.1 NLO PDF set [25]. For each replica, we perform a fit to determine the parameters. In the end, we have 100 sets of parameters, which determine the gluon distributions following the ansatz-1 or the ansatz-2. For the ansatz-3 in which the parameters β, γ and δ are chosen independently for g + (x) and g − (x), we perform a simultaneous fit to the unpolarized [25] and polarized [11] gluon distributions. In this case, the result of the polarized gluon distribution is driven by the global fit. As a result, the polarized gluon distribution associated with ansatz-3 has a better match with the NNPDF global analysis. The results of the unpolarized gluon distribution are compared with the global analysis in FIG. 1, where the central value is evaluated from the average value of the 100 replicas for each ansatz and the uncertainty band is the standard deviation among them. One can observe that the three ansatzes have almost indifferentiable results and match the global analysis well.
From the definition of the polarized gluon distribution in Eq. (5), we now obtain the polarized gluon distribution based on the above fit results of xg + (x) and xg − (x). Unlike the unpolarized distribution, the results of ∆g(x) determined from ansatz-1, 2, and 3 have observable difference, especially in the region 10 −2 < ∼ x < ∼ 0.5, as shown in FIG. 2. However, all these determinations of x∆g(x) are in good agreement within the uncertainty range of the NNPDF global analysis with a noticeable difference between the NNPDF and ansatz-1 distributions in the 0.09 < ∼ x < ∼ 0.2 region. We note that the small uncertainties of ∆g(x) from ansatz-1 and ansatz-2 are biased by the parametrization form, Eqs (1) and (2), where the (1 − x) power difference are fixed between g + (x) and g − (x). For the ansatz-3, we introduce independent parameters for the (1 − x) powers and the polynomial parts of the two helicity basis distributions. Thus, the ansatz-3 is a more flexible parametrization than ansatz-1 and ansatz-2, but it requires a simultaneous fit to both unpolarized and polarized distributions to determine the parameters. Therefore, the result from ansatz-3 is driven by global analysis and less biased. The difference between the result from ansatz-3 and the one from ansatz-1 or ansatz-2 indicates the model uncertainty of imposing the (1 − x) power difference of g + (x) and g − (x), or, in other words, how much the uncertainties of the results from ansatz-1 and ansatz-2 are biased.
Due to the current precision of experimental data, the phenomenological determination of ∆G is sensitive to the parametrization form in the global analysis. If allowing a possible sign-change of ∆g(x) at some x value, one will find large uncertainties of ∆g(x) and thus very poor constraint on ∆G. In our approach, the helicity retention is incorporated in our parametrization of ansatz-1 and ansatz-2, where the polarized gluon distribution is fixed once the unpolarized distribution is determined. As we will show in the next section, the result from the ansatz-3 is also consistent with the helicity retention, although it is not imposed in the parametrization form.
One can observe that the uncertainties of the x∆g(x) 10 -3 10 -2 10 -1 10 0 Polarized gluon distributions from the fit parameters determined from fitting ansatz-1 and ansatz-2 to the NNPDF unpolarized gluon distribution. Ansatz-3 refers to the polarized gluon distribution obtained from a simultaneous fit to the NNPDF3.1 NLO PDF set [25] and NNPDFpol1.1 PDF set [11] using fit paramterization in Eq. (3). The gray band shows the polarized gluon distribution x∆g(x) as given by the NNPDFpol1.1 global analysis [11]. The blue, red, and cyan bands labeled by Ansatz-1, Ansatz-2, and Ansatz-3 show distributions obtained using parameters obtained in the fits of xg + (x) and xg − (x) to the NNPDF gluon distribution.
determined from ansatz-1 and ansatz-2 are highly constrained. This is due to the bias of the parametrization form, which assumes a relation between the two helicitybasis distributions and thus leads to an underestimation of the uncertainties of the polarized distribution. On the other hand, ansatz-3 is more flexible and the uncertainty of x∆g(x) is governed by the global analysis of x∆g(x). An outstanding question is how to distinguish between these three different determinations of x∆g(x) distributions, especially in the large x-region which is of primary interest for the nonperturbative LQCD calculations of PDFs. One answer is, as we will see in Section IV, the gluon helicity ∆G obtained from the Ioffe-time distribution obtained from ansatz-1 parametrization has a magnitude almost twice as large compared to the one resulting from ansatz-2. Similarly, the Ioffe-time distribution of the polarized gluon distribution obtained from ansatz-1 is almost double in magnitude compared to that obtained from ansatz-2. The difference of these two Ioffetime distributions in the small Ioffe-time region ω ≈ 0−6 can be investigated in LQCD calculations to discriminate between these two ansätze. On the other hand, the polarized gluon distribution determined by fitting the NNPDFpol1.1 global analysis [11] is data-driven and has much larger uncertainty compared to that obtained from ansatz-1 and ansatz-2. Therefore, the resulting polarized ITD also has larger uncertainty. It will be a good opportunity to explore the polarized gluon ITD in LQCD calculations and have a possible impact on constraining the uncertainty. We will present detailed discussion of this prospect in Sections IV and IV.

III. GLUON ASYMMETRY DISTRIBUTION
The COMPASS experiment at CERN has measured gluon asymmetry distribution ∆g(x)/g(x) from the cross section helicity asymmetry of photon-gluon fusion (γ * g → qq) in the semi-inclusive deep-inelastic scattering (DIS) of proton-proton collision [49][50][51]. Although the open charm events provide the cleanest signal to the γ * g → qq (q = c) events [52,53], the rate of these events is very small. The high statistics two-jets events with large transverse momentum p T with respect to the virtual photon direction can give access to the photon-gluon fusion subprocess but with a price of significant background which has to be subtracted in a model-dependent way to determine ∆g(x)/g(x).
In FIG. 3, we compare the ∆g(x)/g(x) ratio obtained from our calculation with data at different x-values extracted from high p T hadrons in the leading-order analyses [49,50] and from the open charm production in the next-to-leading order analysis [51] at COMPASS, from high p T hadrons in the leading-order analyses by the Spin Muon Collaboration (SMC) at CERN [54] and at the HERMES experiment [55]. We note that the endpoint values of ∆g(x)/g(x) are fixed in ansatz-1 and ansatz-2. In the limit x → 0 the ratio goes to 0 and as x → 1 the ratio goes to 1, no matter what values are assigned to the parameters in Eqs. (1) and (2). The difference between the results from ansatz-1 and ansatz-2 may be regarded as the model uncertainty, while the uncertainty for either of them is model biased. It is not a surprise to find that the result from ansatz-3 has much larger uncertainty, because the parametrization of ansatz-3 is more flexible and thus less biased. One can also notice that the ratio does not necessarily go to 1 in the limit x → 1 for the ansatz-3, since we introduce independent parameters for the (1 − x) powers of g + (x) and g − without the requirement of any (1 − x) power suppression of g − (x) in comparison with g + (x). However, the result from the simultaneous fit is still consistent with 1 at the x → 1 endpoint.
The two different solutions of ∆g(x) obtained in the COMPASS next-to-leading order analyses [56,57] allow ∆g(x)/g(x) to be positive or negative in the entire xregion, whereas our analyses with ansatz-1 and ansatz-2 give a positive ∆g(x)/g(x) in the entire x-region and a very small ∆g(x)/g(x) in the x < 10 −1 region. The x → 1 value of this asymmetry distribution is consistent with the pQCD prediction of the helicity retention [21,58]. Even for the ansatz-3, where we do not require the helicity retention in the parametrization, the simultaneous fit result driven by the global analysis is still consistent with this prediction.
Comparison between the two determinations of ∆g(x)/g(x) from this calculation with the experimental measurements. The direct measurements of COMPASS [49,50], HERMES [55], and SMC [54] are obtained in leading order from high pT hadrons and from open charm muon production at COMPASS [51] in next-to-leading order at different xvalues are shown. The blue, red, and cyan bands labeled by Ansatz-1, Ansatz-2, and Ansatz-3 show the gluon asymmetry distributions determined using the parameters obtained using ansatz-1, ansatz-2, and ansatz-3 for the xg + (x) and xg − (x) distributions, respectively as discussed in the main text.

IV. GLUON IOFFE-TIME DISTRIBUTIONS
As first discussed by Gribov, Ioffe, and Pomeranchuck [45], the large coherent length distances defined in the target rest frame become important at high energies for the virtual photon-nucleon scattering. Ioffe further demonstrated a connection between the DIS scattering amplitude and the spacetime representation of the correlator of the electromagnetic current [46], establishing a relation between longitudinal coherent distances and Bjorken scaling. Application of the Ioffe-time distribution (ITD) to study parton distributions in coordinate space using nonperturbative method was proposed in [47]. It was argued in [47] that calculations of PDFs in momentum space receive contributions from both small and large longitudinal distances for each value of x and can result in a problem of treating different physics associated with different distances simultaneously. The Ioffetime τ I measures the interval between the absorption and emission of virtual-photon by a hadron in DIS and gives the coherence length of the pair production in the target rest frame, where M is the target mass. Braun et. al. entitled the Lorentz invariant variable ω as the Ioffe-time; and we will use this naming convention for the remainder of this article, where the same language is seen in contempo-6 rary coordinate space LQCD formalisms used to isolate parton distributions [59][60][61]. We note ω = M τ I is defined in the hadron's rest frame, hence the designation "time" despite "ω" itself being neither time nor space in LQCD calculations. It is indeed the ω-dependence of the ITD that converts into the x-dependence of the parton distributions. Recently, a method for obtaining frameindependent, three-dimensional light-front coordinatespace wave functions and its relevance to LQCD calculations of PDFs has also been discussed in terms of frame-independent longitudinal distance (the Ioffe-time) in [62].
One can now write the unpolarized gluon distribution in terms of its Ioffe-time distribution as [47,63] M(ω, at a scale µ 2 . Similarly using the definition of the polarized gluon distribution from [64], x∆g(x, µ 2 ) is related to its Ioffe-time distribution [47,63] via In comparison with Eq. (7), the sin(xω) in the integrand of Eq. (8) leads to one additional power of x suppression when x is small and therefore reduce the small-x region contribution, as well as its uncertainty, to the ∆M(ω, µ 2 ). With knowledge of the polarized gluon ITD, from Eq. (8), one can immediately obtain the gluon helicity contribution to the nucleon spin where we have used the principal value prescription to calculate the integral ∞ 0 dω sin(xω). As seen from Eq. (9), and as will be explored further in Section VI, we shall see that access to the asymptotic region of the ITD ∆M(ω, µ 2 ), namely up to ω ≈ 15 can provide a stringent constraint on the gluon helicity in the nucleon. Given the gauge-invariant and frame-independent definition of the ITD, one can take advantages of calculating the ∆G(µ 2 ) from the ITD to avoid the issues in the spindecomposition [65]. Using Eqs. (7) and (8), we calculate the ITDs of the unpolarized and polarized gluon distributions from the parameterizations of g + (x) and g − (x) and present them in FIGs. 4 and 5, respectively.
Similar to the unpolarized gluon PDF, we see from FIG. 4 that ansätze 1, 2, and 3 produce almost identical ITDs for the xg(x) distributions. However, we see that the magnitude of the x∆g(x) ITD (FIG. 5) in the ω ∼ 5−10 range is almost double when using the ansatz-1 parametrization relative to the ansatz-2 parametrization. Consequently, using Eq. (9), we obtain ∆G = 0.451 (7) from ansatz-1 and ∆G = 0.258(4) from ansatz-2. As mentioned above, the difference between the ∆G values from ansatz-1 and ansatz-2 may be viewed as model uncertainty, while the uncertainty band of either one is biased by the parametrization form. The result from the simultaneous fit with ansatz-3 is ∆G = 0.23 (41). Such a large uncertainty from ansatz-3 indicates the fact that 7 the ∆G is still very poorly known. We note that the lower value of ∆G = 0.258(4) obtained from ansatz-2 is consistent with the previous LQCD determination [16] of ∆G. The most recent calculation of the nucleon spin decomposition at the physical pion mass [66] found the total gluon angular momentum contribution in the proton to be 0.187 (47). According to this calculation, unless the gluons contribute a large and negative orbital angular momentum to the nucleon total angular momentum budget, the ∆G contribution to the proton spin is expected to be less than ∼ 50%. On the other hand, if one excludes the low x contribution of the polarized gluon PDF one obtains ∆G ∼ 40% for x > 0.02 in [11] and for x > 0.05 in [12]. Of course, these values can change in global fits, if gluons are shown to have large positive contributions to the nucleon spin -a prospect to be explored at the EIC [3]. Another possibility is the x∆g(x) distribution can change sign in the low-x region, thereby reducing the contribution of ∆G to the nucleon spin budget.
One important question we encountered in Sec. II was how one can discriminate between ansatz-1 and ansatz-2 and the resulting ∆g(x) distributions. An interesting feature we observe from FIGs. 4 and 5 is the significant difference in the magnitude between these two ITDs in the ω ∼ 0−5 range. This will provide a great opportunity for LQCD calculations to discriminate between these two different determinations for the x∆g(x) distribution. On the other hand, the polarized gluon ITD obtained from ansatz-3 has a much larger uncertainty. It remains to see that if LQCD calculation in the lower Ioffe-time region can be precise enough to provide complementary information about the uncertainty of the polarized gluon ITD. We highlight recent LQCD calculations [67][68][69][70] have obtained precise ITD data in the ω < ∼ 5 region. Resolution to this problem in future LQCD calculations of the polarized gluon ITD hinges on precise LQCD data for ω ∼ 0 − 5, as well as mitigation of higher-twist effects which in turn ensures the validity of the short-distance pQCD factorization of lattice matrix elements into PDFs.

V. CALCULATION OF HIGHER MOMENTS FROM GLUON IOFFE-TIME DISTRIBUTIONS
From the ITDs in Eqs. (7) and (8), one can immediately calculate the associated higher moments of the distributions. Performing a power series expansion of cos(xω) in Eq. (7) (with scale µ 2 omitted ), we can get access to any number of moments depending on the available ω-range of ITDs. It is clearly seen from FIG. 6 that one can reproduce the ITDs with increasingly better accuracy and in the larger ω-region in terms of higher moments. Another way to state this observation is that, access to ITD in increasingly larger ω range would result in access to increasingly higher number of moments. However, even without access to ITD in the region ω → ∞, it is still possible to obtain PDFs with very good accuracy (again assuming the Regge phenomenology is valid in approximately x > 10 −2 region). As we will see in the next Section that access to ITD in the region ω ≈ 0 − 15 for the gluon distributions can be sufficient to extract their x-dependence in the 10 −2 > x > 1 region.
unpolarized and polarized gluon ITDs which are accessible within the reach of present LQCD calc tions and can provide complementary information reconstruct full x-dependence of the unpolarized polarized gluon distributions in the future calculatio The helicity-aligned gluon distribution is expected fall of at least (1 x) faster than the valence quark tribution in the nucleon [? ] as x ! 1. The helic aligned gluon distribution g + (x) is thus expected to h leading (1 x) 4 fall-o↵ as x ! 1 compared to (1 fall-o↵ of the valence quark distribution. The he ity anti-aligned gluon distribution g (x) is expected be suppressed by an additional factor of (1 x) 2 ative to g + (x) as x ! 1. We note that non-lead helicity flip amplitudes and contributions from hig Fock components are power-law suppressed and we o include the next higher Fock component contribution x ! 1 in our calculation. Besides, we also allow varia in the exponents of (1 x) fall-o↵ through the param and instead of fixing the pomeron intercept to exactly equal to 1, we also allow variation in the expon ↵. Instead of strictly imposing the power counting x ! 1 and the pomeron intercept at x ! 0, we o take them as a guidance and phenomenologically in duce two parameters ↵ and to allow the variation of power behavior at the small and large x regions as u ally adopted in global analyses. Since it is expected t the shape of PDFs is not just completely governed by small and large x limiting behaviors and being motiva by the parametrizations in the global fits of PDFs, also allow a somewhat flexible interpolating func unpolarized and polarized gluon ITDs which are not accessible within the reach of present LQCD calculations and can provide complementary information to reconstruct full x-dependence of the unpolarized and polarized gluon distributions in the future calculations.

II. HELICITY-DEPENDENT GLUON DISTRIBUTION FROM HELICITY-INDEPENDENT GLUON DISTRIBUTION
The helicity-aligned gluon distribution is expected to fall of at least (1 x) faster than the valence quark distribution in the nucleon [? ] as x ! 1. The helicityaligned gluon distribution g + (x) is thus expected to have leading (1 x) 4 fall-o↵ as x ! 1 compared to (1 x) 3 fall-o↵ of the valence quark distribution. The helicity anti-aligned gluon distribution g (x) is expected to be suppressed by an additional factor of (1 x) 2 relative to g + (x) as x ! 1. We note that non-leading helicity flip amplitudes and contributions from higher Fock components are power-law suppressed and we only include the next higher Fock component contribution as x ! 1 in our calculation. Besides, we also allow variation in the exponents of (1 x) fall-o↵ through the parameter and instead of fixing the pomeron intercept to be exactly equal to 1, we also allow variation in the exponent ↵. Instead of strictly imposing the power counting at x ! 1 and the pomeron intercept at x ! 0, we only take them as a guidance and phenomenologically introduce two parameters ↵ and to allow the variation of the power behavior at the small and large x regions as usually adopted in global analyses. Since it is expected that the shape of PDFs is not just completely governed by the small and large x limiting behaviors and being motivated by the parametrizations in the global fits of PDFs, we also allow a somewhat flexible interpolating function The helicity-aligned gluon distribution is fall of at least (1 x) faster than the valenc tribution in the nucleon [? ] as x ! 1. T aligned gluon distribution g + (x) is thus expe leading (1 x) 4 fall-o↵ as x ! 1 compared fall-o↵ of the valence quark distribution. ity anti-aligned gluon distribution g (x) is be suppressed by an additional factor of ( ative to g + (x) as x ! 1. We note that helicity flip amplitudes and contributions Fock components are power-law suppressed include the next higher Fock component con x ! 1 in our calculation. Besides, we also all in the exponents of (1 x) fall-o↵ through th and instead of fixing the pomeron inte exactly equal to 1, we also allow variation in t ↵. Instead of strictly imposing the power x ! 1 and the pomeron intercept at x ! take them as a guidance and phenomenolog duce two parameters ↵ and to allow the var power behavior at the small and large x reg ally adopted in global analyses. Since it is ex the shape of PDFs is not just completely gove small and large x limiting behaviors and bein by the parametrizations in the global fits o also allow a somewhat flexible interpolati polarized gluon distributions in the future The helicity-aligned gluon distribution i fall of at least (1 x) faster than the valen tribution in the nucleon [? ] as x ! 1. aligned gluon distribution g + (x) is thus exp leading (1 x) 4 fall-o↵ as x ! 1 compare fall-o↵ of the valence quark distribution ity anti-aligned gluon distribution g (x) i be suppressed by an additional factor of ative to g + (x) as x ! 1. We note tha helicity flip amplitudes and contributions Fock components are power-law suppressed include the next higher Fock component co x ! 1 in our calculation. Besides, we also a in the exponents of (1 x) fall-o↵ through t and instead of fixing the pomeron int exactly equal to 1, we also allow variation in ↵. Instead of strictly imposing the powe x ! 1 and the pomeron intercept at x ! take them as a guidance and phenomenolo duce two parameters ↵ and to allow the va power behavior at the small and large x re ally adopted in global analyses. Since it is the shape of PDFs is not just completely go small and large x limiting behaviors and be by the parametrizations in the global fits also allow a somewhat flexible interpola Similarly, from Eq. (8), one can express the polarized gluon ITD in terms of odd moments as From Eqs. (9) and (11) it is seen that a calculation of ITD of polarized gluon distributions in ω ∼ 15 region will not only give a reliable estimate for gluon helicity ∆G in the nucleon but also several other higher moments which are currently unknown in theoretical calculations and not constrained by the experimental data. We see from FIGs. 6 and 7 that with increasing number of moments according to the Taylor expansion in Eqs. (10) and (11) (5) g h xi (7) g h xi (9) g The helicity-aligned gluon distribution is expected to fall of at least (1 x) faster than the valence quark distribution in the nucleon [? ] as x ! 1. The helicityaligned gluon distribution g + (x) is thus expected to have leading (1 x) 4 fall-o↵ as x ! 1 compared to (1 x) 3 fall-o↵ of the valence quark distribution. The helicity anti-aligned gluon distribution g (x) is expected to be suppressed by an additional factor of (1 x) 2 relative to g + (x) as x ! 1. We note that non-leading helicity flip amplitudes and contributions from higher Fock components are power-law suppressed and we only include the next higher Fock component contribution as x ! 1 in our calculation. Besides, we also allow variation in the exponents of (1 x) fall-o↵ through the parameter and instead of fixing the pomeron intercept to be exactly equal to 1, we also allow variation in the exponent ↵. Instead of strictly imposing the power counting at x ! 1 and the pomeron intercept at x ! 0, we only take them as a guidance and phenomenologically introduce two parameters ↵ and to allow the variation of the power behavior at the small and large x regions as usually adopted in global analyses. Since it is expected that the shape of PDFs is not just completely governed by the small and large x limiting behaviors and being motivated by the parametrizations in the global fits of PDFs, we also allow a somewhat flexible interpolating function where A and B are normalization parameters to be determined. The inclusion of the subleading term in power of (1 x) is to account for the contribution from higher Fock state. For each term, the power of (1 x) di↵ers by 2 as suggested by the pQCD analysis []. We refer to the parametrization form Eq. (??) as the ansatz-1. As a phenomenological exploration, the power di↵erence suggested by the pQCD power counting may also be modified. To take into account the model dependence on imposing the (1 x) power di↵erence between g + (x) and g (x), we consider another parametrization by assuming the power di↵erence is 1 as We also consider the following possibility when the large-x distribution of g (x) is suppressed by one power of (1 x) relative to the g + (x) distribution.
(Tianbo, can you provide a better reasoning here? some careful discussion needed for the choice of ansatzes 1 and 2) which we refer to as the ansatz-2. The polarized gluon distribution g(x) decreases to 0 as x ! 0 for both ansatz-1 and ansatz-2. This indicates the helicity correlation between the gluon and its parent nucleon disappears when x ! 0, which is a natural expectation since the relative rapidity becomes infinity.
With the parametrization of the spin-aligned and the spin-antialigned gluon distributions, one can directly obtain the unpolarized and polarized gluon distributions from the sum and the di↵erence of them, The unpolarized gluon distribution xg(x) can be written as a sum of the helicity aligned and anti-aligned glouon distributions To determine the unknown parameters in g +, (x) in Eqs. (??) and (??), we fit the unpolarized gluon distribution from the NNPDF global analysis []. Our precedure described here can be applied to any other gluon distribution given by global analysis or model calculation. To fit The helicity-aligned gluon distribution is expected to fall of at least (1 x) faster than the valence quark distribution in the nucleon [? ] as x ! 1. The helicityaligned gluon distribution g + (x) is thus expected to have leading (1 x) 4 fall-o↵ as x ! 1 compared to (1 x) 3 fall-o↵ of the valence quark distribution. The helicity anti-aligned gluon distribution g (x) is expected to be suppressed by an additional factor of (1 x) 2 relative to g + (x) as x ! 1. We note that non-leading helicity flip amplitudes and contributions from higher Fock components are power-law suppressed and we only include the next higher Fock component contribution as x ! 1 in our calculation. Besides, we also allow variation in the exponents of (1 x) fall-o↵ through the parameter and instead of fixing the pomeron intercept to be exactly equal to 1, we also allow variation in the exponent ↵. Instead of strictly imposing the power counting at x ! 1 and the pomeron intercept at x ! 0, we only take them as a guidance and phenomenologically introduce two parameters ↵ and to allow the variation of the power behavior at the small and large x regions as usually adopted in global analyses. Since it is expected that the shape of PDFs is not just completely governed by the small and large x limiting behaviors and being motivated by the parametrizations in the global fits of PDFs, we also allow a somewhat flexible interpolating function where A and B are normalization parameters to be determined. The inclusion of the subleading term in power of (1 x) is to account for the contribution from higher Fock state. For each term, the power of (1 x) di↵ers by 2 as suggested by the pQCD analysis []. We refer to the parametrization form Eq. (??) as the ansatz-1. As a phenomenological exploration, the power di↵erence suggested by the pQCD power counting may also be modified. To take into account the model dependence on imposing the (1 x) power di↵erence between g + (x) and g (x), we consider another parametrization by assuming the power di↵erence is 1 as We also consider the following possibility when the large-x distribution of g (x) is suppressed by one power of (1 x) relative to the g + (x) distribution.
(Tianbo, can you provide a better reasoning here? some careful discussion needed for the choice of ansatzes 1 and 2) which we refer to as the ansatz-2. The polarized gluon distribution g(x) decreases to 0 as x ! 0 for both ansatz-1 and ansatz-2. This indicates the helicity correlation between the gluon and its parent nucleon disappears when x ! 0, which is a natural expectation since the relative rapidity becomes infinity.
With the parametrization of the spin-aligned and the spin-antialigned gluon distributions, one can directly obtain the unpolarized and polarized gluon distributions from the sum and the di↵erence of them, The unpolarized gluon distribution xg(x) can be written as a sum of the helicity aligned and anti-aligned glouon distributions To determine the unknown parameters in g +, (x) in Eqs. (??) and (??), we fit the unpolarized gluon distribution from the NNPDF global analysis []. Our precedure described here can be applied to any other gluon distribution given by global analysis or model calculation. To fit h xi (3) g h xi (5) g h xi (7) g h xi (9) g The helicity-aligned gluon distribution is expected to fall of at least (1 x) faster than the valence quark distribution in the nucleon [? ] as x ! 1. The helicityaligned gluon distribution g + (x) is thus expected to have leading (1 x) 4 fall-o↵ as x ! 1 compared to (1 x) 3 fall-o↵ of the valence quark distribution. The helicity anti-aligned gluon distribution g (x) is expected to be suppressed by an additional factor of (1 x) 2 relative to g + (x) as x ! 1. We note that non-leading helicity flip amplitudes and contributions from higher Fock components are power-law suppressed and we only include the next higher Fock component contribution as x ! 1 in our calculation. Besides, we also allow variation in the exponents of (1 x) fall-o↵ through the parameter and instead of fixing the pomeron intercept to be exactly equal to 1, we also allow variation in the exponent ↵. Instead of strictly imposing the power counting at x ! 1 and the pomeron intercept at x ! 0, we only take them as a guidance and phenomenologically introduce two parameters ↵ and to allow the variation of the power behavior at the small and large x regions as usually adopted in global analyses. Since it is expected that the shape of PDFs is not just completely governed by the small and large x limiting behaviors and being motivated by the parametrizations in the global fits of PDFs, we also allow a somewhat flexible interpolating function 2 as suggested by the pQCD analysis []. We refer to parametrization form Eq. (??) as the ansatz-1.
As a phenomenological exploration, the power d ence suggested by the pQCD power counting may be modified. To take into account the model de dence on imposing the (1 x) power di↵erence betw g + (x) and g (x), we consider another parametriza by assuming the power di↵erence is 1 as We also con the following possibility when the large-x distributio g (x) is suppressed by one power of (1 x) relativ the g + (x) distribution.
(Tianbo, can you provide a better reasoning here? s careful discussion needed for the choice of ansatzes 1 2) which we refer to as the ansatz-2. The polarized g distribution g(x) decreases to 0 as x ! 0 for ansatz-1 and ansatz-2. This indicates the helicity relation between the gluon and its parent nucleon d pears when x ! 0, which is a natural expectation the relative rapidity becomes infinity.
With the parametrization of the spin-aligned and spin-antialigned gluon distributions, one can directly tain the unpolarized and polarized gluon distribu from the sum and the di↵erence of them, The unpola gluon distribution xg(x) can be written as a sum o helicity aligned and anti-aligned glouon distribution To determine the unknown parameters in g +, (x Eqs. (??) and (??), we fit the unpolarized gluon dist tion from the NNPDF global analysis []. Our prece described here can be applied to any other gluon dist tion given by global analysis or model calculation. T The helicity-aligned gluon distribution is expected to fall of at least (1 x) faster than the valence quark distribution in the nucleon [? ] as x ! 1. The helicityaligned gluon distribution g + (x) is thus expected to have leading (1 x) 4 fall-o↵ as x ! 1 compared to (1 x) 3 fall-o↵ of the valence quark distribution. The helicity anti-aligned gluon distribution g (x) is expected to be suppressed by an additional factor of (1 x) 2 relative to g + (x) as x ! 1. We note that non-leading helicity flip amplitudes and contributions from higher Fock components are power-law suppressed and we only include the next higher Fock component contribution as x ! 1 in our calculation. Besides, we also allow variation in the exponents of (1 x) fall-o↵ through the parameter and instead of fixing the pomeron intercept to be exactly equal to 1, we also allow variation in the exponent ↵. Instead of strictly imposing the power counting at x ! 1 and the pomeron intercept at x ! 0, we only take them as a guidance and phenomenologically introduce two parameters ↵ and to allow the variation of the power behavior at the small and large x regions as usually adopted in global analyses. Since it is expected that the shape of PDFs is not just completely governed by the small and large x limiting behaviors and being motivated by the parametrizations in the global fits of PDFs, we also allow a somewhat flexible interpolating function where A and B are normalization parameters to be termined. The inclusion of the subleading term in pow of (1 x) is to account for the contribution from hig Fock state. For each term, the power of (1 x) di↵ers 2 as suggested by the pQCD analysis []. We refer to parametrization form Eq. (??) as the ansatz-1.
As a phenomenological exploration, the power di↵ ence suggested by the pQCD power counting may a be modified. To take into account the model dep dence on imposing the (1 x) power di↵erence betwe g + (x) and g (x), we consider another parametrizat by assuming the power di↵erence is 1 as We also consi the following possibility when the large-x distribution g (x) is suppressed by one power of (1 x) relative the g + (x) distribution.
(Tianbo, can you provide a better reasoning here? so careful discussion needed for the choice of ansatzes 1 a 2) which we refer to as the ansatz-2. The polarized glu distribution g(x) decreases to 0 as x ! 0 for b ansatz-1 and ansatz-2. This indicates the helicity c relation between the gluon and its parent nucleon dis pears when x ! 0, which is a natural expectation si the relative rapidity becomes infinity.
With the parametrization of the spin-aligned and spin-antialigned gluon distributions, one can directly tain the unpolarized and polarized gluon distributio from the sum and the di↵erence of them, The unpolari gluon distribution xg(x) can be written as a sum of helicity aligned and anti-aligned glouon distributions To determine the unknown parameters in g +, (x) Eqs. (??) and (??), we fit the unpolarized gluon distri tion from the NNPDF global analysis []. Our preced described here can be applied to any other gluon distri tion given by global analysis or model calculation. To ω = 0 has infinite radius of convergence and therefore even the asymptotic region can be continuously reached from the origin of M(ω) as was also demonstrated in [71]. Therefore along with the advantage that one can perform a reliable extraction of PDFs from the ITDs given the precision, accuracy and availability of ITD over a moderate range of ω < ∼ 15, this also serves as a powerful formalism to calculate higher moments which are not feasible in the LQCD calculations due to power divergent mixing of lower mass dimension operators in the local moments calculations.

VI. ASYMPTOTIC LIMIT OF GLUON IOFFE-TIME DISTRIBUTION
In this section, we calculate analytic expressions for the asymptotic limits of the unpolarized and polarized gluon ITDs associated with the functional forms of xg(x) and x∆g(x) used in Section II. We examine the ω region around which one can approach the asymptotic region of the ITDs. We again emphasize that the purpose of this calculation is not to calculate the gluon PDFs at extremely small x and we therefore consider only the lowx region where the Regge phenomenology is valid and avoid the complication with unresolved issues in the extremely small x < 10 −3 region [36][37][38][39][40][41][42][43][44]. To calculate the asymptotic limits of M(ω) and ∆M(ω), we start with the simplest functional form of PDFs x a (1 − x) b since xg(x) and x∆g(x) determined using functional forms in Eqs. (1) and (2) can be viewed as linear combination of this form with different values of exponents and appropriate normalizations. We first consider the asymptotic expansions (in the limit ω → ∞) of the following integral: (12) where E R (a, b; ω) is the standard aysmptotic expansion of the confluent hypergeometric function, M (a + 1, a + b + 2; iω), and up to order ω −a−R : We also have and We define, Using the above expressions, we can summarize the asymptotic expansions of the unpolarized and polarized gluon ITDs for ansatz-1: +γ C R (α + 1/2, 4 + β; ω) + δ C R (α + 1, 4 + β; ω) and M(ω, µ 2 ) = A S R (α, 4 + β; ω) +γ S R (α + 1/2, 4 + β; ω) + δ S R (α + 1, 4 + β; ω) Similar expressions can be written for ansatz-2 and ansatz-3. We show the asymptotic limits of the unpolarized and polarized gluon ITDs for ansatz-2 in FIGs. 8 and 9, respectively. For the demonstration purpose, we select one arbitrary set of parameters from to the fit to 1 replica of the NNPDF unpolarized gluon distribution described in Section II. The M(ω) and ∆M(ω) approach the asymptotic limits around ω ∼ 15 as can be seen in FIGs. 8 and 9. It is important to note that if future LQCD calculations of gluon ITD can reach the region ω ∼ 15, they will be able to provide nonperturbative information to the Lipatov's pomeron [36,37].  Using the fact that gluon PDF diverges much faster than the valence quark PDF in the limit x → 0, one can show that the asymptotic limit of the ITD corresponding to nucleon valence quark distribution will set in at earlier ω compared to the gluon ITDs, also noted in [47]. This implies that the asymptotic region of the nucleon valence quark ITD can be approached easily in the nonperturbative calculations compared to the gluon ITDs.

VII. APPLICATIONS TO LATTICE QCD CALCULATIONS OF PDFS
In recent years, several LQCD methods have been proposed and developed to probe the light-cone structure of hadrons, including the path-integral formulation of the deep-inelastic scattering hadronic tensor [72], coordinatespace method for the calculating light-cone distribution amplitudes [59], inversion method [73], quasi-PDFs/LaMET [17,18], pseudo-PDFs [60], and good lattice cross sections [61,74]. For the most recent review of LQCD calculations of PDFs, see [75].
The extraction of PDFs from LQCD calculations has received great interest since Ji's proposal in [17,18]. Instead of directly calculating the light-cone correlation functions which define the PDFs, one can extract them from the spatial correlation of parton fields calculable on the Euclidean lattice. We begin this section by acknowledging that any LQCD calculations of PDFs using any of the above formalisms share the common challenge of how best to extract a continuous distribution from discrete data, compounded by a limited number of data points due to a finite range of spatial separations and hadron momenta. As the ITDs in question herein are Fourier transforms of the underlying PDFs, the available data in a discrete and limited domain leads to an ill-posed inverse problem well-known to the LQCD community [69,[76][77][78]. Similar to functional forms used in global fits, different PDF parametrizations obviate the ill-posed inverse Fourier transform at the cost of additional systematic errors in the determination of PDFs. The analyses of the global fitting community have matured over the past several decades to where these systematics can be reliably estimated. In the ideal scenario, precise LQCD data will allow this systematic error to likewise be estimated and corrected by fitting several models and examining relevant figures of merit. Significant investigations are being performed at present to handle this problem better; even the neural network approach for determining PDFs from LQCD data has become feasible now, as we will discuss below.
With the current resources available in LQCD calculations of PDFs, it remains a challenge to obtain precise and accurate data at large ω = z ·p, where p is the hadron momentum and z is the spatial separation between parton fields [17,60] or gauge-invariant currents [59,61]. This problem was actually anticipated in [47]. As a remedy to this problem, it was proposed to consider the low ω region, where LQCD data can be precise and accurate, and the large ω domain of the ITD separately. The LQCD accessible small ω-region and the large ω region could then be matched with the knowledge of Regge phenomenology as described in [47,63,71]. This approach therefore enables one to obtain reliable estimates of PDFs in the entire x-region. Recently, a similar proposal was suggested in [79] to circumvent such a problem in the lattice QCD calculation of quasi-PDFs using the LaMET formalism.
First-principles LQCD determinations of the gluon unpolarized and polarized distributions, with controlled statistical and systematic uncertainties, have been of particular interest recently with significant theoretical developments [80][81][82][83][84] soliciting a synergy of increasing importance between experimental and theoretical efforts. However, convincing LQCD calculation of gluon PDFs has remained very challenging. It is a big challenge to obtain precise ITD at large ω [70] and to reach the asymptotic region of the ITD, especially for the gluonic observables. In this light, the phenomenological knowledge of the asymptotic limits of the gluon ITDs can be considered as an interesting opportunity to match the low-ω ITD from LQCD calculations and large-ω asymptotic limits utilizing the method discussed in Section VI.
LQCD data with larger separations between the quark/gluon fields are seen to be consistent with zero, even change sign, getting affected by various systematic errors due to issues in the renormalization and perturbative matching, and are also expected to be contaminated by higher-twist effects [79]. This can be directly seen from the nucleon ITD calculated in the pseudo-PDF formalism and also in the renormalized matrix elements of the quasi-PDF/LaMET formalism both of which start with the same matrix elements [17]. This leads to the potential issue with the LQCD calculations of PDFs; the ITD falls off much faster compared to the ITD constructed from the phenomenologically determined PDFs, e.g. nucleon valence quark distribution which is very well constrained from different global fits. From the phenomenological point of view, this faster-fall of the LQCD data immediately results in a faster-converging PDF at low x and softer fall-off at large x. This has been demonstrated using neural network framework applied to the extraction of the nucleon valence quark distribution (which is assumed to be a simpler LQCD calculation compared to the gluon PDF, sea-quark distributions, etc.) from LQCD calculated matrix elements using the quasi-PDF [85] and pseudo-PDF [86] approaches. A similar observation can be found in a recent Monte Carlo based analysis of LQCD data for nucleon valence quark PDF [87]. To investigate these issues, for example, one can construct ITDs from the well-determined nucleon valence PDFs of global fits and compare them with the LQCD calculated matrix elements and examine at which ω the LQCD determined ITD starts to deviate from the phenomenological ITD. With the help of the asymptotic limit of the ITD, one can also investigate the possible sources of discrepancy between the LQCD calculation and the ITD derived from global fits, such as higher-twist contributions. Investigation based on these observations is an ongoing subject of a future research project.
Along with the calculation of PDFs as discussed above, an ITD can also be used to determine moments of PDFs in a reliable manner. One can immediately fit the LQCD calculated ITD using Eqs. (10) and (11) to extract moments of PDFs. The accuracy and the number of accessible moments will depend on the available ω-range and on the precision of the ITD. This can provide an alternative approach for extracting local moments of the distribution [88,89] by directly fitting the ITD. In particular, this can be useful for the gluon [70] and sea-quark [90] distributions for which the LQCD data is seen to be much noisier compared to the non-singlet quark distributions and the ill-posed inverse problem of extracting PDFs is much severe.

VIII. CONCLUSION AND OUTLOOK
In this paper, we investigate the unpolarized and polarized gluon distributions and their applications to the Ioffe-time distributions, which are closely related to the extraction of PDFs from lattice QCD calculations. We parametrize the gluon distributions in the helicity basis and construct the functional form motivated by the counting rules based on perturbative QCD analyses at large x and the phenomenological behavior at low x. Once the helicity-basis gluon distributions are determined, one can easily obtain the unpolarized and polarized gluon distributions. By fixing the (1 − x) difference and apply the same polynomial factor, we determine the helicity-basis distributions with the unpolarized gluon distribution and infer the polarized gluon distribution using two different ansätze. Although the results are both in relatively good agreement with the global analysis, the two results show a sizable difference from each other. To see the model uncertainty, we also perform a simultaneous fit to unpolarized and polarized gluon distributions from global analyses using a more flexible parametrization form.
As an application, we calculate the Ioffe-time distributions and discuss the possibility that the asymptotic expansion of the Ioffe-time distributions in the large-ω region which can be combined with the future lattice QCD calculations that are limited in a relatively smaller ωrange and can guide to extrapolate the lattice data in the larger ω-region. Using our calculation, we have demonstrated that the magnitude of the polarized gluon Ioffetime distribution within a moderate Ioffe-time window in a future nonperturbative QCD calculation can provide important constraints on the gluon spin content in the nucleon. We have discussed, as a general application of this method, the discrepancy between the fall-off rate of the phenomenologically well-determined Ioffe-time distributions and those calculated in lattice QCD calculations can serve as an interesting platform to investigate higher twist effects in lattice QCD calculations. Following recent observations, we have discussed that with the present and possibly near future resources and numerical techniques, lattice QCD calculations alone might not be sufficient enough to extract the full x-dependence of PDFs with desired precision and accuracy. In such cases, phenomenologically well-determined Ioffe-time distributions in the large ω-region can provide complementary information to the ongoing efforts of the calculation of xdependent hadron structures within first-principles nonperturbative QCD calculations.

RSS thanks Colin Egerer, Xiangdong Ji, Joseph
Karpie, Nikhil Karthik, Kostas Orginos, and David Richards for useful suggestions and discussions. We thank Jian-Wei Qiu who provided insights and expertise that greatly assisted this research. RSS is supported in part by U.S. DOE grant No. DE-FG02-04ER41302. TL is supported in part by National Natural Science Foundation of China under Contract No. 11775118. This work is also supported by the U.S. Department of Energy contract No. DE-AC05-06OR23177, under which Jefferson Science Associates, LLC, manages and operates Jefferson Lab, and within the framework of the TMD Topical Collaboration.