Nonperturbatively-renormalized glue momentum fraction at physical pion mass from Lattice QCD

We present the first nonperturbatively-renormalized determination of the glue momentum fraction $\langle x \rangle_g$ in the nucleon, based on lattice-QCD simulations at physical pion mass using the cluster-decomposition error reduction (CDER) technique. We provide the first practical strategy to renormalize the glue energy-momentum tensor (EMT) nonperturbatively in the RI/MOM scheme, and convert the results to the $\overline{\textrm{MS}}$ scheme with 1-loop matching. The simulation results show that the CDER technique can reduce the statistical uncertainty of its renormalization constant by a factor of ${\cal O}$(300) in calculations using typical state-of-the-art lattice volume, and the nonperturbatively-renormalized $\langle x \rangle_g$ is shown to be independent of the lattice definitions of the glue EMT up to discretization errors. We determine the renormalized $\langle x \rangle_g^{\overline{\textrm{MS}}}(2\textrm{ GeV})$ to be 0.47(4)(11) at physical pion mass, which is consistent with the experimentally-determined value.

Introduction: A longstanding problem raised by deepinelastic scattering (DIS) and Drell-Yan experiments on the nucleon is that the gluons contribute almost as large a fraction of the nucleon momentum as the quarks [1,2], contradicting the naive quark model. The momentum fractions of the quarks and glue equal the second moments of their respective parton distribution functions (PDFs) f p (x) (p = u,ū, d,d, s, ..., g): where the PDF can be determined from global fits of experimental results with certain assumptions about their functional forms. The recent CT14NNLO global PDF fit [2] yields x MS g (2 GeV) = 0.42 (2), and the value at the TeV scale will be around 0.5 which is irrespective of its value at lower scales. Beside the importance in understanding the nucleon momentum, the value of x g is also an important input to obtain the glue contributions to the nucleon mass and spin [3,4], so calculating it from a first-principle lattice-QCD simulation is of fundamental interest, in addition to providing an independent input and check of the experimental PDF determinations.
Lattice calculations of x g in the nucleon [4][5][6][7] have been significantly refined in the last 10 years. However, values of x MS g (2 GeV) vary widely, two quenched cal-culations found 0.43 (9) and 0.33 (6) [4,5], and recent dynamical N f = 2 calculation obtained 0.267 (22)(30) [6,7]. The recent quenched (Refs. [4,5]) and dynamical (Refs. [6,7]) lattice calculation of x g used different lattice definitions of the gluon energy-momentum tensor (EMT) with the 1-loop renormalization based on the lattice perturbation theory (LPT). It is known that LPT is poorly convergent at 1-loop level without smearing of the glue EMT [8,9], and LPT calculations beyond 1loop level are extremely difficult. Whether smearing of the glue EMT can improve the convergence of LPT remains an open question, but it was found in Ref. [10] that hypercubic (HYP) smearing [11] of the glue operator can change the bare glue matrix element by a factor of ∼3. Nonperturbative renormalization (NPR) of x g is thus essential to check whether different lattice definitions of the glue EMT and smearing can provide a consistent prediction of x g .
In this work, we present the first NPR of the glue EMT using the cluster-decomposition error reduction (CDER). We confirm that, nonperturbatively-renormalized x g is independent of the lattice definition of the glue EMT and whether the HYP smearing is applied to it.
The glue NPR technique that we introduce will be applicable for the quantities beyond the x g . State-of-theart calculations of the glue spin contribution to the proton spin [10] and the glue transversity in hadrons [12] have been presented recently, the renormalization of the glue operators in these calculations are determined at the 1-loop level or neglected entirely. Approaches that target the entire glue PDF instead of the moments, like largemomentum effective theory (LaMET) [13] and the lattice cross-section approach [14], have been explored recently. NPR will be also essential to obtain accurate predictions for those quantities.
NPR Simulation Strategy: At tree level, the glue EMT T g,µν ≡ F µρ F ρ ν − 1 4 g µν F 2 includes 9 Lorentz structures, where µ and ν denote the external Lorentz indices of the EMT, ρ (or τ ) is the Lorentz index of the external gluon state A ρ/τ . But under the physical conditions p ρ = p τ = 0, p 2 = 0 [15], only a single structure survives, The lattice calculation of the gluon propagator is generally off-shell, so the requirement p 2 = 0 cannot be satisfied directly; therefore, some restrictions on the external momentum and indices must be imposed to pick up the coefficient of the Lorentz structure in Eq. (2). More precisely, the RI/MOM renormalization constant of the off-diagonal pieces of the glue EMT at the renormalization scale µ 2 R = p 2 can be defined using the following approach, which is analogous to that commonly used for the quark bilinear operators [16]: where the index ρ is not summed and V is the volume of the lattice. The regularization independent (RI) gluonfield renormalization constant Z RI g is replaced by the gluon propagator on the right-hand-side of the second equal sign and then does not appear in the calculation directly, so it will not be discussed in details here.
The Landau gauge-fixed gluon field A ρ (p) used above is defined from the gauge links U µ (x) as: Note that even though the operator T may be HYP smeared, no smearing will be applied to the gauge field A ρ (p), since the gauge action is not smeared and no reweighting is applied to the configurations. Similarly, the RI/MOM renormalization constants of the traceless diagonal pieces of the glue EMT can be defined by: The bare lattice glue EMT can be defined by the clover definition of the field tensor F µν [4,5], where the plaquette P µ, and P [µ,ν] ≡ P µ,ν − P ν,µ . The bare traceless diagonal component T g,µµ also has a simpler definition (the plaquette definition) [6,7]: Different definitions and choices of smearing on the links U µ (x) in these definitions of T g yield different bare hadron matrix elements, but the renormalized results should agree up to O(a 2 ) correction. After the renormalization constant Z −1 (µ 2 R ) is obtained perturbatively or nonperturbatively under the lattice regularization at µ 2 R = p 2 , the matching factor to convert the result to the MS scheme should be calculated using dimensional regularization. At the µ R used in this work, the 1-loop corrections to match the MS scheme at 2 GeV are at a few percent level [17]. The mixing with the quark EMT is also small [17] and will be considered as a systematic uncertainty; more detailed discussions of the matching and mixing effects can also be found there.
Calculation of the correlation function is numerically challenging, even when the gluon propagator has been determined at better than the 1% level.
The glue operator renormalization constants Z −1 T in MS at 2 GeV with and without CDER (i.e., cutoffs on the distance between the gauge fields/operator). Without CDER the errors are large and the signal cannot be resolved (bands in the background). The errors can be reduced by a factor of ∼300 with r1 = 0.9 fm, r2 = 1.3 fm, shown by the red dots (blue boxes) for Z −1 T with (without) HYP smearing.
volume L 3 × V = 48 3 × 96 (L = 5.5 fm) [18]. The statistical uncertainties are very large and Z −1 T can not be resolved at any scale.
However, we can apply the cluster-decomposition error reduction (CDER) technique to reduce the errors [19]. The cluster-decomposition principle enunciates that correlerators falloff exponentially in the distance between operator insertions, and implies that integrating the correlator over this distance beyond the correlation length will only garner noise not signal. The CDER technique will cut off the volume integral beyond a characteristic length, and then one can gain a factor of √ V in the signal to noise ratio Ref. [19]. Applying CDER to C 3 (p) in Eq. (8) introduces two cutoffs, r 1 between the glue operator and one of the gauge fields, and r 2 between the gauge fields in the gluon propagator, and then leads to the cutoff correlator: For example, with cutoffs r 1 =0.9 fm, r 2 =1.3 fm, the statistical uncertainty can be reduced by a factor of approximately 300. Using these parameters, a very clear signal can be resolved, shown as the red dots and blue boxes in Fig. 1, for Z −1 T with and without HYP smearing, respectively. The values of Z −1 T differ by a factor of ∼3 for the calculations with or without the HYP smearing, at a 2p2 ∼ 1.
A naive cost estimate for the partial triple sum on the volume V in Eq. 9 is O(V r 4 1 r 4 2 ), but the practical cost can be reduced to O(V log V ) by applying the fast Fourier transform several times [19]. Details of the implementation are provided in the supplementary material [20].
Tests on CDER: To test whether Eq. (9) under CDER leads to some systematic bias, we calculated Z −1 T and Z −1 without CDER on 32,254 decorrelated configurations of a quenched Wilson gauge ensemble "24Q" with a = 0.098 fm and L 3 × T = 24 3 × 64 and compared them with those from 300 of the 32,253 configurations with CDER. It is found that CDER results with r 1 ≥ 0.7 fm and r 2 ≥ 1.0 fm agree with the CDER-free results for all a 2p2 . The detailed comparison is given in the supplementary material [20].
The MS renormalization constant Z −1 (2 GeV) on the 24C ensemble as a function of a 2p2 , with different cutoffs on the gluon field-operator correlation (r1) and propagator (r2). A high-statistics calculation without cutoff on a lattice with smaller volume but the same paramters is also presented (red boxes) for comparison.
We also studied the dynamical case: we calculate Z −1 (1 HYP) with CDER on 2,123 configurations of the 2-flavor clover fermion Lüscher-Weisz gauge ensemble "24C" with lattice spacing a = 0.117 fm, m π = 450 MeV and L 3 × V = 24 3 × 64 [21]. For comparison, we repeat the calculation of Z −1 on 19,570 configurations on the "12C" ensemble (with the same lattice setup as 24C except a smaller volume 12 3 × 24) without CDER. Similar R dependence test as in the quenched case is provided in the supplementary material [20]. As in Fig. 2, we should choose r 1 ≥ 0.9 fm and r 2 ≥ 1.3 fm on 24C (black crosses) to get the consistent results with those on 12C without CDER (the red boxes).
If we fit the CDER-result of Z −1 on 24C with a polynomial form including a 2np2n (n ≤2) terms in the range a 2p2 ∈[1.5, 5], the result is 0.89(2) with χ 2 /d.o.f.=0.80. in resolving clean signal of Z −1 T , it is nevertheless important to confirm that the renormalized x g is independent of the lattice definition of T g or whether the HYP smearing is applied, up to O(a 2 ) corrections. Fig. 3 gives the CDER-results on the 48I ensemble as the functions of a 2p2 . The red dots and blue boxes show Z −1 T with and without HYP smearing, respectively, using the clover definition in Eq. (6); the green triangles show the HYPsmeared case using the plaquette definition in Eq. (7), Z −1 T . The a 2p2 dependence and the a 2p2 → 0 limit of the renormalization constants are different between the different definitions, while the presumed rotation symmetry breaking between Z −1 (black triangles) and Z −1 T are consistent with zero within the uncertainties. With the functional form Z −1 T (a 2p2 ) = Z −1 T (0) + C 1 a 2p2 + C 2 a 4p4 , we fit the range a 2p2 ∈ [1. To determine the bare x g , the following ratio is calculated in the rest frame of the nucleon on 81 configurations of the 48I ensemble with partially quenched valence overlap fermion for the pion mass m π ∈ [135, 372] MeV, where χ is the nucleon interpolation field, Γ e is the unpolarized projection operator of the proton and M N is the nucleon mass. When t f is large enough, the derivative of the t-summed ratio R(t f , t) becomes the glue momentum fraction in the nucleon, as applied in the recent high-accuracy nucleon matrix element calculation [22], up to the excited-state contamination at O(e −δm t f ). The calculation setup is the same as for our previous work on the glue spin [10]: a 4 × 4 × 4 smeared grid source with low mode substitution (LMS) [23] is used for the nucleon two-point functions, and all the time slices are looped over to increase statistics. We followed the same strategy in Ref. [19] to apply CDER to the numerator of R(t f , t). With a cutoff around 1 fm which is enough as demonstrated in the the NPR cases studied here, the statistical uncertainties ofR(t f ) can be reduced by an factor of ∼10. The systematic uncertainties in bareR(t f ) due to CDER will be investigated in the future following the strategy in Ref. [19]. The renormalizedR(t f ) at m π = 372 MeV is shown in Fig. 4 as a check with the best signals we have. The errors from Z T and the barē R(t f ) are combined in quadrature. As shown in that figure, even though the renormalization constants with or without HYP smearing differ by a factor of ∼3 as we saw in Fig. 3, the renormalizedR R (t f ) ≡ Z TR (t f ) are consistent within 2σ for t f ≥ 4.

Final result and summary:
We fitR(t f ) to a constant in the range t f > 7a to obtain x g and plot its m 2 π dependence in Fig. 5. With a linear fit to m 2 π for m π < 400 MeV on the 1-HYP smeared data with the clover definition, we obtain x MS g (2 GeV) at the physical pion mass to be 0.42(4) (11). The variance of the values from three definitions, the uncertainties of the renormalization constants, and the mixing effect from the quark momentum fraction x q (which is estimated by 1 − x g times the 1-loop mixing coefficient 0.1528 [17]) are combined in quadrature as the systematic uncertainty. The prediction is consistent with the global fitting result CT14 [2] 0.42 (2) in MS at the same scale. The major systematic uncertainty is the mixing from the quark and can be eliminated with a similar non-perurabtive calculation with the quark external states.
In summary, we have presented a systematic implementation of NPR for the glue momentum fraction x g . We demonstrated that the CDER technique can provide an unbiased improvement on the lattice with the cutoffs r 1 ∼ 0.9 fm and r 2 ∼ 1.3 fm, and that the renormalized x g is insensitive to the lattice definition of the glue EMT or HYP smearing within uncertainties.
Our calculation also shows that HYP smearing can make the a 2p2 dependence of the renormalization constant much stronger than the case without HYP smearing, even though the a 2p2 -extrapolated value can be closer to one. The cases with more steps of HYP smearing is shown in the supplementary material [20]. Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575. We also thank National Energy Research Scientific Computing Center (NERSC) for providing HPC resources that have contributed to the research results reported within this paper. We acknowledge the facilities of the USQCD Collaboration used for this research in part, which are funded by the Office of Science of the U.S. Department of Energy.

SUPPLEMENTARY MATERIAL
A. Implementing CDER with FFT The 3-point glue correlation function can be approximated with applying CDER: where r 1 is the spatial cutoff on the glue propagator, and r 2 is that on distance between the glue operator T µν (z) and the gauge field A p (x). With non-symmetric cutoffs as above, the cost of the computing C 3 can be reduced to V log V using the following strategy: 1. Construct O r µν (x) = |r |<r2 d 4 r T µν (x + r ) by Fourier transforming T µν (x) and f (x) = θ(r 2 − |x|), multiply the transformed functions together in momentum space, and then perform the anti-FT.
3. Apply the cluster decomposition to d 4 xd 4 ye ip·(x−y) B r ρµν (x)A ρ (y) [19]: perform the FT for both A and B, applying the anti-FT to A(p)B(−p), apply the cut g(x) = θ(r 1 − |x|) in coordinate space and then FT the product.
The CDER with symmetric cutoffs can also be efficient if a V log V implementation can be obtained. 6. The cutoff R dependence for r1,2 of the renormalization constant Z −1 (2 GeV) and Z −1 T (2 GeV) on the 24Q ensemble with a 2 p 2 = 2.00. Calculation on 300 configurations with r1 ≥ 0.7 fm and r2 ≥ 1.0 fm are consistent with those using 32,253 configurations without any cutoff. The result is less sensitive to the cutoff r1 than r2; thus, most of the variance reduction comes from reducing r1, while reducing r2 is also useful. The blue and black data are shifted horizontally to enhance legibility.  Since the number of configurations in the 48I ensemble at m π =140 MeV is limited, we turn to three ensembles with smaller volume and larger statistics to check the systematic uncertainties of the CDER approach. To reduce statistical uncertainties then provide a stronger check, we will apply 1 step of HYP smearing on the glue EMTs used in this section.  T results with a 2p2 =2.00, 2.48 and 3.00 respectively. In those figures, the red bands show the results on 32,253 configurations without CDER, and the black boxes show the results with r 1 = 0.7 × r 2 = R agree with the red bands for all the R's not smaller than 0.7 fm. Result with the cutoff on either r 1 or r 2 set to ∞ (the green triangles and purple dots) are also shown in the figures, and it is obvious from the most left data points that the cutoff effects on r 2 are as strong as those on r 1 when r 1 = 0.7 × r 2 . Thus setting the r 1,2 with this relation can be a proper choice to simplify the parameter tuning. The results also demonstrate that cutoffs on either r 1 or r 2 also reduce the statistical uncertainties of Z −1 . The improvement from CDER in this case is only approximately a factor of V 24Q /V 48I × 300 ∼ 13, as the volume is much smaller than that of the 48I ensemble shown in the content.  9. The cutoff R dependence of the renormalization constant Z −1 (2 GeV) on the 24C/12C ensembles with a 2 p 2 = 2.00 and 2.48. Fig. 9 shows similar plots for the dynamical case with 24C and 12C lattices (a = 0.117 fm, m π = 450 MeV, L 3 × V equal to 24 3 × 64 and 12 3 × 24 respectively). The red bands show the results on 19,570 configurations without any cutoff, and the data points show the CDER-results. They are all consistent for all the R's not smaller than 0.9 fm.
C. The discretization error with more steps of HYP smearing In this section, we repeat the NPR and matrix elements calculation on 48I, but with 2 and 5 steps of HYP smearing. As shown in the left panel of Fig. 10, Z −1 T becomes increasing non-linear on a 2p2 when more HYP smearing steps are applied on the glue EMT. Without HYP smearing, the a 2p2 dependence of Z −1 T can be well described by a linear term and the coefficient of the next order a 4p4 term is consistent with zero. With more HYP smearing steps, the coefficients of the a 2p2 and a 4p4 terms increase significantly. Since all momenta p on the external legs of the glue EMT will be integrated in the hadron matrix element, a 2np2n corrections will result in O(a 2n ) discretization errors at finite lattice spacing. From the renormalizedR(t f ) in the right panel of Fig. 10, the results with 2 steps of HYP smearing still agree with the results with 1 step of HYP smearing; but if we jump to the 5-step HYP smearing used by some previous studies, the a 2np2n corrections will be much larger and the renormalized result will be have large systematic uncertainties from determining Z −1 T (green triangles and blue boxes). In the 5HYP case, with the same range a 2p2 ∈ [1.5, 5] and the polynomial form up to a 4p4 term, Z −1 T (5HYP)=0.62(3) is obtained with χ 2 = 1.2 (the default fit). If the a 6p6 term is added and the range is switched to a 2p2 ∈ [1, 4], Z −1 T (5HYP) will jump to 0.98(10) with χ 2 = 0.4 (the tuned fit). The data of Z −1 T (5HYP) (the green triangles) with the band from the default fit (the green band) and tuned fit (the blue band) are plotted in the left panel of Fig. 10, and the renormalizedR(t f ) with both fits of Z −1 T are shown in the right panel. The errors from Z T and the bareR(t f ) are combined in quadrature. It is obvious that the renormalizedR(t f ) with 5-step HYP smearing (green triangles) based on the default fit of Z −1 T is much higher than those with 1,2 steps of HYP smearing. Even though the consistency can be improved if the tuned fit of Z −1 T is applied (the blue boxes), the systematic uncertainties from the fit of Z −1 T will make the final uncertainties in the 5-step HYP smearing case larger than the cases with fewer steps of HYP smearing. and renormalizedR(t f ) with 1,2,5 steps of HYP smearing are shown as the red dots, purple reversed triangles and green triangles, respectively. Both the blue and green bands are the fit of the 5-HYP data, with the regions a 2p2 ∈ [1, 4] and [1.5,5] respectively. The renormalizedR(t f ) in the 2HYP case still consistent with the 1HYP case even though the a 2p2 dependence of Z −1 T are quite different for t f ≥3, but the 5HYP case will be very sensitive to the fit of Z −1