Gluon PDF of the proton using twisted mass fermions

In this paper, we present lattice QCD results for the $x$-dependence of the unpolarized gluon PDF for the proton. We use one ensemble of $N_f=2+1+1$ maximally twisted mass fermions with a clover improvement, and the Iwasaki improved gluon action. The quark masses are tuned to produce a pion with a mass of 260 MeV. The ensemble has a lattice spacing of $a=0.093$ fm and a spatial extent of 3 fm. We employ the pseudo-distribution approach, which relies on matrix elements of non-local operators that couple to momentum-boosted hadrons. In this work, we use five values of the momentum boost between 0 and 1.67 GeV. The gluon field strength tensors of the non-local operator are connected with straight Wilson lines of varying length $z$. The light-cone Ioffe time distribution (ITD) is extracted utilizing data with $z$ up to 0.56 fm and a quadratic parametrization in terms of the Ioffe time at fixed values of $z$. We explore systematic effects, such as the effect of the stout smearing for the gluon operator, excited states effects, and the dependence on the maximum value of $z$ entering the fits to obtain the gluon PDF. Also, for the first time, the mixing with the quark singlet PDFs is eliminated using matrix elements with non-local quark operators that were previously analyzed within the quasi-PDF framework on the same ensemble. Here, we expand the data set for the quark singlet and reanalyze within the pseudo-PDFs method eliminating the corresponding mixing in the gluon PDF.

In this paper, we present lattice QCD results for the x-dependence of the unpolarized gluon PDF for the proton.We use one ensemble of N f = 2 + 1 + 1 maximally twisted mass fermions with a clover improvement, and the Iwasaki improved gluon action.The quark masses are tuned to produce a pion with a mass of 260 MeV.The ensemble has a lattice spacing of a = 0.093 fm and a spatial extent of 3 fm.We employ the pseudo-distribution approach, which relies on matrix elements of non-local operators that couple to momentum-boosted hadrons.In this work, we use five values of the momentum boost between 0 and 1.67 GeV.The gluon field strength tensors of the non-local operator are connected with straight Wilson lines of varying length z.The light-cone Ioffe time distribution (ITD) is extracted utilizing data with z up to 0.56 fm and a quadratic parametrization in terms of the Ioffe time at fixed values of z.We explore systematic effects, such as the effect of the stout smearing for the gluon operator, excited states effects, and the dependence on the maximum value of z entering the fits to obtain the gluon PDF.Also, for the first time, the mixing with the quark singlet PDFs is eliminated using matrix elements with non-local quark operators that were previously analyzed within the quasi-PDF framework on the same ensemble.Here, we expand the data set for the quark singlet and reanalyze within the pseudo-PDFs method eliminating the corresponding mixing in the gluon PDF.

I. INTRODUCTION
As the mediators of the strong force, gluons play a significant role in the internal structure of hadrons.However, color confinement, a key aspect of quantum chromodynamics (QCD), prevents direct observation of quarks and gluons.Instead, both theoretical and experimental approaches to hadronic structure calculations rely on QCD factorization, which separates the perturbatively-calculable hard-scattering part from the non-perturbative part described by form factors and distribution functions, including parton distribution functions (PDFs).PDFs are probability distributions quantifying the likelihood of finding partons with a particular momentum fraction.Precise and accurate calculations of the gluon PDF are necessary for J/ψ photo production at Jefferson lab, the cross-section of Higgs boson production and jet production at the Large Hadron Collider (LHC), as well as providing theoretical input to experiments at the future Electron-Ion Collider (EIC) in the U.S. and the Electron-Ion Collider in China (EicC).
Lattice QCD is a first-principles approach to calculating strong force quantities performed on a discrete 4dimensional Euclidean lattice.While lattice QCD calculations have proven successful in extracting the nonperturbative dynamics of QCD governing hadron structure, the light-like nature of PDFs prevents direct calculation on Euclidean lattices.Several methods have been proposed over the last decade to relate lattice data to physical light-cone distributions.Two notable and most widely used approaches are the quasi-distribution [1,2] and pseudodistribution [3][4][5][6][7] methods.These approaches utilize the same matrix elements of momentum-boosted hadrons coupled to non-local operators containing a Wilson line but differ in the way the Euclidean observable is factorized into its light-cone counterpart directly in coordinate space (pseudo) or after reconstruction of the x-dependence, i.e. in momentum space (quasi).Typically, they are also renormalized differently.By construction, the renormalization for pseudo-distributions employs canceling the divergences by forming an appropriate ratio of matrix elements (ratio scheme).In turn, quasi-distributions are typically renormalized using a dedicated calculation of vertex functions of the operator under study that leads to an RI/MOM type of renormalization.It should be noted that, the ratio scheme is also increasingly utilized for quasi-distributions in hybrid schemes [8] that treat short and long scales differently.Another typical difference is in the x-dependence reconstruction.For quasi-distributions, this step uses Euclidean matrix elements in the full range of the non-local operator lengths, z.Pseudo-distributions, in turn, are matched in coordinate space, which imposes limitations on the value of z, which needs to be kept relatively small so that it remains in the perturbative region.Thus, without access to the full range of z, approaches based on pseudo-distributions typically employ a physically-motivated fitting ansatz for the functional form of the reconstructed function.
In this work, we present our calculation of the unpolarized gluon PDF for the proton using the pseudo-PDF approach.We calculate the Ioffe-time pseudo-distribution function (pseudo-ITD) by taking the ratio of matrix elements and evolving to a common scale.The ITD describes the interaction of the nucleon with the probe in deep inelastic scattering (DIS) interactions.We use a fitting ansatz to reconstruct the pseudo-PDF from the pseudo-ITD.This approach has proved successful for the extraction of the quark pseudo-PDF.The gluon component presents additional difficulties, including the need for an order of magnitude more statistics arising from the noise associated with the purely disconnected diagram.The gluon PDF also mixes with the quark singlet PDF.Previous lattice calculations have neglected this mixing.We present the first analysis incorporating the quark singlet mixing from lattice QCD data.We compare our pseudo-PDF results neglecting mixing with lattice results from the HadStruc collaboration [66].We also compare our results with and without mixing to global analysis from the JAM collaboration [71].This paper is organized as follows.In Secs.II A and II B, we describe the theoretical and lattice setups for the calculation.In Sec.III A, we present our analysis of various smearing and source-sink time separation values of the matrix elements and reduced-ITDs.Sec.III B shows the results of the pseudo-ITD and pseudo-PDF neglecting mixing with the quark singlet, and Sec.III C presents the results addressing the mixing with the quark singlet.

A. Approach
The computationally expensive component of the methodology is the evaluation of matrix elements with momentumboosted proton states, N (P ), that couple to non-local gluon operators; P indicates the proton momentum.The operator is constructed by two gluon field-strength tensors, F µν , located at two lattice points that are spatially separated in the ẑ direction by distance z.The operator also contains two straight Wilson lines, connecting points 0 → z and z → 0, to ensure gauge invariance.The matrix element reads where F µν is the gluon field strength tensor defined as and g is the bare coupling constant.Potential candidates for the gluon operator are given below for different values of the indices µ , ν, i, j, which can be temporal or spatial, that is The various options of the indices lead to the construction of operators with different properties.Here, we use the operator O 4 , which does not exhibit mixing under renormalization.This operator has a non-vanishing vacuum expectation value that must be subtracted.Since the calculation of the gluon loops is computationally very inexpensive, the vacuum expectation value subtraction does not pose any challenges in the calculation.It should be noted that, regardless of the choice of operator, the unpolarized gluon PDF mixes with the unpolarized singlet quark PDF.We take this mixing into account in our analysis and we quantify its effects by comparing to results with the mixing neglected.
The matrix elements of Eq. (1), M µi;νj , are extracted from the ground state contribution to the ratio taken between the three-point and two-point correlation functions.The variables t s , τ , and t 0 indicate the time of the sink, operator insertion, and source, respectively.Without loss of generality, we have taken the source position to be at t 0 = 0.The ground state contribution, M O ≡ M µi;νj , is identified at large enough values of t s and at τ away from the source and the sink.Practically, we seek convergence with a variance of t s and at τ .
For the calculation of gluonic contributions to the proton, C 3pt O correspond to the so-called disconnected contributions, which are constructed by the expectation value of a product of a gluon loop with the proton two-point function.Also, for the unpolarized gluon PDF, the appropriate parity projector is Γ 0 ≡ 1 4 (1 + γ 0 ) for both the three-and two-point functions.
In our analysis, we implement the pseudo-ITD framework, which requires several nontrivial steps to extract the x-dependence of the gluon PDF.For convenience, we use M g to denote the ground state contribution for the operator O 4 .First, the matrix elements at different values of P and z are combined to construct the reduced Ioffe-time distribution (pseudo-ITD), which depends on the Lorentz-invariant quantities ν ≡ z • P (Ioffe time) and z 2 .For multiplicatively renormalizable operators, the reduced ITD acts as a gauge invariant renormalization scheme that removes UV divergences, including the power divergence due to the presence of the Wilson line.The effects of the residual scale 1/z can be accounted for by an evolution term (see below) and data from different scales 1/z can be combined into ITDs defined at a common renormalization scale, µ 2 .Furthermore, it is anticipated that Eq. ( 9) leads to suppressed discretization and higher-twist effects, which are assumed similar in the two single ratios shown above [44].
Another component of this work is the calculation of the unpolarized quark PDF to address the mixing with the gluon case.The matrix element can be written similarly to the gluon case, that is where the fermionic field ) is taken to be the up, down, and strange quark; f indicates the flavor.For a proper flavor decomposition of the up and down quark contributions, we calculate the disconnected diagram in addition to the connected one.Moreover, the strange-quark contribution is purely disconnected for the nucleon case.Forming the quark-disconnected contributions requires the evaluation of quark loops that are combined with the nucleon two-point correlators.The quark loop of the non-local operator reads where D −1 f (x ins ; x ins + z) is the quark propagator, whose endpoints are connected by a Wilson line.More details in the calculation of the disconnected contributions can be found in Ref. [35].Here we combine the connected and disconnected contributions to the matrix element to form the singlet u + d + s combination, M q .The latter is normalized by constructing the pseudo-ITD, M q similarly to the definition of Eq. ( 9).
To extract the light-cone counterpart of M g , indicated as Q gq , one must apply a matching procedure known to one-loop-level accuracy [74,75], where ⟨x⟩ µ g is the gluon momentum fraction renormalized at the scale µ, and with q f (x, µ 2 ) (q f (x, µ 2 )) being the quark (antiquark) PDF of flavor f , and the sum runs over all considered quark flavors (f = u, d, s).This distribution is related to the imaginary part of the double ratio M q , Differentiating this equation with respect to the upper limit of the integral, we get Thus, the singlet quark Ioffe-time distribution appearing in the matching equation, M S (ν, µ 2 ), is purely real and related to the imaginary part of the quark double ratio.
The matching kernels read and the plus prescription is given by ).The matching equations involve evolving the reduced gluon ITD to a common scale (B gg (u) term), converting the expressions to the light-cone gluon ITD in the MS scheme (L(u) term) and taking its mixing with the singlet quarks into account (B gq term).It is convenient to rewrite Eq. ( 12) in three parts so that one can inspect the role of the three terms separately, where M ′ g (ν, z 2 , µ 2 ) is the evolved gluon ITD, which depends on ν, the final scale µ 2 and the initial scale z 2 .The matching and conversion to the MS scheme is given by Finally, we take the mixing with the singlet quark into account, M q , arriving at the final light-cone ITD, The matched gluon ITD still keeps track of the initial scale z 2 at this stage.However, different scales z 2 should lead to the same light-cone ITDs up to higher-twist effects.For data points where this holds, i.e. leading to consistent values of Q g (ν, z 2 , µ 2 ) from different initial z 2 , Q g is averaged over the same values of ν extracted from different combinations of P and z.We denote such ν-averaged ITDs by Q g/gq (ν, µ 2 ), i.e. dropping the argument indicating the initial scale z.
To extract the x-dependent gluon PDF, xg(x), the light-cone ITDs need to be subjected to a cosine Fourier transform, The extraction of xg(x, µ 2 ) poses an inverse problem [45], because one attempts to calculate a continuous distribution from a limited number of lattice data points for a finite range of Ioffe times up to some ν max .Therefore, to determine xg(x, µ 2 ), one requires additional information, which can be chosen in several ways.Here, we reconstruct the gluon PDF by using a fitting ansatz commonly used in the analysis of experimental data sets, that is where the exponents a, b are fitting parameters and N is the normalization that is fixed by the gluon momentum fraction 1 0 dx xg(x) = ⟨x⟩ g .The lattice data are, thus, fitted according to the minimization of We consider the reconstruction in the cases with (Q gq ) and without (Q g ) the mixing taken into account, to assess the effect of this mixing at the level of the x-dependent distributions.The data are weighted by the inverse variance of the light-cone ITDs, σ 2 Q g/gq (ν,µ 2 ) .Q f (ν, µ 2 ) is the cosine Fourier transform of the assumed fitting ansatz.

B. Setup of lattice calculation
The calculation is performed using an N f = 2 + 1 + 1 ensemble of twisted-mass clover-improved fermions and Iwasaki-improved gluons [76].The quark masses are fixed such that the pion has approximately twice its physical mass (m π = 260 MeV).The lattice spacing is a = 0.0938(2)(3) fm, and the lattice volume is 32 3 × 64 in lattice units.The parameters of the ensemble are summarized in Table I.Matrix elements of gluon operators have increased gauge noise, and one needs to (a) obtain high statistics and (b) use smoothing techniques.To this end, we calculate the correlation functions from different source positions on the same configuration, as the cA211.30.32 ensemble has about 1,200 thermalized gauge configurations [76].Utilizing several source positions per configuration combined with the large speed-up achieved with the use of the multigrid [77][78][79][80], leads to an efficient increase in statistics.Here, we analyze a total of 200 source positions for each configuration.To further increase statistics without loss of generality, we calculate the matrix element using six kinematically equivalent setups, where both the Wilson line and momentum boost are in the ±x, ±y, ±z directions.These six matrix elements can be averaged over, leading to total statistics exceeding one million measurements, as shown in Table II.Since the pseudo-ITD utilizes matrix elements at several values of the proton momentum, we use five values, that is, P = 0, 0.42, 0.83, 1.25, 1.67 GeV.Each matrix element is normalized with the P = 0 case and we found non-negligible correlations between the numerator and denominator of the reduced ITD in Eq. ( 9).These are eliminated by calculating all matrix elements at the same configurations and identical source positions.Regarding excited-state contamination, we use the measurements of Table II   The increased gauge noise is addressed by employing the stout smearing smoothing technique [81] on the gauge links entering the gluon field strength tensor and the Wilson line.The stout smearing parameter is ρ = 0.129 [82,83], and the number of smearing steps is chosen independently in the gluon field strength tensor (N F stout ) and the Wilson line (N W stout ).We apply a 4D smearing to the field strength tensor and a 3D smearing to the gauge links of the Wilson line.We have tested a 4D smearing in the Wilson line, obtaining compatible results after the double-ratio renormalization.We calculate 25 combinations of N F stout and N W stout , by using the values 0, 5, 10, 15, and 20 for each.Another technique to decrease the noise-to-signal ratio is to improve the overlap with the proton ground state.We apply momentum smearing [84] for the three highest momentum boosts, P = 0.83, 1.25, 1.67 GeV, which has been proven essential in suppressing the gauge noise in matrix elements with boosted hadrons and non-local operators [12].We found the optimized value at ξ = 0.6 for the momentum smearing parameter.
The evaluation of the quark matrix elements is an extension of the previous work of Ref. [35], which obtained the quark PDFs within the quasi-distributions method with momenta P = 0.42, 0.83, 1.25 GeV.Here, we added P = 0, 1.67 GeV, so that we obtain the reduced ITDs for the quarks, M q , needed in Eq. (19).As for the gluon case, we implement the momentum smearing method and five stout smearing steps on the gluon fields of the Wilson line entering the operator.To reduce the stochastic noise coming from the low modes [85] in the calculation of the quark loops, we compute the first N ev = 200 eigenpairs of the squared twisted-mass Dirac operator.Then, the low-mode contribution to the all-to-all propagator can be exactly reconstructed, and the high-modes can then be evaluated with stochastic techniques, such as hierarchical probing [86].The latter allows for the reduction of the contamination of the off-diagonal terms in the evaluation of the trace of Eq. ( 11), up to a distance 2 k , using Hadamard vectors as basis vectors for the partitioning of the lattice.Here, we use k = 3 in four dimensions leading to 512 Hadamard vectors.In addition to the hierarchical probing, we make use of the one-end trick [87,88] and fully dilute spin and color sub-spaces.More information can be found in Ref. [35], as well as Refs.[83,[89][90][91].
On a large enough lattice, and given that the source positions are selected randomly, the autocorrelations become very small, and the data on multiple source positions on the same configuration can be considered statistically independent.To check for autocorrelations, we analyze different subsets of data for the two-point functions and extract the relative error on the energy, as shown in Fig. 1 for two representative values of the momentum boost, P = 0.83 GeV (p = 2) and P = 1.67 GeV (p = 4).We find that the statistical error of various quantities scales with the inverse square root of the number of source positions indicating uncorrelated data.

A. Gluon matrix elements and reduced ITDs
Before presenting the final bare matrix elements, it is useful to examine the effect of the stout smearing in terms of the signal quality.The stout smearing is extensively used in the calculation of non-local operators of Refs.[12, 16, 17, 20, 27-29, 35, 36] demonstrating the noise reduction.Also, in Ref. [20], we demonstrated the independence of the renormalized matrix elements from the level of smearing.However, the above statements regard quark bilinear operators, so similar tests are imperative for gluonic operators.As mentioned in the previous section, we construct the gluon matrix elements for 25 combinations of stout steps in the gluon field strength tensor and the Wilson line, that is 20] in steps of 5.The bare matrix elements are shown in Fig. 2 for a subset of these combinations, which includes {N F stout , N W stout } = {0, 10, 20}.All presented matrix elements have been evaluated at t s = 9a, which, as we will demonstrate below, is the one used in the final analysis.It is interesting to observe that the smearing on the field strength tensor has a bigger impact on the signal compared to the smearing on the Wilson line.For instance, the signal already improves significantly with N F stout = 10 and N W stout = 0. Comparing the effect of the stout smearing directly at each momentum can offer another qualitative understanding of signal improvement.In Fig. 3 we show selected cases of the N W stout and N F stout combinations, presented as N stout = (N F stout , N W stout ).As previously discussed, the stout smearing applied on the gluon fields of the field strength tensor is crucial to get a signal.In all values of P , further signal improvement is found as N F stout and N W stout increase.We observe a saturation at (N F stout , N W stout ) = (20, 10), which we will use for the remainder of this analysis.In Fig. 6, we will examine the effect of the mixing in the pseudo ITDs.In this work, we also examine excited-states effects using the preferred setup for the stout smearing, N F stout = 20, N W stout = 10.Fig. 4 shows the matrix elements at four values of the source-sink time separation, that is, t s = 8a, 9a, 10a, 11a.For P = 0, 0.42, and 0.83 GeV, there is an indication of excited-states effects at t s = 8a, which differs from t s = 10a and t s = 11a.The effect is visible mainly due to the high statistical accuracy of the data.For higher momenta, all matrix elements are compatible within uncertainties, which are enhanced compared to the lower momenta.Therefore, t s = 9a is favorable, as it is consistent with t s = 10a and t s = 11a in all cases, while good signal is maintained.Below, we will also consider excited-states effects in the reduced ITDs (see, e.g., Fig. 6).L p with p = 0, 1, 2, 3, 4 are shown with blue squares, red circles, green downward-pointing triangles, yellow upward-pointing triangles, and magenta rightward-pointing triangles, respectively.
To summarize the presentation of the bare matrix elements, we compare in Fig. 5 the data for all values of the momentum boost using t s = 9a, N F stout = 20, and N W stout = 10.The P dependence of the data is as observed in the quark case, that is, the signal quality decreases.We find that the relative error at z = 0 for P = 0 is about 6%, while for P = 1.67 GeV, the error becomes close to 9% despite the same statistics.In all cases, we find that the matrix elements decay to zero at about z = 8a.
The matrix elements of Fig. 5 are the core of our calculation and are used to construct the double ratio of Eq. ( 9).We note that systematic uncertainties might affect the matrix elements and pseudo ITDs differently due to possible correlations between the ratios in the numerator and/or denominator.Thus, investigating systematic effects, such as excited states and stout smearing, in the ratio of Eq. ( 9) is important.Since the pseudo ITD, i.e. the double ratio, serves as a renormalization prescription, it should be independent of the number of smearing steps.We examine the validity of this argument, and a summary is shown in Fig. 6.Due to the large uncertainties of certain combinations of smearing steps, their inclusion in the plot is not meaningful, as their errors cover the whole range of the plot.As can be seen, all combinations of N F stout and N W stout are in full agreement within errors, demonstrating that the pseudo ITDs can be extracted from any of these combinations.(15,15), (20,10), (20,20) are shown with blue circles, green squares, red up triangles, black down triangles, and orange left-pointing triangles, respectively.Right: Excited states in the reduced ITDs at N F stout = 20 and N W stout = 10.ts = 8a, 9a, 10a, 11a are shown with blue circles, green squares, red up triangles, and orange down triangles, respectively.We study excited-states effects in M g , as shown in Fig. 6 for four values of the source-sink time separation, that is t s = 8a, 9a, 10a, 11a.The increase of the statistical error is sizeable between t s = 9a and t s = 11a and the signal is lost at t s = 12a; the latter is not shown here.Overall, we find that both t s = 8a and t s = 9a are good options for these data.Nevertheless, we choose t s = 9a for a more conservative estimate.For completeness, we show in Fig. 7 the double ratio for all values of P corresponding to t s = 9a, N F stout = 20 and N W stout = 10.We note that each value of ν is constructed from all possible combinations of available z and P , but we constrain z up to 6a ∼ 0.56 fm.This leads to ν max ∼ 5. Comparing all combinations of P and z at a given ν allows one to comment on the effect of P dependence.We find that dependence on P is within the statistical errors for up to z = 6a.Thus, M g can be described by a smooth function in terms of the Ioffe time, which allows for a controlled interpolation.The latter is needed for the scale evolution and matching procedure, as discussed below.
The reduced ITDs are interpolated in terms of ν at each value of z 2 , so that one obtains a continuous function in ν/z, which is needed for the matching procedure.Having five values of ν at a fixed z 2 , one can test different parametrizations of the ν dependence.Here, we test a linear and a second-order polynomial fit, which can be seen in Fig. 8 for selected values of z.We find that the two fits are compatible and choose the polynomial fit to proceed.0.0 0.5 0.0 0.2 0.4 0.6 FIG. 8. Lattice data of the reduced ITDs for z = 1a − 6a (blue points) and their interpolation at fixed z 2 using a first order (green bands) and second order polynomial fit (red bands).

B. Reconstruction of the gluon PDF
In general, the extraction of light-cone ITDs from reduced ITDs contains combining effects of three functions, the B gg , L and B gq kernels, as given in Eq. (12).That is, one must apply the evolution to a renormalization scale of choice (µ), convert the data to light-cone ITDs in the MS scheme, and eliminate the mixing with the quark-singlet PDF.
Here, we chose 2 GeV for the renormalization scale, as commonly used in global analyses.In all previous calculations of the gluon PDF, the mixing with the quark-singlet case has been ignored due to the lack of lattice results for the latter, as it requires information from disconnected contributions, which are computationally very expensive.Here, we extend the calculation of Ref. [35] to include all values of P implemented in this work, which allows us to eliminate the mixing by considering B gq .To demonstrate the effect of the mixing, we first apply B gg and L, but ignore B gq .The resulting evolved and matched ITDs are shown in Fig. 9.We find that the scale evolution increases the values of the evolved ITDs (M ′ ) relative to those of the reduced ITDs (M), while the matching has the opposite effect than the evolution and brings the light-cone ITDs (Q) closer to the reduced ITDs, making them consistent with the latter within error bars.Such behavior is also observed in the case of quark PDFs (see, e.g., Refs.[50,55]).We note that the dependence on the individual P and z is minimal for all three functions, M, M ′ , and Q, as the values from different (P, z) pairs fall on a universal curve.In the right panel of Fig. 9, we show the matched ITDs, Q(z 2 , µ 2 ), where we average over the (P, z) pairs for a given value of the Ioffe time.
To extract the x-dependence of the gluon PDF, we use the fitting reconstruction and follow the procedure discussed in Sec.II A. As can be seen in Eq. ( 21), one cannot isolate the gluon PDF, because it appears normalized with the gluon momentum fraction.The latter has not been extracted on the ensemble under study, so we use, instead, the lattice results of Ref. [83].The aforementioned calculation used an ensemble that has the same gluon and fermion action as this work, but different lattice parameters.In particular, the lattice spacing is 0.08 fm, and the pion mass is 139 MeV.The reported value for the gluon momentum fraction is ⟨x⟩ MS, 2GeV g = 0.427 (92), which we use below.Another input of the reconstruction procedure is the value of z max .We tested z max = 5a, 6a, 7a, which correspond to z max = 0.47, 0.56, 0.66 fm, respectively.In this subsection, we show results for z max = 6a and we demonstrate the independence of the results on z max in the case where mixing with the quark singlet is eliminated (see next subsection).The conclusions of this test fully pertain also to the case with the mixing neglected.
Using the above value of ⟨x⟩ MS, 2GeV g and z max = 6a, we obtain the gluon PDF, which is given in the left panel of Fig. 10.The corresponding fitted ITDs are shown in the right panel of Fig. 9.We remind the reader that we have not yet considered the mixing with the quark-singlet PDF; this will be addressed in the next subsection.In the right panel of Fig. 10, we compare our final results to the lattice results of HadStruc [66], in which the gluonquark singlet mixing has not been considered.HadStruc used an ensemble of N f = 2 + 1 clover Wilson fermions with stout-link smearing and the Symanzik-improved gauge action.The ensemble has the same volume and lattice spacing as this work.However, their pion is heavier, namely m π = 358 MeV.Their source-sink time separation is also 9a, which is the same as the value used here.In general, our results are consistent with the ones from HadStruc.It is worth noting that the reconstruction performed by HadStruc includes values of Ioffe time up to ν max = 7.07, while our reconstruction includes up to a maximum Ioffe time of ν max = 4.71 (z max = 6a).The smaller statistical error of HadStruc may possibly be attributed to two factors: (a) the use of the distillation method [92]; (b) the higher pion mass compared to this work.In Fig. 10, we also compare the lattice data to the global analysis of JAM20 [70].As can be seen, all results are in full agreement within errors.We note that all comparisons are qualitative, as the lattice results are obtained on a single ensemble with different lattice formulations.Nevertheless, the agreement between lattice results and global analysis is very promising.C. Elimination of mixing with quark-singlet PDF In this section, we provide, for the first time, the quark-singlet PDF using the pseudo-distribution method.This is a continuation of the work of Ref. [35], which used a subset of the data of Table II to obtain the quark PDFs within the quasi-distributions method.The quark-singlet PDF will be used to eliminate the mixing with the gluon contribution using the matching formalism of Ref. [93].In principle, with our data for the quark and gluon PDFs, we can also obtain the quark-singlet PDFs without mixing.However, Refs.[74,75] only provide the components of the mixing kernel that are relevant to the gluon PDF, that is, B gg and B gq .While the complete 2 × 2 matching kernel is presented in Ref. [93], it corresponds to a different definition of the gluon operator than the one we use in this work, so we are not able to apply it here.
First, let us present the bare quark matrix elements for the singlet combination u + d + s.The matrix elements contain all kinematic factors, so they can be compared directly at z = 0 for different values of P .As seen in Fig. 11, the data are consistent at z = 0.This is expected theoretically, because z = 0 is directly related to ⟨x⟩, which is independent of the kinematic frame.As z increases, we find that the behavior with the increase of P is as expected.That is, the real part of the matrix element falls faster, while the imaginary part is enhanced.The corresponding quark reduced ITDs are shown in Fig. 12.In the real part, we find an agreement between the different P and z combinations corresponding to the same value of ν.Some difference is observed in the imaginary part for {p, z/a} = {1, 4} as compared to {p, z/a} = {2, 2} and {p, z/a} = {4, 1} (where P = 2π L p).Similarly, {p, z/a} = {1, 6} deviates from {p, z/a} = {2, 3}, and {p, z/a} = {3, 2}.However, the momenta with p > 1 are in agreement within errors.We use the above quark-singlet reduced ITDs to eliminate the mixing in the light-cone gluon ITDs.In particular, only M S -the ν-derivative of the imaginary part of M q -enters the matching formalism, as explained in Sec.II A. For completeness we show M S in Fig. 13.Finally, the resulting effect of the mixing is shown in Fig. 14  gluon ITDs before (B gq = 0) and after (B gq ̸ = 0) the elimination of the mixing with the quark-singlet.The fitting bands from the x-dependence reconstruction procedure are also shown.The main finding is that the gluon ITDs move slightly towards lower values, with the mixing being well within the occurring statistical uncertainties.As hinted in the previous subsection, we also establish the robustness of the results against the choice of the value of z max , using z max = 5a, 6a, 7a, see Fig. 15.We find a small difference between z max = 5a and z max = 6a, but the effect is significantly smaller than the statistical uncertainties.The difference z max = 6a and z max = 7a is almost negligible, and the two bands cannot be visually distinguished in the figure.That is, the addition of z = 7a points does not influence the reconstruction, mainly due to their large statistical errors.Thus, our choice of z max = 6a is validated at this level of data precision and given the compatibility of results for different z max , it is not necessary to assign a reconstruction-related systematic uncertainty to our results.
For completeness, we also present the effect of the mixing in the x-dependent gluon PDF, as seen in the right panel of Fig. 16.The conclusion is consistent with Fig. 14, as the effect of the mixing is smaller than the statistical uncertainties.In the left panel of Fig. 16, we show our final results together with JAM20 [70], demonstrating full compatibility.As previously mentioned, the statistical uncertainties are currently larger than the ones from global analysis.

IV. SUMMARY
The main component of this work is the calculation of the unpolarized gluon PDF of the proton using numerical simulations of QCD.The calculation is performed using an N f = 2 + 1 + 1 ensemble of clover-improved twisted mass fermions with the quark masses tuned to give a pion mass of 260 MeV.The lattice spacing is 0.093 fm, and the volume is 32 3 × 64.For the calculation, we employ the pseudo-distribution approach that significantly simplifies the renormalization procedure by forming ratios of matrix elements, leading to the reduced pseudo-Ioffe time distributions, expressed in terms of the combination ν = z • P .In our calculation, we use a nucleon momentum boost with values up to 1.67 GeV, and, in the final results we restrict the length of the Wilson line to 0.56 fm.We find that the combination of P and z suffices to extract a continuous dependence on ν and reconstruct the gluon PDF.We explore systematic effects such as excited-states effects, the effect of stout smearing, and the dependence on the maximum value of z entering the fits to obtain the ITD.For the evolution and conversion to the MS scheme at a scale of 2 GeV, we use a one-loop formalism.We use the fitting reconstruction method to address the inverse problem and obtain the x-dependence of the gluon PDF.A novel aspect of the calculation is the elimination of the mixing with the quarksinglet unpolarized PDF, which we extract for the same ensemble.The effect of the mixing brings the gluon ITD to smaller values, but the effect is much smaller than the statistical uncertainties.However, when the precision stage is reached for such lattice calculations, mixing will inevitably become a more effect.Our results are compared with other lattice data obtained using a different lattice formulation, methodology, and setup [66], and we find a very good agreement.In such a comparison, we ignore the quark-gluon mixing for a more appropriate comparison with Ref. [66].Furthermore, a comparison of our final data with the global analysis of the JAM collaboration [70] reveals agreement, with the global analysis being much more accurate than lattice data at this stage.The above-mentioned comparison uses our data after the elimination of the quark-gluon mixing as done in JAM20.An extension of this work is the investigation of other sources of systematic uncertainties, such as volume and discretization effects, as well as pion mass dependence.In near future, we will address the continuum limit by adding two ensembles with smaller lattice spacing.

4 FIG. 1 .
FIG. 1.The relative error of the proton energy at momentum boost P = 2π L p as a function of the source positions analyzed.As examples, we show p = 2 (blue circles) and p = 4 (red squares).The lines correspond to the 1/ √ Nsrc scaling.

4 FIG. 4 .
FIG. 4. Source-sink time separation dependence of the bare matrix elements at each momentum boost.We use N W stout = 10 and N F stout = 20 in all cases.The results for ts = 8a, 9a, 10a, 11a are shown in blue circles, green squares, red up triangles, and orange down triangles, respectively.The data at momentum boost P = 2π L p with p = 0, 1, 2, 3, 4 are shown in the top, middle left, middle right, bottom left, and bottom right panels, respectively.

4 FIG. 5 .
FIG. 5. Matrix elements of Eq. (1) as a function of the length of the Wilson line, z/a.The data at momentum boost P = 2πL p with p = 0, 1, 2, 3, 4 are shown with blue squares, red circles, green downward-pointing triangles, yellow upward-pointing triangles, and magenta rightward-pointing triangles, respectively.

4 FIG. 7 .
FIG. 7. Final results for the reduced-ITD at ts = 9a, N F stout = 20 and N W stout = 10.We show data with all values of P = 2π L p and z up to 6a ∼ 0.56 fm.Data for p = 1, 2, 3, 4 are shown with blue circles, red down triangles, green up triangles, and magenta right triangles, respectively.

FIG. 10 .
FIG.10.Left: The reconstructed gluon PDF without mixing elimination.Right: A comparison of our results from the left panel (red), the lattice results of HadStruc[66] (green), and the global analysis of JAM20[70] (blue).Results are shown in the MS scheme at a scale of 2 GeV.

FIG. 11 .
FIG. 11.Bare matrix elements for the quark-singlet case as a function of the length of the Wilson line, z/a.The data at momentum boost P = 2π L p with p = 0, 1, 2, 3, 4 are shown with blue squares, green circles, red downward-pointing triangles, yellow upward-pointing triangles, and magenta rightward-pointing triangles, respectively.

Q(ν, µ 2 )FIG. 14 .
FIG.14.Comparison of the light-cone ITD before (Bgq = 0, shown in red) and after (Bgq ̸ = 0, shown in green) the elimination of the mixing with the quark-singlet case.The bands correspond to the fits of the lattice data.

TABLE I .
Parameters of the ensemble used in this work.
at multiple t s values.This comes at no additional computational cost, as, by construction, disconnected contributions are evaluated at open sink time.

TABLE II
. Total statistics for the calculation for each value of P .N confs is the number of configurations, Nsrc the number of source positions, N dir is the number of spatial directions for the Wilson line and P , and Nmeas is the number of total measurements (Nmeas = N confs × Nsrc × N dir ).