Proton Mass Decomposition from the QCD Energy Momentum Tensor

We report results on the proton mass decomposition and also on related quark and glue momentum fractions. The results are based on overlap valence fermions on four ensembles of $N_f = 2+1$ DWF configurations with three lattice spacings and three volumes, and several pion masses including the physical pion mass. With fully non-perturbative renormalization (and universal normalization on both quark and gluon), we find that the quark energy and glue field energy contribute 33(4)(4)\% and 37(5)(4)\% respectively in the $\overline{MS}$ scheme at $\mu = 2$ GeV. A quarter of the trace anomaly gives a 23(1)(1)\% contribution to the proton mass based on the sum rule, given 9(2)(1)\% contribution from the $u, d,$ and $s$ quark scalar condensates. The $u,d,s$ and glue momentum fractions in the $\overline{MS}$ scheme are in good agreement with global analyses at $\mu = 2$ GeV.

The H m , H E , and H g in the above equations denote the contributions from the quark condensate, the quark energy, and the glue field energy, respectively: The QCD anomaly term H a is the joint contribution from the quantum anomaly of both glue and quark, All the H are defined by N |H|N / N |N where |N is the nucleon state in the rest frame. Note that H E +H g , H m and H a are scale and renormalization scheme independent, but H E (µ) and H g (µ) separately have scale and scheme dependence. The nucleon mass M can be calculated from the nucleon two-point function. If one calculates further H m and H E (µ), then H g (µ) and H a can be obtained through Eqs. (1) and (2). The approach has been adopted to decompose the S-wave meson masses to gain insight about contributions of each term from light mesons to charmoninums [2]. But the mixing between H E (µ) and H m will be non-trivial under the lattice regularization, when there is any breaking of the quark equation of motion at finite spacing. On the other hand, if we obtain the renormalized quark momentum fraction x R q in the continuum limit, and define the renormalized quark energy H R E in term of x R q and H m with the help of the equation of motion, i.e., arXiv:1808.08677v2 [hep-lat] 6 Sep 2018 then the additional mixing can be avoided. Similarly, the renormalized glue field energy can be accessed from the glue momentum fraction x R g by In the present work, we use the lattice derivative operator for the quark EMT and combination of plaquettes for the gauge EMT and address their normalization in addition to renormalization and mixing. We calculate the proton mass and the renormalized x q,g on four lattice ensembles, and extrapolate the results to the physical pion mass with a global fit including finite lattice spacing and volume corrections. Then we combine previously calculated H m [3] to obtain H a from Eq. (2), and the full decomposition of the proton energy in the rest frame as shown in Eq. (1).
Numerical setup: We use overlap valence fermions on (2 + 1) flavor RBC/UKQCD DWF gauge configurations from four ensembles on 24 3 × 64 (24I), 32 3 × 64 (32I) [4], 32 3 × 64 (32ID) and 48 3 × 96 (48I) [5] lattices. These ensembles cover three values of the lattice spacing and volume respectively, and four values of the quark mass in the sea, which allows us to implement a global fit on our results to control the systematic uncertainties as in Ref. [3,6]. Other parameters of the ensembles used are listed in Table I. The effective quark propagator of the massive overlap fermion is the inverse of the operator (D c + m) [7,8], where D c is chiral, i.e. {D c , γ 5 } = 0 [9] and its detailed definition can be found in our previous work [10][11][12]. We used 4 quark masses from the range m π ∈(250,400) MeV on the 24I and 32I ensembles, and 6/5 quark masses from m π ∈(140,400) MeV on the 48I/32ID ensemble respectively which have larger volumes and thus allow a lighter pion mass with the constraint m π L > 3.8. One step of HYP smearing is applied on all the configurations to improve the signal. Numerical details regarding the calculation of the overlap operator, eigenmode deflation in inversion of the quark matrix, and the Z 3 grid smeared source with low-mode substitution (LMS) to increase statistics are given in [10][11][12][13].
Proton mass: We first calculate the proton mass on these four ensembles and apply the SU(4|2) mixed action 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0. 16  HBχPT functional form [14] to fit the results, where M 0 , C 1,2,3,4 , the axial vector coupling g A and an additional partially quenched one g 1 are free parameters, f π = 0.122(9) GeV is the pion decay constant, m v,sea π is the valence/sea pion mass respectively, m pq π = (m v π ) 2 + (m sea π ) 2 + ∆ mix a 2 is the partially quenched mass with the mixed action term ∆ mix a 2 , and a is the lattice spacing. The O(m 3 π ) logarithm function F in the original functional form is dropped since it turns out to be not useful to constrain the fit. Note that we used C I 3 for the 24I/48I/32I ensembles and C ID 3 for 32ID ensemble as they used different gauge actions. We get the prediction of the proton mass at the physical point as M (m phys π , m phys π , 0, ∞)=0.960(13) GeV with χ 2 /d.o.f.=0.52. From the fit, we can also get the light quark mass sigma term H m,u+d ∂M ∂mπ m π /2 = 52(8) MeV which is consistent with our previous direct calculation 46 (7)(2) MeV [3]. The g A we get from the fit is 0.9(2) which is consistent with the experimental result 1.2723(23) [15] within 2σ. Alternatively, using the experimental value of g A predicts the proton mass as 0.931 (8) with a χ 2 /d.o.f. of 1.5. The results of the proton mass with the partially quenching effect (m sea π = m v π ) subtracted are plotted in Fig. 1 as a function of the valence pion mass, together with the blue band for our prediction in the continuum limit. The difference between the results with different symbols reflects the discretization er-rors and finite volume effects which are reasonably small, as shown in Fig. 1.
Momentum fraction: The quark and gluon momentum fractions in the nucleon can be defined by the traceless diagonal part of the EMT matrix element in the rest frame [16], In practice, we calculated ratios of the three-point function to the two-point function where χ S is the standard proton interpolation field with Gaussian smearing applied to all three quarks, and Γ e is the unpolarized projection operator of the nucleon. All the correlation functions from the source points x in the grid G are combined to improve the the signal-to-noise ratio (SNR) [12]. When t f is large enough, R q,g (t f , t) approaches the bare nucleon matrix element N |T q,g 44 |N .
For each quark mass on each ensemble, we construct R(t f , t) for several sink-source separations t f from 0.7 fm to 1.5 fm and all the current insertion times t between the source and sink, combine all the data to do the two-state fit, and then obtain the matrix elements we want with the excited-states contamination removed properly. The more detailed discussion of the simulation setup and the two-state fit can be found in our previous work [3,6,17].
To improve the signal in the disconnected insertion part of x q,g , all the time slices are looped over for the proton propagator. For x g , the cluster-decomposition error reduction (CDER) technique is applied, as described in Refs. [18,19].
The renormalized momentum fractions x R in the MS scheme at scale µ are where x u,d,s,g is the bare momentum fraction under the lattice regularization, and the renormalization constants in the MS at scale µ are defined through the RI/MOM scheme and . Note that the iso-vector matching coefficient R QQ ( µ µ R ) has been obtained at the 3-loop level [20] while just the 1-loop level results of the other R's are available [21].
We list the renormalization constants for T q,g 44 at MS 2 GeV in Table. II and the details of the NPR calculation are discussed in the supplementary materials [22].
After the renormalization, the total momentum fraction is generally larger than 1 by 20%-30% on the four ensembles. Thus the normalization effect here due to the discretization errors is large and we apply a uniform normalization on both the quark and glue momentum fractions at each quark mass of each ensemble, and plot these normalization factorsZ = x −1 u+d+s+g in the lower-right panel of Fig. 2.
Then the pion mass dependence of the renormalized and normalized x R u,d,s,g are fitted with the following em-pirical form simultaneously, , and 0.482(69)(48) respectively, where the first error is the statistical one and the second error includes the systematic uncertainties from the chiral, continuum, and infinite volume interpolation/extrapolations. The systematic uncertainties from the two-state fit and CDER for x g haven't been taken into account yet and will be investigated in the future. With the normalization factor shown in lower-right panel of Fig. 2. All the predictions of the momentum fractions are consistent with the phenomenological global fit at MS 2 GeV, e.g., CT14 [23] values x R u = 0.348(5), x R d = 0.190(5), x R s = 0.035(9) and x R g = 0.416 (9). The other global fits results [24][25][26][27][28] summarized in Ref. [29] are consistent with CT14. After the partially quenching effect term proportional to D 2 is subtracted, the x R u,d,s,g at different ensembles and valence quark masses are illustrated in Fig. 2 as a function of m 2 π , in the upper-left panel for the u and d cases and the upper-right panel for the s and g cases. The bands on the figures show our predictions in the continuum limit with their uncertainties (blue for the statistics and cyan for the total).
We also predict the iso-vector momentum fraction x R u−d as 0.151 (28)(29), which is consistent with the CT14 result 0.158(6) [23], in the lower-left panel of Fig. 2.
Final proton mass decomposition: With these momentum fractions at MS 2 GeV, we can apply Eqs. (5) and (6) to obtain the quark and glue energy contributions in the proton mass (or more precisely, the proton energy in the rest frame). Combined with the quark scalar condensate and trace anomaly contributions, the entire proton mass decomposition is illustrated in Fig. 3 as a function of the valence pion mass. As shown in the figure, the major quark mass dependence comes from the quark condensate term, and the other components are almost independent of the quark mass. At the physical point, the quark and glue energy contributions are 32(4)(4)% and 36(5)(4)% respectively. With the quark scalar condensate contribution of 9(2)(1)% [3], we can obtain that a quarter of the trace anomaly contributes 23(1)(1)% with N f = 2 + 1.
In summary, we present a simulation strategy to calculate the proton mass decomposition. The renormalization and mixing between the quark and glue energy can be calculated non-perturbatively, and the quark scalar condensate contribution and the trace anomaly are renormalization group invariant. Based on this strategy, the lattice simulation is carried out on four ensembles with three lattice spacings and volumes, and several pion masses including the physical pion mass, to control the respective systematic uncertainties. With nonperturbative renormalization and normalization, the individual u, d, s and glue momentum fractions agree with those from the global fit in the MS scheme at 2 GeV. Quark energy, gluon energy, and quantum anomaly contributions to the proton mass are fairly insensitive to the pion mass up to 400 MeV within our statistical and systematic uncertainties.
where x u,d,s,g is the bare momentum fraction under the lattice regularization, and the renormalization constants at MS scale µ are defined through the RI/MOM scheme at scale µ R , and In the above equations, the RI/MOM renormalization constants Z of T q,g 44 are defined with the following conditions suggested by Ref. [20] for cases with quark external legs: where V is the lattice volume, p is the momentum of the external quark/gluon state,p µ = sinp µ , Γ q µµ = iγ µpµ − ip 2 μ p 2p / as suggested by Ref. [30], the quark propagators areS q (p) ≡ S q (p) ≡ x e ipx S q (p, x) with S q (p, x) = y e −ipy ψ(x)ψ(y), and Z q is defined through the axial-vector vertex correction and Ward identity [31]. Note that the index µ in Eqs. (16)- (18) is not summed while the results with different values of µ can be averaged.
For the case of gluon external legs, the definitions are the following, as inspired by Refs. [19,21], where ξ ≡ µ p 4 µ ( µ p 2 µ ) 2 , and Z a/b are defined by with all the repeated indices summed. The 3-loop level result of the matching coefficient R QQ ( µ µ R ) = 1 + g 2 16π 2 C F [ 8 3 log(µ 2 /µ 2 R ) + 31 9 ] + O(α 2 s ) has been obtained in Ref. [20], while just the 1-loop level results of the other R's are available [21]: Thus we will use the 3-loop matching for the renormalization of the quark momentum fraction and keep all the matching at 1-loop level in the rest of the renormalization calculation.
In the practical lattice calculation, all the Z's suffer from discretization errors and we need to repeat the calculation of Z at different p 2 , match them to the MS scheme at the µ R scale, evolve them from µ R to a fixed scale such as 2 GeV, and then apply the a 2 µ 2 R = a 2 p 2 extrapolation to get the final result of Z. To apply this extrapolation properly, we should calculate the vertex corrections in the quark and gluon states with exactly the same momenta, and combine them first. But it is not necessary for most of the cases except for the quark to gluon mixing, as we will discuss case by case.   For the quark external state, we choose 12 momenta with µ p 4 µ ( µ p 2 µ ) 2 < 0.27 for 30 configurations and 6 valence quark masses at each lattice spacing, and use Landau gauge fixed momentum volume sources to improve the signal. With the given momenta, we extrapolated the result to the chiral limit and then applied the continuum matching to get the value at MS 2 GeV. The values of Z QQ R QQ for three lattice spacings are plotted in the left panel of Fig. 4, and the values at a 2 p 2 = 0 limit are obtained based on linear extrapolation of the data in the range a 2 p 2 ∈ [4,8]. It turns out that the a 2 p 2 dependence is small by usingp µ in the definition of the tree level operator and projection, when the chosen momenta are along the body-diagonal direction. The calculation on the 48I ensemble is skipped as its lattice spacing is almost the same as that of the 24I ensemble.
We followed the same strategy used in Ref. [31,32] to analyze the systematic uncertainties and the error budgets are collected in Tab. III. The systematic uncertainty of m sea s = 0 are estimated as ∼1% on the 24I and 32I ensembles and ∼0.25% on 32ID ensemble, based on our previous study with multiple sea quark masses [31]. The values for the off-diagonal parts of the quark EMT cases at a=0.143, 0.111 and 0.083 fm (at a 2 p 2 = 0 limit) are slightly different. They are 0.786(2)(9), 0.796(1)(11) and 0.793(1)(11) respectively.
In the right panel of Fig. 4, the other two terms needed by the singlet quark EMT renormalization, δZ QQ R QQ and Z QG R GQ for the disconnected contribution, are plotted. Regardless of the simulation details of the Z QG which will be addressed later, the contribution of second term is much smaller than the first term, and thus the a 2 p 2 extrapolation can be carried out with the δZ QQ R QQ term alone in the a 2 p 2 ∈ [4,8] range, considering the value of the other term as a systematic uncertainty which is ∼0.001. Note that the CDER technique [18] is applied to the correlation function of δZ QQ with the cutoff r 0 ∼ 1 fm, In our previous investigation of the glue EMT renormalization [19], we chose the momenta with two transverse components to calculate Z GG : But the discretization errors at the range a 2 p 2 ∈ [4,8] used in the quark external legs are too large to carry out a reliable fit, when we want to combine it with Z GQ . Thus in this work, we introduce two improvements to suppress the a 2 p 2 corrections: 1. The a 2 p 2 correction from the HYP smearing on the glue operator: We found that most of the a 2 p 2 correction in Z GG with the HYP smeared gauge EMT can be removed by the following ratio f (a 2 p 2 ), where is the HYP-smeared gauge potential defined from the HYP-smeared gauge link U ρ , and f (a 2 p 2 ) is the ratio of two propagators S(p) ≡ Tr[A ρ (p)A ρ (−p)] with and without HYP smearing.
Since the scale dependence of f (a 2 p 2 ) is cancelled, only the a 2 p 2 discretization errors exist up to the normalization f (a 2 p 2 → 0). Note that such a normalization corresponds to the tadpole effect of the HYP smearing which should be included in the renormalization constant Z GG . Since the value of S(p) doesn't exist at p 2 = 0, we will have to fit f (a 2 p 2 ) with a polynomial of a 2 p 2 to get the value of f (a 2 p 2 → 0). 2. We also extend the calculation of the Z GG to the momenta along the body-diagonal direction, by considering the projected renormalization constants defined in Eq. 20. They can be rewritten into the combination of the renormalization constant of the traceless diagonal gauge EMT Z GG , and the off-diagonal one Z of f GG by FIG. 5. The gauge EMT renormalization constant, with (left panel) and without (right panel) the a 2 p 2 improvements described in Eq. (26) and (27). The contribution fromZ RI GG RGG is dominant and that from N f ZGQRQG is negligible. The a 2 p 2 → 0 results with and without improvement can provide consistent predictions.
With 949/356/290 configurations respectively, the values ofZ RI GG R GG on 32ID/48I/64I ensembles with 1 step of HYP smearing are illustrated in the left panel of Fig. 5 . Note that the ensemble 64I, which has almost the same lattice spacing as 32I (0.0837(2) vs. 0.0828(3)) but with a larger volume and physical pion mass, is used for the calculation, as a similar calculation on 32I would require ∼5,000 configurations, which are not available, to reach the similar accuracy. After both improvements above are applied, the a 2 p 2 dependences of theZ RI GG R GG are mild. The contributions from the other term N f Z GQ R QG are also illustrated on the same figures. Both the values and their uncertainties (∼0.005) are much smaller than the statistical uncertainty ofZ RI GG and thus can be dropped safely. With linear extrapolation of the data in the range a 2 p 2 ∈ [2, 7] and the polynomial fit of the f (a 2 p 2 ) in the same range, we obtain the renormalization constants at a=0.143, 0.114 and 0.084 fm to be 0.79(3)(6), 0.94(3)(3), and 0.91(3)(4) respectively, with the second error determined from varying the starting/ending points of the a 2 p 2 range by 1.
For comparison, we also illustrate the results obtained from the previous strategy used in Ref. [19] in the right panel of Fig. 5, with the fit of the data in the range of a 2p2 ≡ (sin pa 2 ) 2 ∈ [1. 5,5]. The a 2 p 2 → 0 results of based on the polynomial fit are consistent with the present linear-fit ones.  For the mixing from glue to quark, both the contributions from Z QQ R QG and Z QG R GG are small, as illustrated in the left panel of Fig. 6. The total contribution can be estimated as −0.010(10), −0.005(10), −0.000(10) per flavor at a=0.143, 0.111 and 0.083 fm respectively, when the linear extrapolation is applied in the range a 2 p 2 ∈ [4,8].
The last piece of the jigsaw is the mixing from quark to glue. The CDER technique used for δZ QQ as in Eq. 24 can be applied for the calculation here with the gauge EMT operator. As in the right panel of Fig. 6, the contributions fromZ RI GG R GQ are just a few percent, while those from Z GQ R QQ are sizable. Since the a 2 p 2 dependence ofZ RI GG is small as in the left panel of Fig. 5, we can approximate theZ RI GG (a 2 p 2 ) by its a 2 p 2 extrapolated value and assign ∼ 3% systematic uncertainty from the difference between the values at a 2 p 2 equal to 4 and 8. After combining the correction and applying the linear extrapolation in the range a 2 p 2 ∈ [4,8], we obtain the total mixing at a=0.143, 0.111 and 0.083 fm as −0.34(2)(4), −0.26(2)(3), and −0.13(2)(1) respectively.