Unpolarized gluon distribution in the nucleon from lattice quantum chromodynamics

In this study, we present a determination of the unpolarized gluon Ioffe-time distribution in the nucleon from a first principles lattice quantum chromodynamics calculation. We carry out the lattice calculation on a $32^3\times 64$ ensemble with a pion mass of $358$ MeV and lattice spacing of $0.094$ fm. We construct the nucleon interpolating fields using the distillation technique, flow the gauge fields using the gradient flow, and solve the summed generalized eigenvalue problem to determine the glounic matrix elements. Combining these techniques allows us to provide a statistically well-controlled Ioffe-time distribution and unpolarized gluon PDF. We obtain the flow time independent reduced Ioffe-time pseudo-distribution, and calculate the light-cone Ioffe-time distribution and unpolarized gluon distribution function in the $\overline{\rm MS}$ scheme at $\mu = 2$ GeV, neglecting the mixing of the gluon operator with the quark singlet sector. Finally, we compare our results to phenomenological determinations.

In this study, we present a determination of the unpolarized gluon Ioffe-time distribution in the nucleon from a first principles lattice quantum chromodynamics calculation. We carry out the lattice calculation on a 32 3 × 64 ensemble with a pion mass of 358 MeV and lattice spacing of 0.094 fm. We construct the nucleon interpolating fields using the distillation technique, flow the gauge fields using the gradient flow, and solve the summed generalized eigenvalue problem to determine the gluonic matrix elements. Combining these techniques allows us to provide a statistically well-controlled Ioffetime distribution and unpolarized gluon PDF. We obtain the flow time independent reduced Ioffetime pseudo-distribution, and calculate the light-cone Ioffe-time distribution and unpolarized gluon distribution function in the MS scheme at µ = 2 GeV, neglecting the mixing of the gluon operator with the quark singlet sector. Finally, we compare our results to phenomenological determinations.

I. INTRODUCTION
Gluons, which carry color charge and serve as the mediator bosons of the strong interaction, play a key role in the nucleon's mass and spin. Confinement in quantum chromodynamics (QCD) ensures that no free quarks or gluons have been observed, so analyses of hadrons involving high energy scattering rely on QCD factorization [1]. Factorization separates the perturbatively-calculable hard-scattering quark and gluon dynamics from the non-perturbative collinear dynamics, described by parton distribution functions (PDFs) of the relevant hadrons.
There are long-standing efforts to conduct global analyses [2][3][4][5][6] of data from available deep inelastic scattering (DIS) and related hard scattering processes to explore the nature of the PDFs. It is essential to have a clear and precise understanding of the gluon PDF in order to calculate the cross-section for Higgs boson production [7] and jet production [8] at the Large Hadron Collider (LHC), and J/ψ photo production [9] at Jefferson Lab. Future colliders, such as the Electron Ion Collider (EIC) [10][11][12], which is to be built at Brookhaven National Lab, and the Electron Ion Collider in China (EicC) [13], are expected to make significant impact on the precision of the gluon PDFs. While the precision of the extracted gluon distribution x g(x) has been improved over the last decade, several issues remain unresolved; for example, the suppression in the momentum fraction region 0.1 < x < 0.4 when ATLAS and CMS jet data are included [3] and how to obtain a more precise determination of g(x) are subjects of ongoing efforts.
The determination of PDFs from lattice QCD is of particular theoretical interest to directly explore the nonperturbative sector of QCD from the first principles. To achieve this goal, there have been several proposals for the extraction of the x-dependent hadron structure from lattice QCD calculations, such as the path-integral formulation of the deep-inelastic scattering hadronic tensor [14], the operator product expansion [15], quasi-PDFs [16,17], pseudo-PDFs [18], and lattice cross-sections [19,20]. Lattice QCD is formulated in Euclidean space, so the bilocal light-cone correlators that are necessary to extract the PDFs cannot be evaluated directly, because they require operators containing fields at light-like separations, z 2 = 0, which cannot exist in Euclidean space. The quasi-PDF framework [16] circumvents this drawback by calculating matrix elements associated with equal time and purely space-like field separations with hadron states at non-zero momentum, p z . The corresponding quasi-PDFs can be matched to the light-cone PDFs when the hadron momentum is large, by applying the Large Momentum Effective Theory (LaMET) [17]. These calculation techniques have been explored extensively in numerical lattice calculations. (For recent reviews see [21,22] and the references therein.) There have been significant achievements in lattice QCD calculations of x-dependent hadron structure: the nucleon valence quark distribution using pseudo-PDFs [23], the calculation of the pion valence distribution using the lattice cross section, quasi-PDF and pseudo-PDF frameworks [24][25][26][27][28], the kaon PDF calculation using the quasi-PDF formalism [29], nucleon unpolarized and helicity distributions within quasi-PDF formalism [30][31][32], the unpolarized and

A. Matrix Elements
To access the unpolarized gluon PDF, we calculate the matrix elements of a spin-averaged nucleon for operators composed of two gluon fields connected by a Wilson line, which have the general form Here, z µ is the separation between the gluon-fields, p µ is the 4-momentum of the nucleon, W [z, 0] is the standard straight-line Wilson line in the adjoint representation, for the gauge field A µ , where P indicates that the integral is path-ordered. The matrix elements can be decomposed into invariant amplitudes, M pp , M zz , M zp , M pz , M ppzz and M gg using the four-vectors, p µ and z µ , and the metric tensor g µν [56]. These amplitudes are functions of the invariant interval z 2 and the Ioffe-time p · z ≡ −ν [44]. The light-cone gluon distribution is obtained from where z is taken in the light-cone "minus" direction, z = z − , and p + is the momentum in the light-cone "plus" direction. The PDF is determined by the M pp amplitude, The density of the momentum carried by the gluons, G(x) = x g(x) is the natural quantity in this definition of the gluon PDF, rather than g(x). The field-strength tensor G µα is antisymmetric with respect to its indices and g −− = 0, so the left hand side of Eq. (3) reduces to a summation over the transverse indices i, j = x, y; perpendicular to the direction of separation between the two gluon fields. The matrix element M ti;it decomposes into the invariant amplitudes [56] where M gg is a contamination term. The matrix element cancels the contamination term from M ti;it [56]. Thus, the proper combination of the matrix elements to extract the twist-2 invariant amplitude, M pp is For spatially-separated fields, the gauge link operator has extra ultraviolet divergences not present for light-like separated fields. The combination of matrix elements M ti;it is multiplicatively renormalizable [57]. And, because of the antisymmetry of the gluon fields, the combination M ji;ij can be written as which contains only one set of indices {µα; λβ}, making explicit the fact that this matrix element is multiplicatively renormalizable too [58]. Furthermore, both M ti;it and M ji;ij have the same one-loop UV anomalous dimension [56], making the whole combination in Eq. (7) multiplicatively renormalizable at the one-loop level, at least.

B. Reduced Matrix Elements
Similar to space-like separations, the extended gluon operator has additional link-related ultraviolet (UV) divergences which are multiplicatively renormalizable [59][60][61]. These UV divergences can be cancelled by taking appropriate ratios. We combine the matrix elements from Eq. (7) which we denote by M(ν, z 2 ) for the rest of the paper, and take the ratio [46] of the combination to its rest-frame value, keeping the separation same. This ratio cancels out the ν-independent UV factor Z(z 2 /a 2 ), making the ratio UV-finite. The kinematic factors remaining in the ratio can be removed by taking the ratio of the non-zero separation to the zero separation matrix elements, at fixed Ioffe-time, in both the numerator and denominator [48].
The resulting reduced matrix element, the reduced pseudo-ITD, can be written as: Taking the ratio, we also eliminate z 2 -dependent, but ν-independent, non-perturbative factors that M(ν, z 2 ) may contain. The residual polynomial "higher twist" dependence on z 2 , if visible, should be explicitly fitted in order to separate it from the twist-2 contribution.

C. Position-space Matching
The reduced pseudo-ITD has a logarithmic z 2 dependence. We relate the reduced pseudo-ITD, M(ν, z 2 ), to the gluon and singlet quark light-cone ITDs, I g (ν, µ 2 ) and I S (ν, µ 2 ) in the MS scheme through the short distance factorization relationship with z 2 as the hard scale. Here, I g (ν, µ 2 ) is related to the gluon PDF, g(x, µ 2 ), by The product x g(x, µ 2 ) is an even function of x, so the real part of I g (ν, µ 2 ) is given by the cosine transform of x g(x, µ 2 ), while its imaginary part vanishes. Neglecting the higher twist terms of M zz , M zp , M pz , M ppzz , and keeping just the M pp term, the one-loop matching relation is [56,62], The singlet quark Ioffe-time distribution I S (ν, µ 2 ) is related to the singlet quark distribution, summed over quark flavors. The Altarelli-Parisi kernel, B gg (u), is given by and the quark-gluon mixing kernel is given by where the plus-prescription is andū ≡ (1 − u). Here, γ E is the Euler-Mascheroni constant and C F is the quadratic Casimir operator in the fundamental representation. Determining the singlet quark Ioffe-time distribution requires evaluation of the disconnected diagrams, which involves the computationally demanding calculation of the trace of the all-to-all quark propagator [63], but contribute only a little to the matching. We neglect quark-gluon mixing in this calculation and implement the matching relation

A. Gluonic Current Calculation
The gluonic currents, inserted into the nucleon to calculate the matrix elements, are not connected to the nucleon state by any quark propagator, so the currents are largely decoupled from the nucleon part of the calculation itself. As a result, on the lattice, we can calculate the gluonic currents and the nucleon two-point correlators separately and combine them together to obtain the three-point correlators from which we extract the matrix elements. On the lattice, the gluonic current can be written with the Wilson line in the fundamental representation as The field-strength tensor can be expressed in terms of the (1 × 1) plaquette operator, U (1×1) µν , as [64] −i 2 where a is the lattice spacing and β = 6/g 2 s . One-third of the trace is subtracted here to enforce the traceless property of the Gell-Mann matrices. The (1 × 1) plaquette operator is defined as the product of the link variables forming a (1 × 1) loop on the lattice, To reduce statistical fluctuations, we take the average of the four possible plaquette operators that can be constructed by changing the signs of µ and ν. Finally, we combine the gluonic currents O(G ti , G it , z) and O(G ji , G ij , z) to calculate M pp . Accounting for the sign change of the gluonic current with the "temporal" index in Euclidean spacetime, the total gluonic current becomes

B. Gradient Flow
In our calculation, we apply the gradient flow [65][66][67] to reduce ultraviolet fluctuations and improve the signal-tonoise ratio for the gluon observables. To implement this technique, the flowed gauge field, B µ (τ, x), is defined by following the procedure in [65],Ḃ where the flowed gauge field is subjected to the boundary condition B µ (τ = 0, x) = A µ (x). Here τ is the flow time and we abbreviate differentiation with respect to τ by a dot. The flow equation of the gauge field is a diffusion equation and the evolution operator in the momentum space acts as an UV regulator for τ > 0. As a result, the gradient flow exponentially suppresses the UV field-fluctuations, which corresponds to smearing out the original degrees of freedom in coordinate space. The operators constructed using flowed gauge fields with positive flow time enter into the relevant theories at length scales of ∼ √ 8τ . On the lattice, the gradient flow is implemented by defining the flowed link variable, V µ (τ, x) as [65]: where g 0 is the bare coupling, S(V µ (τ, x)) is the flowed action, V µ (τ = 0, x) has the boundary condition of being equal to the link variable, U µ (x), and ∂ x,µ stands for the natural SU(3)-valued differential operator with respect to V µ (τ, x). The action, S(V µ (τ, x)) is a monotonically decreasing function of τ , and the gradient flow corresponds to a continuous stout-link smearing procedure [68]. We use unimproved Wilson flow and calculate the gluonic currents with flow times from τ = a 2 to τ = 3.8a 2 . We construct the double ratio of Eq. (9) using the flowed matrix elements, which further reduces UV fluctuations and suppresses the flow time dependence. The residual τ -dependence is removed by fitting the flowed reduced matrix elements to an appropriate functional form which, in turn, gives us the reduced pseudo-ITD at zero flow time.

C. Nucleon Two-point Correlator
We calculate the nucleon two-point correlators by applying interpolators at the source time-slice and the sink time-slice on the lattice. We apply distillation [69], a low-rank approximation to the gauge-covariant Jacobi-smearing kernel, J σ,nσ (t) = 1 + σ∇ 2 (t) nσ nσ [70]. The tunable parameters {σ, n σ } ensure that, in the large iteration limit, the kernel approaches that of a spherically-symmetric Gaussian. The quark fields are smeared using the distillation smearing kernel where where N c is the dimension of the color space, N x , N y , N z are the extents of the lattice in the three spatial directions, and N D is the dimension of the distillation space. The k th column of V D (t), ν (k) x (t) is the k th eigenvector of the second-order three-dimensional differential operator, ∇ 2 , evaluated on the background of the spatial gauge fields of time-slice t, once the eigenvectors have been sorted by the ascending order of the eigenvalues. The two-point correlator for the nucleon can be written as where, and Here, Φ i (t) and P(t m , t n ) are referred to as elementals and perambulators, respectively; D is the lattice representation of the Dirac operator; α,ᾱ, β,β, γ,γ are the spin indices; a, b, c are the color indices. The Φ i (t) encodes the structure of the interpolating operator as well as has a well-defined momentum, while P(t m , t n ) encodes the propagation of the quarks, and does not have have any explicit momentum projection. Elementals can be decomposed into terms that act only within coordinate and color space, like Γ, and only within spin space, like S αβγ . We adopt distillation for two reasons. First, the computationally demanding parallel transporters of the theory, the perambulators, depend only on the gauge field, and not on the interpolators. Therefore, we can calculate the perambulators on an ensemble of gauge fields once, and then reuse them for an extended basis of interpolators, thus reducing the computational cost to a great extent. This extended basis of interpolators is the key to perform a successful summed generalized eigenvalue problem (sGEVP) analysis [71], enabling us to attain a clear signal for the ground state nucleon.
Second, distillation admits a momentum projection both at the source interpolating operator, and at the sink interpolating operator, in contrast to the more usually adopted methods. Thus for the gluonic three-point functions computed here, we are able to impose momentum projection at all three time-slices, ensuring the most complete possible sampling of the lattice. Moreover, the low-lying spectra of the nucleon can be faithfully captured with a relatively small number of distillation eigenvectors [72], thus lowering the cost of the calculation further. The expectation is that N D should scale as the physical volume, and the cost of computing the corresponding correlation functions scales as N 4 D for the case of the nucleon. In this calculation, we employed N D = 64 eigenvectors. The efficacy of distillation for the calculation of nucleon charges was demonstrated in ref. [73], and subsequently extended to the case of the nucleon in motion [74]. Recently, the unpolarized, isovector PDF of the nucleon has been computed using the same ensemble within the distillation framework [55].

D. Interpolators
The lattice regulator explicitly breaks the continuum SO(3) rotational symmetry, so the associated symmetry group reduces to the double-cover octahedral group, O D h for the nucleon at rest. Although there are six irreducible representations (irreps.) available in O D h , we focus on G 1g , because the states with continuum spin 1 2 , such as the ground state nucleon, are subduced onto this irrep. Here, the subscript g stands for positive parity. At non-zero spatial momenta, the O D h group breaks into further little groups depending on the direction of the boost. We consider boosts only along the z-direction, so the associated little group is the order-16 dicyclic group or Dic 4 .
To calculate the low-lying spectra of the nucleon, we include interpolators with zero orbital angular momentum, which have the largest overlaps with the ground state of the nucleon. For the lowest excited-states, we include Spatial momentum Interpolators interpolators with gauge-covariant derivatives acting on the quark fields to capture the effect of the non-zero angular momenta between the quarks [75]. All these interpolators are "non-relativistic", in the sense that they feature only the upper components of the Dirac spinors. We also include the interpolators that have derivatives of second order and form combinations corresponding to the commutation of two gauge-covariant derivatives acting on the same quark field. These interpolators, also referred to as hybrid interpolators [76], vanish in the absence of a gauge-field and correspond to the chromomagnetic components of the gluonic field-strength tensor. We tabulate our choice of interpolators for the nucleon at rest as the first row in Table I, using the spectroscopic notation of: X 2S+1 L π J P where X is the nucleon, N ; S is the Dirac spin; L = S, P, D, . . . is the orbital angular momentum; π = S, M or A is the permutation symmetry of the derivatives; J is the total angular momentum; and P is the parity. For the construction of the three-point correlators needed for the unpolarized distributions, we take the sum of the spin = + 1 2 and spin = -1 2 nucleon two-point correlators.
For the case of the correlation functions at non-zero spatial momentum, parity is no longer a good quantum number and further operators are classified according to their helicity. We therefore include operators corresponding both to higher spin, and to negative parity, in our basis within the little group Dic 4 . We choose the direction of momenta to be in the same direction of the polarization to ensure longitudinal polarization. We access the unpolarized gluon PDF by taking the sum of helicity = + 1 2 and helicity = -1 2 nucleon two-point correlators. The basis of interpolators is tabulated as the second row in Table I.

E. Momentum Smearing
To access a wide range of Ioffe-times, we perform the lattice calculation at multiple spatial momenta. On the lattice, the spatial momentum is discretized and expressed as Here, L = 32, is the spatial extent of the lattice. For p, where l > 3, we enhance the overlap of the interpolators onto the lowest-lying states in motion by applying momentum smearing [77]. We follow the procedure introduced in [74] and add a phase to the distillation eigenvectors for higher momenta to preserve translational invariance, which is essential for the projection onto the states of definite momenta. The "phased" distillation eigenvector becomes, It is sufficient to modify the previously computed eigenvectors to perform calculation at the higher lattice momenta, though the perambulators and the elementals need to be recalculated with these "phased" eigenvectors. For our calculation, choosing gives the momentum smearing needed for boosts up to p = 6 × 2π aL .

IV. LATTICE DETAILS
We perform our calculation on an isotropic ensemble with (2 + 1) dynamical flavors of clover Wilson fermions with stout-link smearing [68] of the gauge fields and a tree-level tadpole-improved Symanzik gauge action, with approximate lattice spacing, a ∼ 0.094 fm and pion mass, M π ∼ 358 MeV, generated by the JLab/W&M collaboration [78]. The rational hybrid Monte Carlo (RHMC) algorithm [79] is used to carry out the updates. One iteration of four-dimensional stout-smearing with the weight ρ = 0.125 for the staples is used in the fermion action. After stout-smearing, the tadpole-improved tree-level clover coefficient, C SW , is very close to the non-perturbative value. This is confirmed using the Schrödinger functional method for determining the clover coefficient non-perturbatively [78]. The tuning of the strange quark mass is done by first setting the quantity, 1678. This quantity is independent of the light quark masses to the lowest order in χPT, depending only on the strange quark mass [80]. So, it can be tuned in the SU(3) symmetric limit. The resulting value of the strange quark mass is then kept fixed as the light quark masses are decreased in the (2+1) flavor theory to their physical values.
We use 64 temporal sources over 349 gauge configurations, with each configuration separated by 10 HMC trajectories. The two light quark flavors, u and d are taken to be degenerate and the lattice spacing was determined using the w 0 scale [81]. We summarize the parameters of the ensemble in Table II.

V. VARIATIONAL ANALYSIS
To check whether the two-point correlators give us the expected results, we investigate the associated principal correlators and extract the energy spectra by performing a variational analysis for the nucleon at rest in the G 1g channel and for all the boosted frames in the Dic 4 little group with the interpolators in Table I. This fitting procedure is discussed in detail in [72,73,75]. We only summarize the procedure here. We solve the GEVP of Eq. (A9) over a range of t 0 . We then define optimal interpolators, in the variational sense, for the energy eigenstates, |n through i u i nŌN,i . Here,Ō N,i are the interpolators used in the calculation and u i n are the weights of these interpolators that define the optimal interpolator. The energy associated with each state |n is obtained by fitting its principal correlator according to In our fitting procedure, we aim to ensure that the principal correlators are dominated by the leading exponential. Thus in each of our fits, we choose t 0 such that we obtain an acceptable χ 2 /d.o.f., that the value of A n is small, typically less than 0.1, and that, for each principal correlator, λ n (t, t 0 ), the subleading energy E n is larger than than the leading energies for all the principal correlators. This indicates that the matrix of two-point correlators is to a large degree, saturated by the lowest-lying states.
In Fig. 2 and 3, we show fits to the leading principal correlators for the nucleon subduced onto the little group, Dic 4 for p = 2 × 2π aL = 0.82 GeV, and p = 6 × 2π aL = 2.46 GeV, respectively. For each panel, the blue band is the reconstruction from the fitted parameters. The approach of the plateaux close to unity at large times is indicative of the small value of A n in the fits, and the small contribution of the other states to each principal correlator.
In Fig. 4, we plot the ground state nucleon energies extracted using the variational analysis with respect to the spatial momentum, together with the expectations from the continuum dispersion relation. Fig. 4 shows that for lower momenta, the unphased ground state nucleon energies are in excellent agreement with the continuum dispersion relation. At p = 3 × 2π aL = 1.23 GeV, the ground state energy starts to deviate, but from p = 4 × 2π aL = 1.64 GeV, after phasing, the ground state energy starts to align with the continuum dispersion curve, indicating that adding a phase to the distillation eigenvectors with ζ = 2 2π L resulted in a significant increase in the overlap of the interpolators onto the lowest-lying states in motion.  46 GeV, subduced onto the little group, Dic4, on the ensemble a094m358, for t0 = 6. The plots show λn(t, t0) e En(t−t 0 ) data on the y-axes and the lattice time-slices on the x-axes; the blue bands are the two-exponential fits as described in the text. The top, middle and bottom panels show the principal correlators for the ground state, the first excited-state and the second excited-state respectively. In each panel, the energy corresponding to the leading exponential state is labelled by En.

A. Three-point Correlator
We calculate the matrix elements by first computing the three-point correlators by inserting gluonic currents between the source and the sink of the two-point correlators. The three-point correlator can be expressed as where C 2pt (t) is the nucleon two-point correlator with source-sink separation t in lattice units and t g is the time-slice on which the gluonic current is inserted.

B. sGEVP Method
We implement the sGEVP method [71,82] to extract the matrix elements from the three-point correlators, a combination of the summation method [83] and GEVP [75] method which begins with the formation of the summed three-point correlation functions formed from our basis of interpolating operators We provide details of the method in appendix A, but the salient feature is that for sGEVP, the systematic error decays as t exp(−∆E t) , which is much faster than the exp(−∆E t) decay for GEVP [75]. This allows us to access the matrix elements at a much smaller temporal separation than would be possible with GEVP. This is crucial for hadron structure calculations, since the signals tend to be heavily contaminated by noise as the temporal separation is increased. sGEVP utilizes the lowest-lying spectra, conveniently calculated using distillation, by rotating the threepoint correlator matrix by a suitable angle, removing much of the excited-state contaminations, and therefore performs better than the summation method [83], which involves only the ground-state nucleon.
In principle, increasing the number of states, N , in the sGEVP analysis should lead to a larger ∆E, which enables matrix elements to be extracted from even smaller temporal separations. This, however, also increases the computational cost, because the N × N correlator matrix needs to be constructed, and makes solving the GEVP for the nucleon two-point correlator matrix more challenging.

C. Bare Matrix Elements
Our calculation requires the extraction of the matrix elements at multiple flow times, multiple nucleon momenta and multiple separations between the gluon fields. We perform the calculation for flow times τ /a 2 = 1.0, 1.4, 1.8, 2.2, 2.6, 3.0, 3.4 and 3.8. For each flow time, we calculate the matrix elements for nucleon momenta, p = 2πl aL where l = 0 to 6, and for field separations, z = s a where s = 0 to 6; a being the lattice spacing. We construct the effective matrix element, M eff (t, z, p, τ ) for each flow time, nucleon momentum and field separation, using the formulation described in appendix A and fit the matrix elements using the functional form in Eq. (A12), which can be written in simplified notation and arguments as Here, A is the matrix element we wish to extract. To perform the fit of Eq. (33) for a particular nucleon momentum, p, we first fit the matrix element for z = 0 using a Bayesian analysis and determine the corresponding fitted value of the parameter, ∆E. As the hadronic spectrum is determined by the two-point correlators, we use the value of ∆E obtained from the fit to the matrix element for z = 0 as the prior for our subsequent fits to the matrix elements for z > 0 at that particular nucleon momentum. We set the prior-width of ∆E for z > 0 to be three times larger than the uncertainty in ∆E and allow for random priors in XMBF [84]. The priors are chosen randomly from normal distributions with the given prior-widths. We perform a simultaneous and correlated fit to the matrix elements for z = {1, 2, 3, 4, 5, 6} × a = 0.094 fm, 0.188 fm, 0.282 fm, 0.376 fm, 0.470 fm, 0.564 fm respectively, where i = 1, 2, · · · 6 and the ∆E is assumed to be the same for matrix elements at a fixed nucleon momentum and flow time. This procedure is particularly helpful for a well-controlled fit to the large momentum matrix elements for which the signal-to-noise ratio is poor, especially at flow times τ /a 2 < 1.6. In Fig. 5, we illustrate our fits to the matrix elements for τ /a 2 = 1.0, in the upper row and for τ /a 2 = 3.0 in the bottom row. Here, we compare the fitted matrix elements among the momenta, p = {1, 6} × 2π aL = 0.41 GeV, 2.46 GeV respectively; and the separations, z = {0, 1, 6} × a = 0, 0.094 fm, 0.564 fm respectively, and list the fitted parameters in Table III. One can immediately see that the ∆E values determined for the non-zero separations are almost identical compared to that obtained for the matrix elements at z = 0 where no prior is assigned on the fit parameter ∆E. This, along with the goodness of the fit in the extraction of the matrix elements for the non-zero separations, indicates the validity of our fitting procedure.
From Fig. 5 and the corresponding fit parameters in Table III we see that the lattice data are described well by our fit procedure. The χ 2 /d.o.f. shows that the choice of prior-width for ∆E at z > 0 is an appropriate one. We notice from Fig. 5 that the matrix elements for z = 6a = 0.564 fm, have a flat behavior as a function of the source-sink separations. This can also be understood from the smallness of B-parameters listed in Table III, with relatively larger uncertainties. The nucleon two-point correlators have quite good signal-to-noise ratios up to the source-sink separation t = 9a = 0.846 fm at p = 6 × 2π aL = 2.46 GeV, as can be seen from Fig. 3. Fig. 5 shows, however, that the matrix elements almost lose any statistical signal around source-sink separation t = 6a = 0.564 fm, which is expected as the nucleon momentum increases. As shown in [85], the optimized interpolators reduce the excited-state contributions allowing us to start the fit at significantly earlier source-sink separations. In support of this, we indeed see from Fig. 5 that the matrix elements for p = 1 × 2π aL = 0.41 GeV reach a plateau around the source-sink separation, t = 4a = 0.376 fm. We note that lattice QCD calculations of the gluonic observables are, in general, much noisier than quark matrix elements. Measures of the goodness of the fits do not necessarily reflect all the systematic uncertainties in our extractions of the fit parameters A, B, and ∆E. However, by using N interpolators within a variational approach, we are better able to sample the Hilbert space in a particular irrep. in finite volume. This has been proven successful in nucleon structure calculation in [73]. The crucial insight is that projecting to the definite finite volume states via the variational solutions allows us to take advantage of the orthogonality of the states in the Hilbert space [82]. There are clearly residual excited-states present as constructing the ideal basis is unrealistic. However, a significant improvement is achieved by incorporating a moderate number of interpolators and applying distillation, one of the most computationally cost-effective methods for implementing a large number of interpolators. Therefore, by adding multiple interpolators we have attempted to systematically improve the determination of A, B, and ∆E in this calculation. Further investigation with larger statistics will be necessary for complete estimate of all the systematic uncertainties associated with excited-state contamination at large nucleon momenta.

D. Reduced Matrix Elements and Zero Flow time Extrapolation
From the bare matrix elements, we calculate the reduced matrix elements using the double ratio in Eq. (9) for different flow times, nucleon momenta and field separations. We present the reduced matrix elements for four different values of τ /a 2 in Fig. 6. We expect the higher twist contributions, discretization effects, and flow time dependence to be minimized through this double ratio.
From the reduced matrix elements at different flow times, we calculate the reduced pseudo-ITD distribution by extrapolating to zero flow time. At fixed values of the field separation, z, and nucleon momentum, p, we find that the τ -dependence is best fit by a linear form, M(τ ) = c 0 + c 1 τ , which we use to determine the reduced pseudo-ITD matrix elements for the subsequent analyses. The values of the fitted parameters are tabulated in appendix B. Out of 36 different fits, we present six examples of such extrapolation in Fig. 7 and for all extrapolations, we find χ 2 /d.o.f. < 1.0. Finally, we present the reduced pseudo-ITD in the zero flow time limit in Fig. 8.

VII. DETERMINATION OF GLUON PDF AND COMPARISON WITH PHENOMENOLOGICAL DISTRIBUTION
Determining PDFs from lattice calculations involves the challenge of how best to extract a continuous distribution from the discrete lattice data, compounded by a limited number of data points due to a finite range of field separations and hadron momenta, and therefore a finite range of ν. By performing a phenomenological analysis of the NNPDF unpolarized gluon PDF [4], it has been found in [86] that a ν-range that is much larger than the present calculation, or any available lattice QCD determination of the gluon ITD [40,41], is necessary to determine the gluon distribution in the entire x-region from the ITD data. Therefore, we do not expect a proper determination of the gluon distribution in the entire x-region, especially in the small-x domain. However, given our lattice data in a limited region, namely ν ∈ [0, 7.07], we extract the gluon PDF from the reduced pseudo-ITD using the Jacobi polynomial parameterization proposed in [23]. The details of this procedure are presented in [23,55]; here we start with the simplest form for the PDF containing the matching kernel and the leading PDF behavior, which we label as 2-param (Q) Here, K(xν, µ 2 z 2 ) is the matching kernel that factorizes the reduced pseudo-ITD directly to the gluon PDF and the beta function, B(a, b) = 1 0 r a−1 (1 − r) b−1 dr . To assess our fit model, and the associated systematic uncertainties, we add terms to the model. We consider the effect of adding one transformed Jacobi polynomial to the functional form of the PDF and label this model 3-param (Q) , Finally, we consider a model that we denote 2-param (Q) + P 1 for which we add a nuisance term to capture possible O a/|z| effects. This nuisance term can be parametrized by a transformed Jacobi polynomial [23] M(ν, where The transformed Jacobi polynomials, J (α,β) n (x) are defined as, with ω (α,β) n,j = n j (−1) j n! Γ(α + n + 1)Γ(α + β + n + j + 1) Γ(α + β + n + 1)Γ(α + j + 1) .
The transformed Jacobi polynomials form a complete basis of functions in the interval [0,1], making it possible to parameterize the PDF. We use Bayesian analysis to extract the PDF from the reduced pseudo-ITD. We denote the set of fit parameters, which includes the exponents α, β, and the linear coefficients of the Jacobi series for the PDF and additional terms, by θ. Bayes' theorem gives the posterior distribution, P [θ|M, I], which describes the probability distribution of a given set of parameters being the true parameters for a given set of data, M(ν, z 2 ), and prior information, I, as Here, P [M|θ] is the probability distribution of the data for a given set of model parameters. The prior distribution, which describes the probability distribution of a set of parameters given some previously held information, is P [θ|I] and P [M|I] is the marginal likelihood or evidence that describes the probability that the data are correct given the previously held information.
In our parameterization, the PDF is dominated by the leading behavior x α (1 − x) β and the other terms should be small corrections to this. Therefore, in the 3-param (Q) model, our prior for the PDF model parameter, d (α,β) 1 is given by a normal distribution, with a mean and width of d 0 and σ d , respectively. Similarly, in the 2-param (Q) + P 1 model, we expect the parameter for the additional P 1 term to be a small correction to the dominant PDF and use a normal distribution as a prior. The mean and width of the distribution are given by e 0 and σ e .
Guided by phenomenological fits of PDFs, we set α and β to be positive and their prior distributions are set to be log-normal distributions, where µ l is the mean and σ 2 the variance of the distribution of log(x − x 0 ), and x 0 is the lower bound of the log-normal distributions. The most likely parameters of the model are found by maximizing the posterior distribution. This is performed by minimizing the negative log of the posterior distribution, where C is the normalization of the posterior, which is independent of the model parameters. In Fig. 9, we compare the light-cone ITDs obtained from these three models. Adding more terms to the functional form of the PDF or adding more nuisance terms does not improve the quality of the fits and the limited Ioffe-time range does not allow us to add an arbitrary number of parameters to the fit models. Fig. 9 demonstrates that the ITDs do not differ among the three models and the resulting PDFs remain quantitatively the same. We list the L 2 /d.o.f. and χ 2 /d.o.f. of the models in Table IV and find no significant change. The χ 2 /d.o.f. and L 2 /d.o.f. values are also in the acceptable range and their proximity shows that the prior distributions on the PDF parameters do not have a significant effect on the fit. Therefore, for our following discussion, we focus on the 2-param (Q) model.  In Fig. 10, the reduced pseudo-ITD calculated is shown for different separations, z, along with its fitted bands obtained from the 2-param (Q) model. In Fig. 11, we plot the light-cone Ioffe-time distribution with the lattice data modified by the matching kernel from the short distance factorization. SDF removes the logarithmic z 2 dependence of the reduced pseudo-ITD, and introduces the µ 2 dependence on the light-cone Ioffe-time distribution. This effect can be observed in Fig. 11, where after applying the matching kernel, the lattice data points with different field separations shift upward, depending on their field separations, and the data points fall on a regular light-cone Ioffe-time distribution for all z 2 . In previous pseudo-PDF calculations such as the pion valence quark distribution determination [49], the PDF moments extracted by implementing SDF show the logarithmic z 2 dependence removed for z up to 1 fm. Similar results can be found in [47], where the moments of quark distribution in the nucleon calculated through SDF are found to be independent of a logarithmic z 2 effect for z as large as 0.93 fm. On the other hand, if SDF breaks down, we should see a non-polynomial z 2 dependence in the lattice data, especially for large z 2 . We do not see such behavior within the current statistics. Instead, the lattice data, after modification by the matching kernel, aligns with the light-cone Ioffe-time distribution band, including the large z 2 data points, indicating that SDF is quite successful in extracting the Ioffe-time distribution.
In Fig. 12, we present the unpolarized gluon PDF (cyan band) extracted from the 2-param (Q) model (fit Eq. (35)) and compare this with the gluon PDFs extracted from the phenomenological data sets CT18 [3], NNPDF3.1 [4], and JAM20 [87] at µ = 2 GeV. A similar comparison can be made with the other global fits of the gluon PDF, such as with CJ15 [5], HERAPDF2.0 [88], MSHT20 [2]. To determine the normalization of the gluon PDF according to Eq. (15), we need to normalize the extracted PDF with the gluon momentum fraction. There has been a number of lattice calculations to extract the gluon momentum fraction [35,89], as well as phenomenological calculations [3,4]. We take the results from [35], which is x g =0.427 (92) in the MS scheme at renormalization scale µ = 2 GeV, and apply this normalization to our gluon PDF. One could similarly adopt the normalization from the x g determination in [89]. We consider the uncertainties of our extracted gluon PDF and the gluon momentum fraction from [35] to be uncorrelated and determine the total uncertainty in the PDF. The statistical uncertainty of the gluon PDF determined from the fit Eq. (35) and the uncertainty from the normalization using x g are added in quadrature and the final uncertainty is shown as the outer band in Fig. 12.
As discussed in [86], from the fitting of the ITD constructed from the NNPDF x g(x) distribution, one needs the lattice data beyond ν ∼ 15 to evaluate the gluon distribution in the small-x region. In the present calculation, we can  extract the ITD up to ν ∼ 7.07. Therefore, the larger uncertainty and difference in the small-x region determined from the lattice data is expected. As a cautionary remark, we also remind the readers that we have not included the mixing of the gluon operator with the quark singlet sector in the present calculation. Moreover, this calculation is performed at the unphysical pion mass and in principle, physical pion mass, continuum, and infinite volume extrapolation should be performed for a proper comparison with the phenomenological distribution. Therefore, it remains a matter of future investigation to draw a more specific conclusion about the x g(x) distribution extracted from the lattice QCD calculation in the large-x region. We also note that the shrinking of the statistical uncertainty band in the PDF near x ∼ 0.15 results from the correlation of the PDF fit parameters. This feature has also been seen in previous works [28,40,49,51].
However, within these limitations, we find the large-x distribution is in reasonable agreement with the global fits of x g(x) distribution, as can be seen from Fig. 12. The value of β = 5.85 (72) determined in this calculation is statistically in good agreement with the leading (1−x) β behavior obtained in [86] from the fit to the NNPDF3.1 gluon distribution and a recent phenomenological calculation [90]. The I S (ν, µ 2 ) distribution, which we have not included in the present work, is expected to have an increasingly larger effect as ν increases and is expected to have an observable effect in the small-x gluon distribution. However, in the present lattice calculation at heavier up-and down-quark masses, one expects the singlet distribution to increase at a slower rate compared to the phenomenological singlet distribution,   [4], and JAM20 [87]. The normalization of the gluon PDF is performed using the gluon momentum fraction x MS g (µ = 2 GeV)=0.427 (92) from [35]. The figures on left and right are the same distributions with different scales for x g(x) to enhance the view of the large-x region.
therefore having a smaller effect on the Ioffe-time distribution in the 0 ≤ ν ≤ 7.07-range.

VIII. CONCLUSION AND OUTLOOK
In this paper, we present the unpolarized gluon parton distribution using the pseudo-PDF approach. We employ the distillation technique, combined with momentum smearing in our lattice. Distillation allows us not only to improve the sampling of the lattice but also to construct the nucleon two-point correlators with an extended basis of interpolators, which is necessary for the implementation of the sGEVP method. By using momentum smearing, momentum as high as 2.46 GeV is achieved. The sGEVP method combines the features of the summation method and GEVP technique, suppressing the excited-state contributions to the matrix elements significantly. Gradient flow reduces the UV fluctuations from the flowed matrix elements. The combination of these techniques enable us to control the signal-to-noise issues to a great extent. The reduced pseudo-ITD is calculated from the flowed reduced matrix elements by fitting the τ -dependence using a linear form and extrapolating to τ → 0 limit. Using the Jacobi polynomial parameterization, the gluon parton distribution is extracted directly from the reduced pseudo-ITD. Although systematics like higher-twist contributions, lattice spacing errors, infinite volume effects, unphysical pion mass effects are not refined from the parton distribution, and quark-gluon mixing is excluded from the calculation, the resultant ITD has a well-regulated signal-to-noise ratio. The gluon PDF extracted is remarkably consistent with that extracted from the phenomenological distributions. Future endeavors include performing the calculation with a larger number of gauge configurations on the same ensemble and also perform a lattice calculation of the gluon momentum fraction, which will enables us to address the systematic uncertainties more completely along with better statistics. Incorporating the quark-gluon mixing to the calculation is another task we are aiming to undertake. When all the systematic uncertainties are properly quantified and the mixing with the isoscalar quark PDF are included, the lattice calculations will help constrain the gluon PDF at large-x, where the PDF is less constrained by experimental data.

IX. ACKNOWLEDGEMENT
We would like to thank all the members of the HadStruc collaboration for fruitful and stimulating exchanges. TK and RSS acknowledge Luka Leskovec and Archana Radhakrishnan for offering their generous help, which greatly assisted this research. TK is support in part by the Center for Nuclear Femtography grants C2-2020-FEMT-006, C2019-FEMT-002-05. TK, RSS, and KO are supported by U.S. DOE Grant #DE-FG02-04ER41302. AR and WM are also supported by U.S. DOE Grant #DE-FG02-97ER41028. JK is supported by U.S. DOE grant #DE-SC0011941. This work is supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under contract DE-AC05-06OR23177. Computations for this work were carried out in part on facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy. This work was performed in part using computing facilities at The College of William and Mary which were provided by contributions from the National Science Foundation (MRI grant PHY-1626177), and the Commonwealth of Virginia Equipment Trust Fund. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. Specifically, it used the Bridges system, which is supported by NSF award number ACI-1445606, at the Pittsburgh Supercomputing Center (PSC) [91,92]. In addition, this work used resources at NERSC, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract #DE-AC02-05CH11231, as well as resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. #DE-AC05-00OR22725. The software codes Chroma [93], QUDA [94,95] and QPhiX [96] were used in our work. The authors acknowledge support from the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research and Office of Nuclear Physics, Scientific Discovery through Advanced Computing (SciDAC) program, and of the U.S. Department of Energy Exascale Computing Project. The authors also acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources, like Frontera computing system [97] that has contributed to the research results reported within this paper. We acknowledge PRACE (Partnership for Advanced Computing in Europe) for awarding us access to the high performance computing system Marconi100 at CINECA (Consorzio Interuniversitario per il Calcolo Automatico dell'Italia Nord-orientale) under the grant Pra21-5389. JLAB-THY-21-3469.