Gluon Parton Distribution of the Nucleon from 2+1+1-Flavor Lattice QCD in the Physical-Continuum Limit

We present the first physical-continuum limit $x$-dependent nucleon gluon distribution from lattice QCD using the pseudo-PDF approach, on lattice ensembles with $2+1+1$ flavors of highly improved staggered quarks (HISQ), generated by MILC Collaboration. We use clover fermions for the valence action on three lattice spacings $a \approx 0.9$, 0.12 and 0.15~fm and three pion masses $M_\pi \approx 220$, 310 and 690~MeV, with nucleon two-point measurements numbering up to $O(10^6)$ and nucleon boost momenta up to 3~GeV. We study the lattice-spacing and pion-mass dependence of the reduced pseudo-ITD matrix elements obtained from the lattice calculation, then extrapolate them to the continuum-physical limit before extracting $xg(x)/\langle x \rangle_g$. We use the gluon momentum fraction $\langle x \rangle_g$ calculated from the same ensembles to determine the nucleon gluon unpolarized PDF $xg(x)$ for the first time entirely through lattice-QCD simulation. We compare our results with previous single-ensemble lattice calculations, as well as selected global fits.


I. INTRODUCTION
Many precision phenomenology and theoretical predictions for hadron colliders rely on accurate estimates of the uncertainty in Standard-Model (SM) predictions.Among these predictions, the parton distribution functions (PDFs), the nonperturbative functions quantifying probabilities for finding quarks and gluons in hadrons with particular momentum fraction, are particularly important inputs in high-energy scattering [1][2][3][4][5][6][7][8][9][10][11].The gluon PDF g(x) needs to be known precisely to calculate the cross section for these processes in pp collisions, such as the cross section for Higgs-boson production and jet production at the Large Hadron Collider (LHC) [12,13], and direct J/ψ photoproduction at Jefferson Lab [14].The future U.S.-based Electron-Ion Collider (EIC) [15], planned to be built at Brookhaven National Lab, will further our knowledge of the gluon PDF [16][17][18].In Asia, the Electron-Ion Collider in China (EicC) [19] is also planned to impact the gluon and sea-quark distributions.Although significant efforts to extract the gluon distribution g(x) have been made in the last decade, there are still problems in obtaining a precise g(x) in the large-x.
Recently, progress has been made in the mostcalculated isovector quark distribution of nucleon by MSULat [49], ETMC [51] and HadStruc Collaborations [103], who studied lattice-spacing dependence.MSULat studied three lattice spacings (0.09, 0.12 and 0.15 fm) and pion masses (135,220,310 MeV) and performed a simultaneous continuum-physical extrapolation using a third-order z-expansion on renormalized LaMET matrix elements [49] with nucleon boost momenta around 2.2 and 2.6 GeV.ETMC also uses three lattice spacings, 0.06, 0.08, and 0.09 fm, but with heavier pion mass (370 MeV) and investigated the continuum extrapolation of the data on renormalized LaMET matrix elements with boost momentum around 1.8 GeV [51].Had-Struc Collaboration studied three lattice spacings, 0.048, 0.065, and 0.075 fm with two-flavor 440-MeV lattice ensembles using the continuum pseudo-Ioffe-time distribution (ITD) [103].Most of the works above found mild nonzero dependence on lattice spacing (varying with the Wilson-link displacement) in the nucleon case for LaMET or pseudo-ITD matrix elements.
In contrast with the quark PDFs, the gluon PDFs calculations are less calculated, due to their notoriously noisier matrix elements on the lattice.To date, there have only been a few exploratory gluon-PDF calculations for unpolarized nucleon [36,92,95], pion [96] and kaon [104], and polarized nucleon [97] using the pseudo-PDF [105] and quasi-PDF [38,106] methods.Most of these calculations, like many exploratory lattice calcula-tions, are done only using one lattice spacing at heavy pion mass.
In this work, we report the first continuum-limit unpolarized gluon PDF of nucleon study using three lattice spacings: 0.09, 0.12 and 0.15 fm with pion mass ranging from 220 to 700 MeV using the pseudo-PDF method.The remainder of this paper is organized as follows: In Sec.II, we present the procedure for how the lattice correlators are calculated and analyzed to extract the ground-state matrix elements for the pseudo-PDF method.We then study the lattice-spacing dependence of matrix elements at 310 and 700-MeV pion mass in Sec.III, checking both O(a) and O(a 2 ) forms using multiple continuum-extrapolation strategies.We perform a physical-continuum extrapolation to obtain continuum reduced pseudo-ITDs (RpITDs) matrix element, before final determination of the nucleon unpolarized gluon PDF xg(x) is obtained from the xg(x)/ x g and x g results.Using the gluon momentum fraction calculated on the same ensemble, we obtain the gluon PDF and compare with the phenomenological global-fit PDF results.We consider the quark-mixing systematics, but they are found to be small.The final conclusion and future outlook can be found in Sec.IV.
On each lattice configuration, we calculate the nucleon two-point correlators using multiple sources: with the nucleon interpolation operator χ as lmn [u(y) l T iγ 4 γ 2 γ 5 d m (y)]u n (y) (where {l, m, n} are color indices, u(y) and d(y) are quark fields), the projection operator Γ = 1 2 (1 + γ 4 ), t is lattice Euclidean time, and P z is the nucleon boost momentum along the spatial z-direction.We use Gaussian momentum smearing [127] on the quark field to improve the signal for nucleon boost momenta up to 3.0 GeV.Hundreds of thousands of measurements are made, varying for different ensembles.Compared to our previous nucleon gluon PDF calculation on one a12m310 ensemble with 10 5 measurements [92], this study uses more measurements and varies the lattice spacing.We then calculate the three-point gluon correlator by combining the gluon loop with nucleon two-point correlators, where t is the gluon-operator insertion time, t sep is the source-sink time separation.O g (z, t) is the gluon operator introduced in Ref. [105]: where the operator , and z is the Wilson link length.To extract the ground-state matrix element, we use a two-state fit on the two-point correlators and a two-sim fit on the three-point correlators: two-sim fit two-sim where the |A N,i | 2 and E N,i are the ground-state (i = 0) and first excited state (i = 1) amplitude and energy, respectively.
To visualize our fitted matrix-element extraction, we compare to ratios of the three-point to the two-point correlator The left-hand side of Fig. 1 shows example ratios for the gluon matrix elements from all four ensembles at pion masses M π ∈ {220, 310} MeV at selected momenta P z and Wilson-line length z.The left column shows the ratio plots with data points of R from different source-sink separation, t sep , along with the reconstructed bands from the fit, showing how well the fit describing the data in Eq. 6; the final ground-state matrix elements are shown in grey bands.We observe that the ratios increase with increasing source-sink separation t sep and continuously to approach the ground-state matrix elements obtained from the simultaneous two-state fit to three-point correlators with five inputs of t sep .The middle and right two-sim fit two-sim Example ratio plots (left), and two-sim fits (last 2 columns) of the strange nucleon correlators at pion mass Mπ ≈ 700 MeV from the a09m310, a12m220, a12m310 and a15m310 ensembles.The gray band shown on each plot is the extracted ground-state matrix element from the two-sim fit that we use as our best value.From left to right, the columns are: the ratio of the three-point to two-point correlators with the reconstructed fit bands from the two-sim fit using the final tsep inputs, shown as functions of t − tsep/2, the one-state fit results for the three-point correlators at different tsep values, the two-sim fit results using tsep ∈ [t min sep , t max sep ] varying t min sep and t max sep .
columns of Fig. 1 show how the ground-state matrix elements vary with t min sep and t max sep with the same t max sep and t min sep used in the left column, respectively.The grey bands in each plot of the middle and right columns are the same ground-state matrix elements extracted from the fit shown in the corresponding plots in left column; these demonstrate how stable our ground-state matrix element extractions are.We observe that overall the ground-state matrix elements are consistent with each other within one standard deviation with different t min sep choices.There seems to be some hint of t max sep -dependence in the groundstate matrix elements, but they do converge.Taking the a12m220 ensemble as an example, we observe larger fluctuations in the matrix-element extractions when small t min sep = 3, or small t max sep = 6 and 7, are used.The groundstate matrix element extracted from two-sim fits comes into reasonable agreement when t min sep > 5 and t max sep > 8. Similar results are observed with the nucleon when calculated at strange-quark mass, as shown in Fig. 2.

A. Lattice-Spacing Dependence of RpITDs
Using the nucleon ground-state matrix elements, we can now compute the reduced Ioffe-time pseudodistribution (RpITD) [81,82,106,128] where Ioffe time ν = zP z , and M(ν, z 2 ) are the nucleon matrix elements at boost momentum P z and gluon operators with Wilson displacement z.By construction, the renormalization of O(z) and kinematic factors are cancelled in the RpITDs, and the ultraviolet divergences are removed.The RpITD double ratios employed here are normalized to one at z = 0, and the lattice systematics are reduced due to the double ratio.These RpITDs will be input into the pseudo-PDF framework detailed in Ref. [81] to obtain the unpolarized nucleon gluon PDFs.
We first examine the pion-mass and lattice-spacing dependence of the nucleon gluon RpITDs.The top panel of Fig. 3 shows the RpITDs at boost momentum around 1.3 GeV as functions of Ioffe time ν for the a12m220, a09m310, a12m310, and a15m310 ensembles.Note that given the different size of the lattice dimensions and lattice spacings, there is no easy way to use the same boost momentum for all four ensembles.For the a12m310, a12m220, and a09m310 ensembles, we are able to find close momenta, but the boost momentum on a15m310 is not so close, 1.54 GeV.Nevertheless, this allows us to study the lattice-spacing and pion-mass dependence without any interpolation; we find in both cases mild dependence.The a = 0.12 fm ensembles seem to prefer slightly higher central values for the RpITDs, but they are consistent with those from a09m310 and a15m310 within one standard deviation.The RpITDs at boost momenta Pz ≈ 1.3 GeV (top) and 2 GeV (bottom) as functions of Ioffe-time ν obtained from the fitted bare ground-state matrix elements for Mπ ≈ {220, 310, 310, 310} MeV on a12m220, a09m310, a12m310, a15m310 ensembles, respectively.There is no visible lattice-spacing or pion-mass dependence.
We then choose a fixed pion mass and Ioffe time ν to explore the lattice-spacing dependence in more detail.We examine the RpITDs at three values of z and set n z ∈ {4, 3, 2} in lattice units (P z = 2π L n z ) for the a09m310, a12m310 and a15m310 ensembles, respectively, to check the lattice-spacing dependence, where the Ioffe time ν = zP z are the same for all ensembles.We assume a linear O(a) or O(a 2 ) dependence and fit the RpITD data: for n = 1 and 2, respectively.The fit results are shown in Fig. 4 for both light and strange nucleons at ν = {π/4, π/2, π}.All the plots show the RpITD points at a fixed ν are consistent within one-sigma error; however, there seems to be a trend that the latticespacing dependence is stronger at larger ν, making the continuum-extrapolated results larger.The χ 2 /dof are all within one standard deviation of 0.5, showing that there is no preference in choice of ν in this respect.In both O(a) and O(a 2 ) extrapolation, the light and strange nucleons appear to have opposite trends with lattice spacing; however, the slopes of the fits are consistent with zero within one standard deviation for each pion mass over all the selected ν.Corresponding to a reading of Fig. 4 29), and 0.837 (84).The deviation of the continuum-limit RpITD increases as ν increases, with the continuum-limit error at ν = π consistently being about a factor of ten larger than those at ν = π/4.The continuum-limit RpITDs using a and a 2 extrapolation are consistent with each other, but the former has larger error, due to extrapolating over a larger distance.The error in the light-nucleon extrapolated RpITD values is consistently double that of the strange nucleon.

B. Continuum-Physical Extrapolation
We extrapolate the nucleon gluon PDF at the physical pion mass and continuum limit based on the RpITDs for the a12m220, a09m310, a12m310 and a15m310 ensembles.Given that we observe small differences between a and a 2 extrapolation at fixed M π and ν in the previous subsection, we focus on a 2 extrapolation here.We use the following ansatz, linear in pion mass and lattice-spacing squared, to extrapolate to the continuum-physical limit: Given the noisiness of the gluon RpITD data, k max = 2 is used in this study.We found that in all the ensembles in our calculation, the z 2 -dependence, c z (a, M π ) is consistent with zero within two standard deviations.The example reconstructed fitted bands for a09m310 and a12m220 are shown in the top plot of Fig. 5.We drop data points at a ν value if they have errors more than 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.7 0.8 0.9 1.0 1.1 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.7 0.8 0.9 1.0 1.1 twice as large another data point at that ν; this improves the clarity of the diagram showing the extrapolation effects.Our fit ansatz is able to describe our data reasonably with ν up to about 7; in the future, when larger P z is used to reach larger ν, a better interpolation form will likely be needed to describe the small-and large-ν regions well.We show the continuum-physical RpITD band on the top plot of Fig. 5 with all the data points from the four ensembles and a 2 extrapolation to the continuum-physical band.The open symbols indicate the strange-mass nucleon calculation from the ensemble.With the a15m310, a12m310 and a09m310 ensembles, since we have the same number of measurements for both strange and light nucleons, within each ensemble we bootstrap the light and strange renormalized matrix elements to keep the correlations.Across the ensembles, the data are independent and the typical bootstrap treatment is used.For comparison purposes, we also replace the a 2 term in Eq. 9 with a; its physical-continuum results are shown as the band with dashed center line in Fig. 5. Similar to what we observe in Fig. 4, both extrapolations give us a consistent continuum-physical limit (within one standard deviation) with Ioffe time ν up to 7, but the O(a)-extrapolated continuum-physical RpITD has larger errors, especially in the larger-ν region.

C. Gluon PDF Results
With the physical-continuum RpITD obtained in the previous section, we can now extract the gluon PDF distribution using the pseudo-PDF matching condition [105] that connects the RpITD M to the lightcone gluon PDF g(x, µ 2 ): where µ is the renormalization scale in the MS scheme and x g = 1 0 dx xg(x, µ 2 ) is the gluon momentum fraction of the nucleon.R gg is the gluon-in-gluon matching kernel originally derived in Ref. [105] and has been used in previous gluon PDF lattice works [92,95,96,104].We use the RpITD extrapolated using the a 2 term (in Eq. 12) at physical pion mass and continuum limit with a fit range of ν ∈ [0, 7], corresponding to the region where we have data in all ensembles.We ignore the quark-PDF contribution to the RpITDs in this calculation; it is likely to be small, based on our past study of the pion gluon PDF [96].We will later estimate the quark contribution as systematic effect.One can obtain the gluon PDF g(x, µ 2 ) by fitting the RpITD through the matching condition in Eq. 10.We adopt the phenomenologically motivated form commonly used in the global analysis for x ∈ [0, 1] and zero elsewhere.The beta function C is used to normalize the area to unity.Such a form is also used in global fits to obtain the global PDF fits, such as nucleon gluon PDF by CT18 [7] and the nucleon and pion gluon PDF by JAM [129][130][131].We then fit the lattice physical-continuum RpITDs M lat (ν, z 2 , a, M π ) obtained in Eq. 9 to the parametrization form M fit (ν, µ, z 2 , a, M π ) Examples of the RpITDs M reconstructed bands from fits in Eq. 9 for a09m310 (blue points and light blue band), a12m220 (green) lattice ensembles.The fit ansatz is able to describe the data well.(Bottom) Collected data for all ensembles with a (dashed band) and a 2 (solid band) continuum extrapolation at the physical pion mass.Open symbols indicates the data point from the same-symbol ensemble but at the heavier quark mass.
in Eq. 10 by minimizing the χ 2 function, Our results for the continuum-physical unpolarized gluon PDF xg(x, µ)/ x g are shown in Fig. 6, along with the same determination from the smallest lattice-spacing ensemble obtained in this work, and selected global-fit gluon PDFs from CT18 [7] and NNPDF3.1 [6] NNLO analysis.The gluon distribution in continuum-physical limit has much larger errors by a factor of 3-5 than those obtained from single-lattice-spacing analysis, due to the continuum extrapolation.Overall, the results from single-ensemble calculations on a09m310 are consistent with the continuum-physical one (which has much larger uncertainties).To reduce the errors in the continuumphysical distribution will be difficult, since it requires re-duced errors in all ensembles, increasing the calculation cost by at least another order of magnitude.Both of our lattice distributions agree with the global-fit gluon distribution at mid to large x but deviate for x < 0.3.This is likely due to lack of large-ν lattice data in the input, which has higher sensitivity to the smaller-x data.Future calculations to push for even larger P z will be needed to improve the small-x gluon distribution.The unpolarized gluon PDF, xg(x, µ)/ x g as a function of x and its zoomed in plot, obtained from the fits to the smallest-lattice-spacing ensemble data compared with the fit to the data obtained from extrapolation to physical pion mass and continuum limit.The black solid line is the central value of the fit to the continuum-physical PDFs, including the gluon-in-quark term in the matching, using CT18 for the quark PDF contributions.The results from the global fits by CT18 [7] and NNPDF3.1 [6] NNLO gluon PDFs are also shown in the plots, and our gluon PDF results are consistent with the global fits for x ∈ [0.3,1] We now consider the systematic uncertainty coming from neglecting the contribution of the quark term, x g R gq (xν, z 2 µ 2 ) in Eq. 10.We ignored this contribution initially based on the assumption (motivated by global fits) that the nucleon total quark PDF q S (x) is smaller than the gluon PDF.We can estimate the systematic due to omitting the q S (x) contribution by using the nucleon flavor-dependent quark PDFs from CT18 at NNLO [7].Following a similar procedure to Ref. [96], we add the quark-gluon mixing term to the extraction of xg(x)/ x g from the RpITD in Eq. 10.The central value of the updated xg(x)/ x g including both gluon-in-gluon (gg) and gluon-in-quark (gq) contributions is shown in Fig. 6 (black solid line).The difference is much smaller than the current statistical errors, so we will ignore this quark contribution in this calculation.However, in the future, when the continuum-physical gluon PDF precision are improved, we should re-examine this contribution more carefully.xg(x, µ)/ x g at µ 2 = 4 GeV 2 as function of x, extracted from continuum-physical RpITDs using the two-(green) and three-parameter (orange) fits described in Eqs.11 and 13, respectively.The results are consistent within statistical errors.
We investigate the systematic uncertainty introduced by the choice of parametrization form used for f g (x, µ).We consider a three-parameter form used in PDF global analysis and some lattice calculations, .
(13) We fit our continuum-physical-limit RpITDs to this form up to maximal Ioffe time ν max = 7.A comparison of the fit-form choice is shown in Fig. 7.We find the goodnessof-fit improves slightly due to the introduction of a new free parameter in the fit form, but the gluon PDF results are noisier and are consistent with the two-parameter fit.We will use xg(x, µ)/ x g from the two-parameter fit as our main result for this work.
constant [133,134].The gluon momentum fraction was extrapolated to the continuum-physical limit and found to be consistent with other recent lattice-QCD results at physical pion mass.Our final unpolarized nucleon gluon PDF xg(x) extrapolated to physical pion mass M π = 135 MeV and the continuum limit a → 0 is shown as green bands in Fig. 8; once again, we found reasonable agreement with the global fits from CT18 [7] and NNPDF3.1 [6] NNLO analysis for x ∈ [0.25, 1], even thought the the gluon momentum fraction obtained from the global fits is about two-sigmas lower than the lattice calculations.We do observe tension with gluon PDF from JAM20 [129] analysis for x < 0.6 regions but its gluon PDF also behave quite different from the the CT18 and NNPDF results, even with smaller errors; we look forward to updates on the global-fit community on resolving these discrepancy.
We also compare our results with other previous lattice-QCD calculation on xg(x).The light cyan bands in Fig. 8 shown the first pseudo-PDF calculation done using clover-on-HISQ with 0.12-fm lattice spacing and 310-and 700-MeV pion mass using 898 lattice configuration with 32 sources per configuration for nucleon two-point correlators [92].The results are extrapolated to physical pion mass using naive two valence pion mass extrapolation with xg(x) reconstructed by multiplying the gluon momentum fraction taken from Ref. [55].The blue bands in Fig. 8 show a followup calculation performed by HadStruc collaboration using 2+1 dynamical flavors of clover fermions with stout-link smearing on the gauge fields, 0.09-fm lattice spacing, 358-MeV pion mass, and 64 source measurements on 349 lattice configurations with gradient-flow improved gluonic operators [95].They used multiple nucleon interpolating fields, allowing them to use generalized eigenvalue method to determine the best overlap with ground-state nucleon gluonic matrix elements.They used the gluon momentum fraction obtained from an independent lattice work (2+1+1-flavor at physical pion mass) to determine xg(x).The outer blue bands indicate their uncertainty estimated from x g .We also show our result on the ensemble with 0.09-fm lattice-spacing and 310-MeV pion mass as a purple band in Fig. 8; it used about 300k measurements spread out over 1000 lattice configurations.Our single-ensemble results have errors comparable to (in some regions, smaller than) CT18 and NNPDF.The lattice-spacing and pion-mass here is similar to those used in the HadStruc calculation [95] but without the additional uncertainties due to continuumphysical extrapolation (shown as a green band).There are noticeable deviations from the HadStruc results, especially in the larger-x region; their large-x gluon PDF is much smaller than ours.However, given that multiple methodological aspects are done quite differently (for example, we used the momentum fraction from the same lattice ensemble and different gluon-operator smearing), it may require the full calculation, including continuumphysical extrapolation, to have meaningfully compare them.All the prior single-ensemble lattice results (without the systematics from lattice discretization) agree with our continuum-physical xg(x) due to the larger total errors from the continuum-physical extrapolation.Future work to include finer lattice-spacing and 220-MeV or lighter pion masses in the extrapolation will help to improve the continuum-physical determination of the lattice gluon PDF.

IV. SUMMARY AND OUTLOOK
We extracted the nucleon x-dependent gluon PDFs xg(x) using clover fermions as valence action and 310-MeV 2+1+1 HISQ configurations generated by the MILC Collaboration at three pion masses and three lattice spacings.We found their dependence to be weak at the current statistics of hundreds of thousands of nucleon twopoint correlator measurements.We removed the excitedstate contributions to the ground-state matrix elements using a two-state fitting strategy and studied the stability of the extraction of the ground-state matrix elements with various fit ranges.We then calculated the reduced pseudo-ITD using the fitted matrix elements and studied their pion-mass and lattice-spacing dependence, which are also mild.We then extrapolated the reduced pseudo-ITDs to physical-continuum limit before extracting the gluon parton distribution xg(x)/ x g in the MS scheme at 2 GeV.Using the nonperturbatively renormalized nucleon momentum fraction calculated on clover-on-HISQ ensembles, we were able to compare our single-ensemble xg(x) calculations at lattice spacing 0.09 fm with prior lattice calculations.We found both our 0.09-fm and continuum-physical limit xg(x) to be in good agreement with CT18 and NNPDF3.1 NNLO global-fit results in the range x ∈ [0.2, 1] in MS scheme at 2 GeV.We have used the CT18 quark PDFs to estimate the quark-gluon mixing and checked the gluon PDF fit-form dependence, but these errors are not included in our continuum-physical results, since the statistical errors are large in comparison.Future work should include finer lattice spacings, even higher statistics to improve the signal and larger boost momentum to expand the range of ν, which will improve the reliability of the results at small x.

FIG. 1 .
FIG. 1.Example ratio plots (left) and two-sim fits (last 2 columns) of the light nucleon correlators at pion masses Mπ ≈ {310, 220, 310, 310} MeV from the a09m310, a12m220, a12m310 and a15m310 ensembles.The gray band shown on each plot is the extracted ground-state matrix element from the two-sim fit that we use as our best value.From left to right, the columns are: the ratio of the three-point to two-point correlators with the reconstructed fit bands from the two-sim fit using the final tsep inputs, shown as functions of t − tsep/2, the one-state fit results for the three-point correlators at different tsep values, the two-sim fit results using tsep ∈ [t min sep , t max sep ] varying t min sep and t max sep .
FIG. 5. (Top)Examples of the RpITDs M reconstructed bands from fits in Eq. 9 for a09m310 (blue points and light blue band), a12m220 (green) lattice ensembles.The fit ansatz is able to describe the data well.(Bottom) Collected data for all ensembles with a (dashed band) and a 2 (solid band) continuum extrapolation at the physical pion mass.Open symbols indicates the data point from the same-symbol ensemble but at the heavier quark mass.
FIG. 6.The unpolarized gluon PDF, xg(x, µ)/ x g as a function of x and its zoomed in plot, obtained from the fits to the smallest-lattice-spacing ensemble data compared with the fit to the data obtained from extrapolation to physical pion mass and continuum limit.The black solid line is the central value of the fit to the continuum-physical PDFs, including the gluon-in-quark term in the matching, using CT18 for the quark PDF contributions.The results from the global fits by CT18[7] and NNPDF3.1[6]NNLO gluon PDFs are also shown in the plots, and our gluon PDF results are consistent with the global fits for x ∈ [0.3,1] FIG. 7.xg(x, µ)/ x g at µ 2 = 4 GeV 2 as function of x, extracted from continuum-physical RpITDs using the two-(green) and three-parameter (orange) fits described in Eqs.11 and 13, respectively.The results are consistent within statistical errors.