A Demonstration of Hadron Mass Origin from QCD Trace Anomaly

Quantum chromodynamics (QCD) claims that the major source of the nucleon invariant mass is not the Higgs mechanism but the trace anomaly in QCD energy momentum tensor. Although experimental and theoretical results support such conclusion, a direct demonstration is still absent. We present the first Lattice QCD calculation of the quark and gluon trace anomaly contributions to the hadron masses, using the overlap fermion on the 2+1 flavor dynamical Domain wall quark ensemble at $m_{\pi}=340$ MeV and lattice spacing $a=$0.1105 fm. The result shows that the gluon trace anomaly contributes to most of the nucleon mass, and the contribution in the pion state is smaller than that in others nearly by a factor $\sim$10 since the gluon trace anomaly density inside pion is different from the other hadrons and the magnitude is much smaller. The gluon trace anomaly coefficient $\beta/g^3=-0.056(6)$ we obtained is consistent with its regularization independent leading order value $(-11+\frac{2N_f}{3})/(4\pi)^2$ perfectly.

Quantum chromodynamics (QCD) claims that the major source of the nucleon invariant mass is not the Higgs mechanism but the trace anomaly in QCD energy momentum tensor. Although experimental and theoretical results support such conclusion, a direct demonstration is still absent. We present the first Lattice QCD calculation of the quark and gluon trace anomaly contributions to the hadron masses, using the overlap fermion on the 2+1 flavor dynamical Domain wall quark ensemble at mπ = 340 MeV and lattice spacing a =0.1105 fm. The result shows that the gluon trace anomaly contributes to most of the nucleon mass, and the contribution in the pion state is smaller than that in others nearly by a factor ∼10 since the gluon trace anomaly density inside pion is different from the other hadrons and the magnitude is much smaller. The gluon trace anomaly coefficient β/g 3 = −0.056 (6) we obtained is consistent with its regularization independent leading order value (−11 + Introduction: The most essential and non-trivial feature of the quantum field theory, is the regularization. It introduces a tiny and artificial modification on the short distance part of the original theory, to make the contributions from the virtual particle loops to be finite and calculable, while sensitive to the details of the modification. Most of the regularization efforts can be absorbed into the renormalization of fields and the coupling constants, and then irrelevant to the long distance physics. Thus the exceptional cases are named as "quantum anomaly". Taking the trace anomaly of the QCD energy momentum tensor (EMT) with quark field ψ and gluon field strength tensor F µν as an example, the gluon part of the original EMT is traceless, but it is unavoidable to introduce a non-vanishing trace term on the EMT during the regularization, likes the dimensional regularization case discussed in Ref. [1]. The renormalization eliminates most of the regularization dependence in the additional trace terms and convert them into something proportional to the strong coupling constant, and then the renormalized EMT trace can be expressed as where the quark mass term H m = q m qq q is the classical trace of EMT, and both the other terms are the trace anomaly. Note that the F 2 here is defined in the euclidean space and its vacuum expectation value is positive. The anomalous terms are proportional to the anomalous dimension of the quark mass m q , γ m = − µ m ∂m ∂µ = 2 π α s + O(α 2 s ), and also that of the strong cou- The trace anomaly leads to the most non-trivial feature of QCD: the quantum particles likes nucleon can have positive masses, even though the classical waves of quark and gluon travel at the speed of light. The argument is simple: the hadron mass is proportional to its matrix element of the EMT trace [2][3][4], with O H ≡ H|O|H ; thus the gluon trace anomaly β 2g F 2 H can contribute a positive hadron mass even in the chiral limit where m q → 0 and H m H vanish. The only exception is pion: its mass will be zero in the chiral limit and then the percentage of trace anomaly contribution in pion mass should be less than 50% [5]. One can also use this picture to understand the QED trace anomaly effect in the lamb shift, while it turns out to be small [4,6].
The above argument can be verified indirectly through the sum rule in Eq. (2): The nucleon mass M N can be measured at 10 −10 accuracy, the O H from three light flavors can be extracted from the SU(3) flavor breaking of the baryons and/or explicit Lattice QCD calculations. Then their difference (∼ 850 MeV [7,8]) should come majorly from the trace anomaly. The asymptotic freedom and confinement of QCD require the QCD coupling to be stronger at larger distance and then β is negative, thus arXiv:2101.04942v2 [hep-lat] 7 Jun 2021 F 2 H should also be negative to satisfy the sum rule, even though the matrix element F 2 in the vacuum will be positive. Note that the heavy quark contribution is canceled by the flavor dependence of β 2g at leading order based on the heavy quark approximation [2,9], and then is effectively decoupled. The experimental measurement of trace anomaly arouses great interest, and is considered as one of the major scientific goals of the future Electron-Ion Collider (EIC) [10] and EicC [11]. One approach to detect trace anomaly is measuring the cross section of exclusive photo-production of heavy quarkonium which depends on the gluon condensate at non-relativistic limit [12][13][14]. The gluon gravitational form factors which can relate to J/Ψ production amplitude in the large momentum limit, is also a possibility worth further investigation [15]. There are also other theoretical discussion on the related experimental possibility [16,17].
In theoretical side, the central challenge is verifying whether the difference between hadron mass and the its quark mass contribution actually comes from the quantum anomaly effort. Such a calculation is highly nonperturbative and then can only be done with Lattice QCD, but the chiral symmetry breaking in the quark mass term of the Wilson-like action can mixed with the original trace anomaly. Thus the expensive chiral fermion is essential and there is no explicit verification yet.
In this work, we present the first demonstration of the quantum trace anomaly contribution to kinds of the hadron masses, after this mass generation mechanism has been proposed for more than 40 years. We confirm that Eq. (2) is satisfied in all hadrons we calculated in this work with the γ m and β determined non-perturbatively. In addition, we also find that the gluon trace anomaly density in the pion turns out to be much smaller than that in the other hadron likes nucleon and vector meson, due to the significant difference on the gluon trace anomaly distribution inside the hadrons. Numerical setup: We preform the calculation on a 24 3 × 64 2+1 flavor Domain-wall fermion ensemble from the RBC collaboration [18], with the pion mass m π = 340 MeV which is heavier than the physical one but still much smaller than the corresponding nucleon mass. The other information of the ensemble we use summarized in Table I. For the valence quark, we use the chiral fermion through the overlap approach [19] to avoid the additional term in the trace of EMT under the lattice regularization.
In order to determine the γ m and β precisely and make an accurate test of the sum rule, we consider the partially quenched QCD [20] which allows the valence quark mass to be different from that in the gauge ensemble. In such a case, Eq. (2) becomes where H g a = β 2g F 2 , H m v H includes the connected quark diagram with the operator m vqv q v only, and the index i just includes the flavors exist in the gauge configurations: degenerated light up/down and also strange quark.
In this work, we considered the cases with 5 quark masses m v a=0.0160, 0.0576, 0.1, 0.2, and 0.3, and the lightest two quark masses correspond to the light and strange quark masses in the gauge ensemble we used.
First of all, the hadron mass can be extracted from the two point correlation function with wall source and wall sink where H is the hadron interpolation field. As implemented in Ref. [21], when we calculate the contribution of valence quark mass, we use the Feynman-Hellmann theorem inspired method [22] to construct the summed three point for the connected insertion case, where O qv = m v qqv q v is the valence quark operator. For the disconnected insertion case, we modified the expression into where the sum of the current time slides in the light sea quark operator O q = m sqs q s (m s is sea quark mass) and gluon operator O g = F 2 are limited to t ∈ (t f , 0) to remove the statistical uncertainty from the unphysical region t < 0 and t > t f , and the gluon field tensor F µν is defined through the standard clover definition. The detailed expression of C 2 and SC 3 and how to construct them with propagators, are shown in the supplemental materials [23].
For each hadron, we carry out a joint correlated fit of C 2 (t f ; H) and SC q,g 3 (t f ; H) to extract the M H , H m H and F 2 H simultaneously. As mentioned above, SC qv,qs,g 3 (t f ; H) is summed three point correlation function which has summed the total contribution of matrix element of the interpolated quark and gluon operator between 0 and t f . The fit expression can be written as [21]: where O qv,qs , the e −δmt f terms are introduced to account for the contamination from higher states, δm, B 0 , B 1 and B qv,qs,g 2,...,4 are free parameters. The contribution from transition between ground state and higher excited state are included in B qv,qs,g 2,...,4 . The above form is equivalent to extract the desired quantities in the large t f limits [21,22], for the PS meson case, to include the loop around effect and describe the data around t f ∼ T a/2.
Result: In this work, we calculated the quark propagator at five different valence quark masses m v a=0.0160, 0.0576, 0.1, 0.2, and 0.3 respectively, and construct the two and three point correlation functions for the nucleon (N), pseudoscalar (PS), and vector (V) mesons. We also looped over all the T = 64 time slides to suppress the statistical uncertainty, and applied 5 steps of the HYP smearing on the gauge operator.
FIG. 1. The differential ratio ∆R q,g P S,V (t f ) of the pseudoscalar and vector mesons with mva = 0.3, which should approaches to the ground state matrix elements Hm P S,V and − F 2 P S,V respectively at the t f → ∞ limit. The bands shows the joint fit prediction using the form defined in Eq. (8) and they agree with the data well.
In Fig. 1, we plotted the effective ratios ∆R q,g (t f ) defined in Eq. (9) for the PS (blue) and V (red) mesons with m v a=0.3. The disconnected light quark contributions are not shown here as their contributions are small. The M H , H m H and F 2 H are obtained using the joint fit defined in Eq. 8 with χ 2 /d.o.f.∼1, and and the bands of ∆R q,g (t f ) predicted by the joint fit are also shown in Fig. 1 and agree with the data perfectly. We can see that the gluon trace anomaly matrix element in the V meson is more than two times of that in the PS meson, even though their quark mass terms H m H = ∆R q (t f )| t f →∞ just differ from each other by 10%. The values of M H , δ m , excited state mass difference δm, fitted matrix elements and χ 2 /d.o.f. of different hadrons with kinds of quark masses can be founded in the supplemental materials [23].
Since the anomalous dimension γ m and β should be independent to the hadron states, we solve the equations and obtain the bare γ m = 0.38(3) and β 2g = −0.08(1). Both γ m and β 2g are significantly smaller than one due to the suppression of the effective coupling constant α s , while the F 2 P S/V is large enough to reverse the naive power counting. At the same time, since g 2 F 2 H is the function of the gauge link without the bare coupling dependence explicitly, it is more natural to divide g 2 on the coefficient β 2g and consider β g 3 = −0.056(6) which is independent of α s at the leading order. Such a value is perfectly agree with the regularization independent leading order β (0) g 3 = (−11+ in Fig. 2. We can see that the R H in all the cases are consistent with one within the uncertainties. Note that the PS and V meson cases with m v a = 0.3 are exactly one since they are the input cases. Such a result verifies the trace anomaly sum rule in Eq. (2), and is consistent with our expectation that both γ m and β are universal. The resulting gluonic trace anomaly contribution H g a H = β 2g F 2 H in the 2+1 flavor ensemble are plotted in Fig. 3. We can see that H g a H in the pion state is generally smaller than that in the other states, especially at the unitary point with 340 MeV pion mass. In such We can see that it is is always small in the PS meson, while approaches to ∼ 800 MeV for the nucleon and vector meson in the chiral limit mv → 0. a case, the gluon trace anomaly contributes about 100 MeV of the pion mass which is ∼30%, but ∼ 800 MeV in the ρ meson and also nucleon. It is exactly what we expected from the trace anomaly sum rule: The trace anomaly contributes most of the hadron masses, except the pion case.
The difference will be much more significant at the physical quark mass. If we use the GMOR relation m 2 π ∝ m q and the Feynman-Hellman theorem H m H = m q ∂mπ ∂mq = 1 2 m π + O(m 3 π ) [5,24], we can estimate the gluon trace anomaly contribution in the physical pion state to be 1 2 (1 − γ m )m π = 43(2) MeV. On the other hand, that in the nucleon at the physical light quark mass will be 816(10) MeV if we just consider the quark mass contribution from three light quarks [7,8].
In order to uncover the origin of this difference, we also investigate the gluon trace anomaly density inside  (12) as ρ H (|r|) can be related to the gluon trace anomaly matrix element H g a H = r d 3 rρ H (|r|) and its squared charge radius , with proper integration like the pion change radius case studied in Ref. [25]. A naive guess is that ρ H would be insensitive to the hadron at small r, while it should decay much faster in the pion than the other hadrons to make the integration and change radius to be significantly smaller. But our lattice QCD calculation has excluded such a possibility, as shown in Fig. 4. First of all, ρ H is much smaller than the leading order estimate of the vacuum expectation value F 2 = 0.012(4) [26], even at small r. At the same time, one can see that the density of pion is much smaller than that of nucleon and vector meson; even more, the density tends to be negative and then have the same sign as that in the vacuum. This means that H g a receives both suppression from the magnitude of density and also short range interaction, and it can make the pion gluon trace anomaly radius to be larger than the other hadrons. More details of distribution calculation can be found in the supplemental materials [23].
Summary and outlook: In this work, we calculated the quark mass term mψψ H and gluon action terms F 2 H in the hadron. Based on the EMT trace anomaly sum rule in the PS and V meson states with m v a=0.3, we determined the bare anomalous dimensions of the quark mass and gluon coupling constant as γ m = 0.38(3) and β g 3 = −0.056(6) (with 5 steps of HYP smearing on the gluon operator) respectively, and confirmed that they are independent on hadron states and quark mass up to the statistical uncertainties. With such γ m and β g 3 , we find that the gluon trace anomaly contribution in the PS meson mass is always much smaller than that in the other hadrons, especially around the chiral limit.
A more accurate check on the trace anomaly mechanism can be carried out by renormalizing the F 2 H nonperturbatively [27] and convert it to that under the MS scheme, and/or deducing the conserved EMT with O(a 2 ) trace term directly from the discretized QCD action. The calculation with different gauge actions with O(a 2 ) difference on the bare g 2 at the same lattice spacing (likes Symanzik action used by the MILC configurations and Iwasaki action used here), and the lattice spacing dependence is also desirable. Besides, a direct calculation of the heavy quark contribution m QψQ ψ Q H can also provides a non-varification through Eq. (3).
A non-perturbative determination of γ m and β g 3 at given lattice spacing opens kinds of the interesting possibility on further studies. First of all, Currently the QCD equation of state (EOS) at finite temperature uses either the lattice spacing dependence of the coupling constants [28], or the gradient flow [29] to determine γ m and β g 3 , and then Eq. 2 can only be hold after the continuum extrapolation. The scheme we demonstrated here allows an accurate determination of EOS at any given lattice spacing and then suppresses the discretization errors. It also allows us to calculate the proton gravity form factor and radius directly at the physical pion mass but not only the chiral limit [30]. At the same time, it suggests that the coupling between hadrons and the bilinear heavy quark operator will be very case sensitive based on the relation Eq. (3). For example, the total trace anomaly contribution including both the γ m and β 2g terms are generally ∼ 800 MeV for most of hadrons except light pseudoscalar mesons. But with γ m ∼ 0.3, we will have F 2 N ∼ 0.84 GeV and F 2 ηc ∼ 0.12 GeV. And then m bψb ψ b N will be larger than m bψb ψ b ηc by a factor of ∼ 5 and should be verified with direct calculations. Another possibility is related to the quark and gluon O(Λ QCD ) "dynamical mass" at small off-shell re-gion [31,32] under the Landau gauge, whose origin and the relation with the hadron mass are still unclear yet. In term of the EMT matrix elements in the parton state, the gauge independent trace anomaly and/or the gauge dependent ones should be the origin of such a mass. Thus practical calculation of the trace anomaly matrix element using the coefficient β g 3 obtained here can unravel the role of quantum anomaly efforts in the dynamically generated parton masses.
Another interesting thing is the gluon trace anomaly densities in the pion and nucleon are comparable in their extent, while their magnitude are very different. Even more, the density in the center of pion is negative, and then have the same sign as that in the vacuum. It open a new gate to understand the relation between pion and vacuum, and why the other hadrons are different from them. ACKNOWLEDGEMENT We thank Heng-Tong Ding, Xiangdong Ji, Luchang Jin, Keh-Fei Liu, Jianhui Zhang, and Jian Zhou for valuable discussions, and the RBC and UKQCD collaborations for providing us their DWF gauge configurations. The calculations were performed using the GWUcode [33,34] through HIP programming model [35]. The numerical calculation has majorly been done on CAS Xiaodao In this section, we will give the expressions of two point and summed three point correlation functions constructed by propagator. In our calculation, we use the coulomb gauge fixed wall source propagator S w ( y, t 2 ; t 1 ) = x S( y, t 2 ; x, t 1 ) and the Feynman-Hellmann propagator [22]S c ( y, t 2 ; t 1 ) = m q x,t S( y, t 2 ; x, t)S w ( x, t; t 1 ) to construct the two point and connected three point correlation functions, where C M are defined as For the meson M(Γ) with the interpolation fieldψΓψ, where I is the unitary matrix, and S( y, t 2 ; x, t 1 ) is the quark propagator from ( x, t 1 ) to ( y, t 2 ). For quark loop and gluon operator, we need calculate the following disconnected three point correlation function, The definition of O qs and O g is in the following content of Eq. (7). The nucleon case with the interpolation filed (u TC d)u is similar and can be obtained without much modifications, where S is defined by (CSC −1 ) T ,C = γ 2 γ 4 γ 5 , and Γ m = 1 2 (1 + γ 4 ) is the unpolarized projector. The summed three point correlation functions for light quark loop and gluon are defined as Note that during the calculation with the overlap fermion, deflating the long-distance subspace of the Dirac operator using its eigenvectors v(λ) with |λ| < Λ QCD are essential to obtain the light quark propagator efficiently [36]. At the same time, we can build the light quark loop operator O q (t) only via those eigenvectors v(λ) with little cost, and the systematic uncertainty from the rest is negligible [7], In practical, 300 pairs of the low lying eigenpairs of the Dirac operator with λ ≤ 0.246 GeV are solved to speed up the propagator calculation and construct the light quark loops. Besides, the clover definition of F µν used for the gluon operator F 2 is the following, where U −ν (x) = U † ν (x − aν) and P [µ,ν] ≡ P µ,ν − P ν,µ .
B. Details of fits to eliminate the excited state contamination Using the power counting of e −δmt f , the summed three function can be rewritten into the following form, SC qv,qs,g And then one can obtain that O qv,qs,g H = C2(t f −1;H) + O(e −δmt f ) in the large t f limit. Using the joint fit of the summed 3pt and 2pt, we extract the matrix element of the ground state. The fit results of χ 2 /d.o.f., mass of the ground state hadron M H , the expectional vale of quark mass q m qq q H and gluon operator -F 2 H and δ m are listed in Table II for all the hadrons used in this work.

C. The trace anomaly density and radius in the hadrons
We calculate the trace anomaly density using the following summed three point correlation function, and then ρ H can be obtained by In practical, we use the joint fit defined in Eq. (8)(9) to eliminate the excited state contamination with the data at finite t f . Generally, the charge radius can be defined through the integration of the normalized densityρ H (r) = ρ H (r) d 3 r ρ H (r ) , But in a finite volume L 3 × T , such a definition can only be accurate whenρ H (r) is negligible at r ∼ L/2, otherwise it will receive additional contribution from ρ H (L − r). Thus access the charge radius through the form factor F g (Q 2 ) as the function of the momentum transfer Q 2 , is a more widely used choice. As shown in the Fig. 5,ρ H is still non-zero at r = L/2, especially for the pion case. Thus integrate it with r 2 to obtain the charge radius can have huge systematic uncertainty. But as in the figure,ρ π (r) is smaller thanρ N (r) for all the r < 1.1 fm. Thus gerenrally we can rewrite the pion charge radius into r 2 g π = d 3 r r 2ρ π (r ) = r 2 g N + d 3 r (r 2 − r 2 0 )(ρ π (r ) −ρ N (r )), (25) where r 0 is the intersection ofρ π (r 0 ) andρ N (r 0 ) and we haveρ π (r 0 ) =ρ N (r 0 ) at r 0 . Then the integrand of the second term in the right hand side will be always positive, and one can expect r 2 g π > r 2 g N .