Towards Systematic Control of the Nucleon Gluon PDF with the Pseudo-PDF Method

Introduction

Quantum chromodynamics (QCD), the non-Abelian gauge theory of quarks and gluons, derives much of its richness and complexity from the gluon sector, whose self-interactions distinguish it from the relative simplicity of electrodynamics. Although absent from the simplest valence-quark picture of the nucleon, gluons play a central role in determining its internal structure. They carry approximately half of the nucleon’s longitudinal momentum at typical high-energy reference scales and are intimately connected to confinement and the origin of visible mass in the universe [1, 2]. The leading-twist, unpolarized gluon parton distribution function (PDF), $g(x,Q^2)$, characterizes the distribution of gluon longitudinal momentum fractions $x$ within the nucleon at a resolution scale $Q^2$. Understanding this distribution is therefore essential to explicating the role of gluons in nucleon structure, including their contributions to its mass and spin as described by matrix elements of the QCD energy-momentum tensor [3, 4].

Despite its importance, the gluon PDF remains a challenging quantity determine precisely. Since gluons do not couple directly to the virtual photon in electromagnetic deep-inelastic scattering (DIS), their distribution is constrained indirectly through the $Q^2$ evolution of quark structure functions and through complementary measurements of gluon-sensitive processes, including inclusive jet, heavy-quark, and prompt-photon production [5, 6]. Furthermore, accurate knowledge of the gluon PDF is particularly important for hadron-collider phenomenology, where it governs predictions for processes such as Higgs boson production through gluon fusion, top-antitop pair production, and high-energy dijet production. Uncertainties in the gluon distribution consequently propagate into precision tests of the Standard Model and searches for new physics [7, 8]. These uncertainties are especially pronounced at small $x$, where the gluon distribution rises rapidly and nonlinear phenomena such as gluon saturation may become relevant, and at large $x$, where experimental constraints remain limited. Improving the determination of the gluon PDF is thus a continuing priority for both experimental and theoretical studies [9–16].

Within the last decade, lattice QCD has provided a complementary, first-principles approach to determining parton distributions. In particular, large-momentum effective theory (LaMET) [17] and the pseudo-PDF method [18] have enabled access to the longitudinal momentum dependence of PDFs through spatially separated operators evaluated in Euclidean space-time. While these approaches have yielded substantial progress in the determination of quark distributions, gluon PDFs remain comparatively difficult to extract due to the poor signal-to-noise ratio of gluonic matrix elements. Nevertheless, pseudo-PDF calculations have established a growing body of results for the gluon distributions of the pion [19], kaon [20, 21], and nucleon [22, 23]. More recently, advances in lattice calculations have also enabled the first extractions of gluon PDFs using LaMET [24–27]. These developments demonstrate the growing potential of lattice QCD to constrain gluon structure, while highlighting the need for improved systematic control at various stages of the analyses.

In particular, the extraction of nucleon gluon PDFs requires careful control of discretization effects, excited-state contamination, and systematic uncertainties in the PDF reconstruction. These challenges are particularly pronounced for gluonic matrix elements, which are disconnected from the valence–quark lines of the nucleon interpolating fields and are extracted from correlations between the nucleon two-point function and a gluonic operator containing a spatial Wilson line. The resulting signal-to-noise ratio deteriorates with increasing source-sink separation and nucleon boost momentum, limiting the Euclidean time range over which statistically reliable data can be obtained. Excited-state contamination can therefore remain significant, complicating the isolation of the ground-state matrix element. Moreover, conventional multi-state fits can exhibit sensitivity to the chosen fit window and number of states, making robust control of these effects essential for reliable PDF extraction. In this work, we address these challenges through a combination of fine lattice spacing, explicit control of excited-state contamination, and systematic treatment of the PDF reconstruction. To constrain excited-state effects, we employ the Lanczos method [28], recently applied to a variety of lattice QCD spectroscopy and matrix elements extractions [29–35]. Building on these developments, we combine this excited-state control of fine-lattice-spacing matrix elements and a systematic treatment of PDF extraction using parameterizations previously used in quark pseduo-PDF extractions [36, 37] to assess the reliability of the resulting nucleon gluon PDFs. The remainder of this paper is organized as follows. In § 2, we describe the lattice setup and matrix-element analysis, including our implementation of the Lanczos method. In Sec. III, we present the gluon PDF extraction and discuss our results. Finally, in § 4, we summarize our findings and outline directions for future improvements.

Lattice Details and Analysis Method

We carry out a lattice QCD calculation of the unpolarized gluon PDF of the nucleon at a physical pion mass on an $N_f = 2 + 1 + 1$ ensemble of highly improved staggered quarks (HISQ) [38] generated by the MILC Collaboration [39–41]. The ensemble has a lattice spacing of $a \approx0.0582(4)$ fm and a volume of $L^3 \times T = 48 \times144$. The results presented in this work are obtained using 676, 512 nucleon two-point correlator measurements. We use Wilson–clover fermions in the valence sector, with the valence-quark masses tuned to reproduce the masses of the light and strange pseudoscalar mesons in the sea sector, corresponding to pseudoscalar masses of approximately 310 and 690 MeV, respectively. Such a mixed-action setup has been extensively studied in previous works [43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56] [42] and found to be free of exceptional configurations, which could otherwise render the Dirac operator ill-conditioned or lead to anomalously large correlation functions.

The nucleon two-point (2pt) correlation functions are computed on each gauge configuration using multiple source locations and are defined as

where $P_z$ denotes the nucleon boost momentum along the $\hat{z}$-direction and $t$ is the Euclidean time separation between the source and sink. The nucleon interpolating operator is given by $\chi= \epsilon_{abc}(d_a^T C\gamma_5\gamma_{\mathrm{NR}}u_b)\mathcal{P}_+u_c$, where $C$ is the charge-conjugation matrix and the projector $\mathcal{P}_+$ is applied to isolate the desired positive-parity nucleon states. The so-called non-relativistic projector $\Gamma_{\mathrm{NR}} = (1+\gamma_4)/2$ improves the signal [57]; though more intentional interpolator choices to improve large-momentum signal have been recently developed [58] [34]. To improve the overlap with boosted nucleon states, we employ Gaussian momentum smearing [59] with momentum-smearing parameter $ak = 5.6$, which introduces a momentum-dependent phase into the spatial smearing kernel to bias the quark fields toward the target momentum $P_z$, enhancing the signal quality of the corresponding correlation functions up to $P_z = 10 \times2\pi/L \approx4.4$ GeV. Using these momentum-smeared interpolating operators, we construct the corresponding three-point (3pt) correlation functions by inserting a nonlocal gluon operator between the nucleon source and sink. These correlators provide access to the gluon structure of the boosted nucleon and are defined as

where $t_{\mathrm{sep}}$ is the source–sink separation and $t$ is the operator insertion time. In practice, the gluonic insertion is evaluated through contractions of the nucleon two-point functions with the corresponding gluon-loop contributions. The unpolarized gluon PDF operator, first introduced in Ref. [18], is defined as

where $F^{\mu\nu}$ is the gluon field-strength tensor and $W(z,0)$ is a Wilson line connecting the field-strength insertions along the spatial $\hat{z}$-direction. The operator indices $(i,j)$ are summed over the relevant spatial directions. To reduce ultraviolet fluctuations in the gluonic matrix elements, we apply Wilson flow [60] to the gauge links before evaluating the gluon-loop contributions. The flowed gauge links $V_\mu(x,\tau)$ evolve from the original links $U_\mu(x)$ according to the differential equation

where $\tau$ denotes the Wilson flow time, $S_W$ is the Wilson gauge action, and $\partial_{x,\mu}$ is the gauge-covariant derivative with respect to the link $V_\mu(x,\tau)$. The flow smooths the gauge fields by suppressing short-distance fluctuations, thereby reducing ultraviolet noise in the gluon-loop contributions. We integrate the flow equation numerically using 100 steps at flow times $\tau= 3a^2$, denoted as WILSON3.

Analysis Methods

We extract ground-state matrix elements from the nucleon two-point correlator $C_{2\mathrm{pt}}(P_z;t)$ and the disconnected three-point correlator $C_{3\mathrm{pt}}(z,P_z;t_{\mathrm{sep}},t_{\mathrm{ins}})$, where $t_{\mathrm{sep}}$ is the source–sink separation and $t_{\mathrm{ins}}$ the insertion time of the gluon operator. Their spectral decompositions are

with $A_n$ the overlap of the nucleon interpolator with state $n$ of energy $E_n(P_z)$. The bare ground-state matrix element $M(z,P_z) \equiv\langle0(P_z)|O_g(z)|0(P_z)\rangle$ represents the target of the correlator analysis. Traditionally, MSULat has extracted these matrix elements by fitting model parameters defined by the two-state truncation of the spectral decompositions in Eqs. (5) and (6). In this study, we use the Lanczos method for the first time and compare with determination of these matrix elements from a two-state fits to validate the newly applied methodology.

a. Two-State Fits. As an independent cross-check of the Lanczos spectrum, we also fit $C_{2\mathrm{pt}}(P_z;t)$ directly to the two-state truncation of Eq. (5) (shaded bands in Figure 1 and open symbols in Figure 2),

where $E_0$ and $E_1 = E_0 + \Delta E$ are the ground- and excited-state energies, respectively, while $A_0$ and $A_1$ denote the corresponding amplitudes. The four fit parameters $(A_0,E_0,A_1,\Delta E)$ are treated as free parameters, with the energy gap constrained to $\Delta E \equiv E_1-E_0 > 0$. For each momentum $P_z$, we consider candidate fit windows $W = [t_{\min},t_{\max}]$ consisting of contiguous time slices within $[1,t_{\mathrm{snr}}]$, with at least eight points. Here, $t_{\mathrm{snr}}$ denotes the last time slice in the contiguous interval starting at $t=1$ for which the two-point correlator satisfies $|C_{2\mathrm{pt}}(P_z;t)|/\sigma_{C_{2\mathrm{pt}}}(P_z;t) \ge2$. Thus, $N_{\max}=t_{\mathrm{snr}}$ is the total number of time slices satisfying this signal-to-noise criterion in the initial interval. Each window is fit by an uncorrelated $\chi^2(W)$ built from the diagonal jackknife uncertainties of $C_{2\mathrm{pt}}(P_z;t)$. Fits with $\chi^2(W)/\mathrm{dof}>3$, or with $E_0$ or $\Delta E$ pinned at the lower bound, are discarded.

The surviving windows are scored with the Bayesian Akaike information criterion (BAIC) [61],

where $2[N_{\max}-N_{\mathrm{pts}}(W)]$ penalizes each time slice excluded from $W$. We select the highest-BAIC-weight window for which a leave-one-out jackknife refit satisfies the stability criteria:

We then quote the jackknife mean and uncertainty of the fit parameters obtained from the selected window.

b. Lanczos. We apply Lanczos as described in Refs. [62] [29, 30] to the two- and three-point correlators. Regarding the two-point correlator, we start the Krylov space at $t_0=2a$, discarding the contact-dominated time slices $t/a \in\{0,1\}$, giving a maximal Krylov dimension of $m_{\max}=(n_t-t_0/a)/2=17$. Spurious Ritz values are removed using the ZCW overlap test of Ref. [62]. The threshold $\varepsilon_{\mathrm{ZCW}}$ is defined as the minimum ZCW overlap at the last iteration for which at least 90% of bootstrap samples have entirely real spectra, divided by $\varepsilon_{\mathrm{ZCW}}=10$. Surviving states are then ordered by energy. For the three-point function, we truncate the filtered Ritz basis to its two lowest states and extract $M(z,P_z)$ by projecting $C_{3\mathrm{pt}}$ onto the ground state within this two-state subspace. Requiring both the source–operator and operator–sink separations to be at least $2a$ restricts the total source–sink separation to $t_{\mathrm{sep}}/a \in[4,34]$.

All Lanczos estimates are computed on correlated nested bootstrap resamplings of the two- and three-point functions with $N_b^{\mathrm{inner}}=N_b^{\mathrm{outer}}=200$ samples. Central values are nested medians and uncertainties are central 68% bootstrap confidence intervals [29, 30]. The ground-state energies are quoted over the final $N_m=3$ iterations, while the matrix elements are quoted at $m=17$.

Two-point results

We first check the consistency of the Lanczos and two-state fit extractions of ground state energy from the 2pt correlators by inspecting the effective mass

where $E(P_z)$ is the ground state energy of the nucleon with momentum $P_z$. Figure 1 shows the effective masses for the nucleon constructed with both light (right) and strange (left) valence quarks at lattice momenta corresponding to $aP_z=2\pi n_{p_z}/L$ for integer $n_{p_z}\in[0,10]$. At early times, $aE_{\mathrm{eff}}(t)$ decreases as excited-state contributions decay, eventually approaching a plateau consistent with the Lanczos ground-state energy at each momentum.

As $P_z$ increases and the signal-to-noise ratio decays earlier in Euclidean time, the BAIC-selected two-state fit window shifts toward smaller times: from $[11, 35]$ to $[5, 15]$ for the strange nucleon and from $[10, 35]$ to $[4, 14]$ for the light nucleon as $P_z$ increases from 0 to 10. In contrast, Lanczos uses the full correlator from $t_0 = 2a$ at every momentum, with no fit range to choose. The two determinations of $aE_0(P_z)$ have comparable statistical uncertainties and agree within $1.5\sigma$ at all eleven momenta for the strange nucleon. For the light nucleon, they agree within $0.6\sigma$ at every momentum, except $P_z = 5$, where the Lanczos value $aE_0 = 0.793(6)$ lies $2.1\sigma$ above the two-state fit value $0.777(4)$. At zero momentum, the Lanczos first-excited-state gaps are consistent with the two-state-fit results: $a\Delta E = 0.225(13)$ versus $a\Delta E = 0.247(4)$ for the strange nucleon, and $0.258(18)$ versus $0.276(14)$ for the light nucleon.

Figure 2 compares the ground-state energies from both extractions against the continuum dispersion curve

$E_{\mathrm{disp}}(P_z) = \sqrt{(M)^2 + (P_z)^2}$ for both the strange (left) and light (right) nucleons, evaluated using Lanczos rest masses $aM = 0.45349(49)$ and $0.3331(15)$ respectively. For the strange nucleon, both extractions follow the continuum curve over the full range of momenta: the ratio $E_0(P_z)/E_{\mathrm{disp}}(P_z)$ remains within 3% of unity for both Lanczos and the two-state fits. The light nucleon deviates further from the continuum curve. Both extractions lie 6–8% above it at $P_z = 4$–6, with the largest ratios 1.079(8) for Lanczos at $P_z = 5$ and 1.061(6) for the two-state fit at $P_z = 6$, before returning within 1% of unity at $P_z = 10$. This larger deviation, which is common to both extraction methods, may reflect lattice artifacts such as discretization or finite volume effects. We do not attempt to disentangle these in the present work. In either case, the deviation from the continuum dispersion relation remains at the level of a few percent over the full range of momenta, indicating that such effects are small at boost.

Three-Point Results

Next, we can compare the two-state and Lanczos extraction of matrix elements from the three-point data. A common way to visualize the quality of this extraction is to look at the ratio between the three- and two-point data

Figure 3 displays this ratio for representative values of $z$ and $P_z$ for the strange (left column) and light (right column) nucleons. At low boosts (top row), the ratios exhibit a plateau near $t_{\mathrm{ins}} = t_{\mathrm{sep}}/2$ for the displayed source–sink separations, and the Lanczos and two-state-fit results agree within $1\sigma$. For the strange nucleon at $(P_z,z) = (2,3)$, we obtain $M = 0.256(25)$ and $0.268(42)$, respectively, while for the light nucleon at $(P_z,z) = (1,1)$, we obtain $M = 0.291(26)$ and $0.322(41)$. At higher boost (bottom row), the signal-to-noise ratio decays much earlier, restricting the usable separation times to $t_{\mathrm{sep}}/a \lesssim12$. Over this range, the ratio continues rising with $t_{\mathrm{sep}}$, with even the largest displayed separation lies 20–30% below the extracted matrix element, so both methods must extrapolate to the ground state. The two extractions again agree within $1\sigma$, with $M = 0.430(44)$ and $0.465(54)$ for the strange nucleon at $(8,1)$ and $0.393(32)$ and $0.453(54)$ for the light nucleon at $(4,0)$.

Results and Discussions

The bare matrix elements extracted in the previous section are multiplicatively renormalizable as [63]

where $Z_L(1/a)$ renormalizes the logarithmic UV divergences the lattice cutoff $1/a$, and the exponential term renormalizes the linear divergence associated with the Wilson line self-energy, $\delta m$. In the pseudo-PDF method, the most convenient renormalization scheme is to construct the reduced pseudo-Ioffe time distribution (RpITD):

The term $M^{B}(0,0)/M^{B}(0,P_z)$ removes kinematic factors coming from state normalization and the gluon operator definition, while division by $M^{B}(z,0)$ cancels the $z$-dependent renormalization factors. The RpITD is then a function of the Lorentz invariant Ioffe time $\nu= zP_z$ and separation distance squared, $z^2$. The pseudo-PDF methodology is based on matching between the pseudo-ITD and the light-cone ITD, which is related to the gluon PDF just through a Fourier transform. In practice, the two steps of this matching can be written as a single relation [64] [18]:

where $g(x)$ and $\Sigma(x)$ are the $\overline{\mathrm{MS}}$ scheme gluon and quark singlet PDFs respectively. The first moments of these two distributions are represented by $\langle x_g\rangle$ and $\langle x_\Sigma\rangle$. Higher twist effects come into play in this matching at $\mathcal{O}(m^2z^2,z^2\Lambda_{\mathrm{QCD}}^2)$, and are suppressed here by ensuring that we use separation distances within the perturbative range $z\lesssim0.3-0.4$ fm. The terms $R_{gg}$, $R_r$, and $R_{gq}$ are perturbatively calculable matching kernels. The first kernel $R_{gg}$ is a Fourier transform at leading order and contains next-to-leading order (NLO) contributions from gluon diagrams. The $R_r$ and $R_{gq}$ kernels contribute at NLO from quark mixing terms, and both are needed to maintain the exact normalization of $\mathcal{M}(\nu=0,z^2)=1$. The quark mixing contribution has been shown to be negligible compared to current statistical and systematic uncertainties in the lattice data [65] [19, 22], as the mixing is suppressed at NLO by the ratio of color factors $C_F/C_A$. Hence, we will ignore this term for this study. Instead, we focus on reducing higher-twist terms and discretization effects by introducing the matching relation as

The terms $z^2B_1(\nu)$ and $(1/z)P_1(\nu)$ are designed to control systematic higher-twist and discretization effects [36, 37], respectively, and are parametrized by

where $\sigma_{0,n}$ are weighted cosine transforms of Jacobi polynomials, shifted to the domain $0\leq x\leq1$,

In this work, we retain only the first two basis functions, $N=2$, and impose $\sigma_{0,n}(0)=0$ to ensure that the systematic corrections preserve the normalization of the reduced pseudo-Ioffe-time distribution (RpITD).

We construct the gluon PDF using the parameterization

where $\langle x\rangle_g$ is the gluon momentum fraction, $B(\alpha+1,\beta+1)$ is the Euler beta function, and the parameters $\alpha$ and $\beta$ control the small- and large-$x$ behavior of the distribution, respectively. This parameterization determines the leading-twist contribution RpITD, which is fitted simultaneously with the higher-twist and discretization corrections defined in Eq. (17). The resulting model depends on the six fit parameters $\theta=\{\alpha,\beta,b_1,b_1,p_1,p_2\}$, which are determined through a correlated fit to the RpITD data. Specifically, we minimize the correlated chi-squared,

where $\mathcal{M}(\nu,z^2)$ and $\mathcal{M}(\nu,z^2;\theta)$ denote the vectors of measured and modeled RpITD values, respectively. The covariance matrix $C$ encodes correlations among the data points and is used to account for these correlations in determining the best-fit parameters and their uncertainties. We repeat the fitting procedure for the bootstrap samples to estimate the uncertainties and correlations of the extracted parameters.

To investigate the constraints imposed by the correlated fits, we examine the joint distributions of the fitted parameters across the bootstrap samples. Figure 4 presents the resulting parameter distributions and pairwise correlations for the light and strange nucleons. The diagonal panels show the marginal bootstrap distributions of the six fit parameters, while the off-diagonal panels display their pairwise relationships. For clarity, the lower triangle shows correlations for the light nucleon (blue points), and the upper triangle shows those for the strange nucleon (orange points). These distributions provide insight into the uncertainties of the extracted parameters and the correlations between the PDF shape and the systematic corrections. The marginal distributions exhibit broadly similar behavior across the two pion masses, suggesting that the fit parameters are constrained in a qualitatively consistent manner for both ensembles. The pairwise distributions reveal a pronounced correlation between $\alpha$ and $\beta$ which reflects the ill-conditioned nature of the inverse problem, whereby different combinations of endpoint exponents can produce similar predictions for the lattice data, leaving the individual parameters less constrained. An even stronger correlation is observed between $p_1$ and $p_2$, which govern the discretization correction. Since both coefficients enter through a shared $\nu$-dependent contribution proportional to $a/|z|$, their effects can partially compensate for one another over the fitted kinematic range. These correlations illustrate how the lattice data constrain combinations of fit parameters rather than necessarily determining each parameter independently. Accounting for these dependencies is therefore essential when propagating the bootstrap uncertainties to the reconstructed RpITDs and gluon PDF. Overall, we expect better control of the PDF parameters to come as the precision for large-$\nu$ data is improved. The $p_i$ parameters could be better constrained by the inclusion of multiple lattice spacings to uniquely separate $a/z$ discretization effects from the $\nu$-dependence of $P_1(\nu)$.

This comparison allows us to assess how well the fitted model describes the lattice data across the available kinematic range. For the lighter-pion ensemble ($M_\pi\approx310$ MeV), the reconstructed predictions generally describe the lattice data well within their $1\sigma$ uncertainties, particularly at larger momenta $P_z$. However, at smaller momenta, the lattice RpITDs exhibit a more rapid decay than the reconstructed predictions. Similar behavior has been reported in previous gluon PDF studies [68] [21, 24], where the lowest-momentum data were found to deviate from the expected matching behavior. The overall fit quality remains reasonable, with mean reduced chi-squared values of $\chi^2/\mathrm{d.o.f.} = 1.128(228)$ and $1.067(206)$ for the strange and light nucleons, respectively. Thus, the current statistical precision is insufficient to determine whether this difference reflects a systematic deviation from the fitted model. The heavier-pion ensemble ($M_\pi\approx690$ MeV) shows broadly similar agreement, but the reconstructed bands tend to lie above the lattice data at larger separations and smaller momenta. Overall, the reconstructed RpITDs provide a comparison between the fitted model and the lattice data, illustrating both the general level of agreement and the variations observed across the available momenta and spatial separations.

Figure 5 compares the lattice RpITD data with the reconstructed model predictions as functions of the Ioffe time for the two pion-mass ensembles, with different colors denoting the spatial separations $z$. The reconstructed predictions incorporate both the leading-twist gluon PDF and the higher-twist and discretization corrections, with the shaded bands representing the uncertainties propagated from the bootstrap distributions of the fitted parameters. This comparison allows us to assess how well the fitted model describes the lattice data across the available kinematic range.

Finally, we present the extracted unpolarized nucleon gluon PDFs for the two valence pion-masses in left panel of Figure 6. Comparing the PDFs at $M_\pi\approx310$ and $690$ MeV, we observe statistically consistent distributions across $x$, as expected from the extracted parameters described earlier. Minimal pion mass dependence in the PDFs are consistent with observations from earlier pseudo-PDF studies [69] [19, 22] and recent LaMET extractions [70] [24, 25, 27]. Nonetheless, lighter pion masses, should be explored to ensure that this trend is consistent closer to the physical pion mass. In the right panel of Figure 6, we show the lightest pion mass result compared to select phenomenological results: CJ22 [66], JAM24 [67], and CT18 [9]. Our results agree with different global fit behaviors at various ranges of $x$. At $x \le0.2$, the phenomenological results sit right at the $1\sigma$ boundary of the lattice results, but we do not expect the RpITD to have much sensitivity the PDF in this range [71]. In particular, our extraction supports the highly suppressed gluon distribution at $x \ge0.5$ of the JAM24 and CJ22 global fits, while presenting tension of greater than $1\sigma$ with the CT18 extraction in the range $0.3 \lesssim x \lesssim0.6$. Appreciating that there is still much flexibility in the correlated values of the $\alpha$ and $\beta$ parameters presented in Figure 4, it seems likely that inclusion of the RpITD data from this analysis into global fits could produce PDFs with good agreement both the experimental and lattice data, as was seen from the pion PDF study presented in Ref. [72].

Conclusion and Outlook

We have presented the most systematically controlled extraction of the nucleon gluon PDF using the pseudo-PDF methodology to date. Our analysis emphasizes the use of a finer lattice spacing than in previous studies, improved control of excited-state contamination through the recently developed Lanczos-based matrix element analysis [62], and explicit parameterizations of systematic effects in the PDF extraction [36, 37]. We compared a rigorous implementation of our previous two-state fit approach with the Lanczos method, finding reasonable agreement between the two. The Lanczos method, however, exhibits substantially reduced sensitivity to analysis choices, particularly the selection of fit windows. The fine lattice spacing of $a \approx0.06$ fm enabled the use of six discrete separation distances within the perturbative region $z \lesssim0.3$–$0.4$ fm, together with ten nucleon momenta reaching approximately 4.4 GeV. These data provide constraints on the gluon PDF at valence pion masses $M_\pi\approx310$ and 690 MeV.

We explicitly parameterized higher-twist contributions and discretization effects and examined correlations among the resulting fit parameters, finding that the data constrain certain correlated combinations of PDF parameters. We observe minimal dependence on the valence pion mass and identify nonzero discretization effects that can be accounted for within our fitting framework. Further control over these effects could be achieved through simultaneous fits to data from multiple lattice spacings. The resulting range of PDF parameterizations describes the reduced RpITD data well and exhibits reasonable agreement with phenomenological extractions over various regions of $x$. In particular, we find a strongly suppressed gluon PDF at large $x$.

The broad range of PDF parameterizations compatible with our lattice data, together with the improved control of systematic effects, suggests that these RpITD results could provide useful constraints in global PDF analyses. In particular, a simultaneous fit to the lattice RpITD and experimental data could test their mutual consistency and further constrain the gluon PDF, analogous to the joint analysis performed for pion PDFs in Ref. [72].

Acknowledgments

We thank the MILC Collaboration for sharing the lattices used to perform this study. The LQCD calculations were performed using the Chroma software suite [73]. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 through ERCAP; facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy, and supported in part by Michigan State University through computational resources provided by the Institute for Cyber-Enabled Research (iCER). The work of AN is partially supported by U.S. Department of Energy, Office of Science, under grant DE-SC0024053 “High Energy Physics Computing Traineeship for Lattice Gauge Theory”. KE is partly supported by the Early-Start Fellowship from the College of Natural Science at Michigan State University. The work of WG and HL is partially supported by the US National Science Foundation under grant PHY 2209424 and 2514533.

References