Gravitational form factors of the pion from lattice QCD

The two gravitational form factors of the pion, $A^{\pi}(t)$ and $D^{\pi}(t)$, are computed as functions of the momentum transfer squared $t$ in the kinematic region $0\leq -t<2~\text{GeV}^2$ on a lattice QCD ensemble with quark masses corresponding to a close-to-physical pion mass $m_{\pi}\approx 170~\text{MeV}$ and $N_f=2+1$ quark flavors. The flavor decomposition of these form factors into gluon, up/down light-quark, and strange-quark contributions is presented in the $\overline{\text{MS}}$ scheme at energy scale $\mu=2~\text{GeV}$, with renormalization factors computed nonperturbatively via the RI-MOM scheme. Using monopole and $z$-expansion fits to the gravitational form factors, we obtain estimates for the pion momentum fraction and $D$-term that are consistent with the momentum fraction sum rule and the next-to-leading order chiral perturbation theory prediction for $D^{\pi}(0)$.


I. INTRODUCTION
Quantum chromodynamics (QCD) [1][2][3][4] provides a rigorous description of hadrons as composite particles made up of quarks and gluons interacting via the strong force.The complexity of QCD, however, is such that constraints on quark and gluon contributions to many aspects of hadron structure are difficult to obtain, and the landscape of experiment and theory efforts on this front continues to evolve.For example, recent years have seen considerable progress in understanding how the physical properties of hadrons, like their mass and internal forces, are generated by their fundamental constituents [5,6].
These aspects of hadron structure can be addressed through determinations of gravitational form factors (GFFs) [5,7].These quantities contain nonperturbative information about the coupling of a hadron state to the energy-momentum tensor (EMT) of QCD, the symmetric part of which can be decomposed as [8][9][10][11] T µν = T µν g + T µν q , where Here F µν is the gluon field strength tensor, ψ f is a quark field of flavor f , D µ = ∂ µ + igA µ , A µ are the gluon fields, g is the strong coupling, and γ µ are the Dirac matrices.
The repeated indices are contracted with the Minkowski space-time metric g µν , the trace is over color space, and a {µ b ν} = (a µ b ν + a ν b µ )/2.The GFFs of a hadron are defined from the matrix elements of the EMT in the hadron state, and can be decomposed into quark and gluon contributions.They encode the distribution of energy, spin, pressure, and shear forces 1 within hadrons [13].The total EMT of QCD is conserved, i.e. ∂ µ T µν = 0, but the individual quark and gluon terms are not, and therefore 1 The physical significance of this interpretation as mechanical forces is debated [12].
the resulting quark and gluon GFFs mix under renormalization and depend on the renormalization scheme and scale.
For a pion state specifically, the matrix element of the symmetric EMT can be decomposed in terms of two GFFs, A π (t) and D π (t), as where |π(p)⟩ is a pion state carrying four-momentum p, P = (p + p ′ )/2, ∆ = p ′ − p, and t = ∆ 2 .In the forward limit, the contributions to the A π GFF define the momentum fraction carried by the quark and gluon constituents, and therefore A π (0) = A π g (0) + A π q (0) = 1, while D π (0) = D π g (0) + D π q (0) is the so-called D-term, which is related to the internal forces of hadrons [13], and is predicted to be −1 for the pion up to chiral-symmetry breaking effects [14][15][16].
The importance of the GFFs in characterizing hadron structure has driven a targeted experimental program in recent years, with the first extractions of proton quark [17] and gluon [18] GFFs achieved from deeply virtual Compton scattering and J/ψ photoproduction measurements respectively.Progress towards the determination of the pion GFFs has been more limited, with the first phenomenological constraints of the pion quark GFFs attained using data from the Belle experiment at KEKB [19][20][21].Further constraints on various hadron GFFs can be expected from current and future facilities, including the JLab 12 GeV program [22][23][24] and the Electron-Ion Collider (EIC) [25].
The theoretical determination of GFFs from first principles is possible through the computation of EMT matrix elements using lattice QCD [26], a numerical framework that defines the QCD path integral on a discrete Euclidean space-time lattice, allowing the calculation of nonperturbative hadronic properties.The pion quark GFFs have previously been calculated using lattice methods in Refs.[27,28].Lattice calculations have also provided predictions for the gluon GFFs of the pion [29,30], for the pion quark momentum fraction [31,32], and for the complete flavor decomposition of the pion momen-  [42], is used for the calculation of the bare matrix elements presented in Sec.II.The calculation of the renormalization coefficients, presented in Sec.III, is performed on ensemble B.
In this work, we present the first lattice QCD calculation of the full flavor decomposition of the pion GFFs A π (t) and D π (t) in the kinematic region 0 ≤ −t < 2 GeV 2 on a single ensemble with N f = 2 + 1 quark flavors, and quark masses corresponding to a close-tophysical pion mass of m π ≈ 170 MeV.The extraction of the bare matrix elements is presented in Sec.II, while the nonperturbative renormalization is discussed in Sec.III.Our final results for the renormalized GFFs are given in Sec.IV, along with a discussion of the forward limits of the GFFs.

II. EXTRACTION OF BARE MATRIX ELEMENTS
In this section, we discuss the extraction of the bare matrix elements of the EMT, which are combined with the renormalization coefficients presented in Sec.III to produce the renormalized GFFs presented in Sec.IV.The bare matrix elements are calculated on a single (2+1)-flavor lattice QCD ensemble [42] of volume L 3 × T = 48 3 × 96, with light quark mass tuned to produce pion mass m π ≈ 170 MeV, and lattice spacing a = 0.091(1) fm [43,44].The ensemble was generated using the Lüscher-Weisz gauge action [45] and cloverimproved Wilson quarks [46] with the clover coefficient set to the tree-level tadpole-improved value and constructed using stout-smeared links [47].The specifics of this ensemble, referred to as ensemble A, are summarized in Table I.The configurations were taken to be independent, and the number used was different for the calculations of the bare quark and gluon contributions, and will be specified in the corresponding subsections below.

A. Two-point functions
We construct momentum-projected two-point correlation functions as where the pion interpolating operator is chosen to be In this expression, ψ u and ψ d are smeared quark fields obtained by applying gauge-invariant Gaussian smearing to radius 4.5a, constructed using links stout-smeared [47] in the spatial directions only.Correlation functions are computed for 1024 source positions (x 0 , t 0 ) on each of a total of 2511 configurations, separated by at least four trajectories.The 1024 sources are arranged in two 4 3 × 8 grids offset by (6,6,6,6) lattice units, with an overall random offset for each configuration.Correlation functions are constructed for all three-momenta p = 2πn/L with |n| 2 ≤ 10, i.e., 147 distinct momenta.
To extract the pion spectrum, we average the twopoint correlation functions on each configuration over all sources, shifting each such that (x 0 , t 0 ) = (0, 0), and form 500 bootstrap [48,49] ensembles from the 2511 measurements, sampling with replacement and sample size for each bootstrap equal to the total number of configurations.We extract the energies by considering the spectral decomposition where t s denotes the sink time and E n p denotes the energy of the nth state with the same quantum numbers as the interpolating operator, |n(p)⟩.Overlap factors are defined as The lowest state in the spectrum (n = 0) is expected to be the pion |π(p)⟩ with energy E π p .To extract E π p , we perform correlated multiexponential fits incorporating up to three states to the two-point functions averaged over all sink momenta with the same magnitude |p|.The effective values of c obtained by taking m π = E π 0 and inverting the dispersion relation E π p = m 2 π + |cp| 2 deviate by 1% at most from unity, as shown in Fig. 1.Both the uncertainties on the determination of the energies E π p and their deviation from the dispersion relation values are negligible in comparison to the statistical uncertainties of the matrix elements that define the GFFs.Thus, the dispersion relation with the central value of am π = 0.0779 and c = 1 is used to set the energies E π p used in the remainder of the analysis.

B. Operators
In discrete Euclidean space-time, symmetric quark and gluon EMT operators can be expressed as where f denotes the quark flavor, the repeated Greek indices are summed over Euclidean components, and the quantities with an E superscript are the Euclidean versions of those defined in Sec.I.The symmetrized lattice covariant derivative and the gluon field strength tensor can be expressed up to discretization effects using where g 0 is the bare lattice coupling 2 , U µ (x) are the lattice gluon link fields, and is the clover term.The Euclidean EMT components are related to the Minkowski ones by for j, k ∈ {1, 2, 3}.To determine the flavor decomposition of the pion GFFs into light and strange quark components, we consider the isosinglet Tq,µν and nonsinglet Tv,µν EMT operators, defined as Tq,µν = Tu,µν + Td,µν + Ts,µν , Tv,µν = Tu,µν + Td,µν − 2 Ts,µν .
The quark isosinglet EMT mixes with the gluon one under renormalization, while the quark nonsinglet EMT does not.Lorentz symmetry is broken on a discrete hypercubic lattice, and traceless linear combinations of diagonal and off-diagonal components of the EMT transform under different irreducible representations (irreps) of the hypercubic group.Two different irreps protected from mixing with lower-dimensional operators are available, τ and τ (6) 3 [50,51].A basis of operators for τ and a basis for τ is Tτ (6)   3,1 = −i √ 2 ( T10 + T01 ), Tτ (6)   3,4 = −i √ 2 ( T20 + T02 ), Tτ (6)   3,6 both in Minkowski space.These bare operators must be renormalized, accounting for mixing between the quark and gluon operators within the same irrep.In the rest of this section, we discuss the computation of matrix elements of the bare lattice operators, while the renormalization and mixing is discussed in Sec.III.

C. Three-point functions
The three-point correlation functions needed in order to isolate the bare matrix elements of the operators of Eqs. ( 12) and ( 13) are where R ∈ {τ 3 } and ℓ runs over the corresponding irrep operator basis.Using translational invariance to set (x 0 , t 0 ) = (0, 0), the spectral representation of the three-point function in the limit where (t s − τ ) ≪ T and t s ≪ T can be expanded as where p = p ′ − ∆.The three-point function contains light-quark connected, light-quark disconnected, strange (disconnected), and gluon terms, The spectral decomposition, Eq. ( 15), holds for each individual piece, 3 allowing the corresponding matrix elements ⟨π(p ′ )| TiRℓ (∆) |π(p)⟩ to be considered separately, where i ∈ {l conn , l disco , s, g}.The factor of 2 multiplying the l conn and l disco terms is due to the identical contributions of the up and down quarks in the case of the pion.
We apply the summation method [52][53][54] to extract the bare matrix elements from the three-point functions.We first form the ratios in which the overlap factors and the time dependence of the ground state terms in Eqs. ( 15) and (5) cancel, up to finite-T effects.When the source, operator, and sink lRℓ always admits a spectral decomposition, and C 3pt,disco lRℓ does because it is identical to a fully disconnected three-point function in a partially quenched theory with an additional light valuence quark, so their difference C 3pt,conn lRℓ does as well.
are well-separated in Euclidean time, the ratio can be expanded as lim ts,(ts−τ )→∞ where ∆E and ∆E ′ are energy differences which depend on ∆ and p ′ .We then sum the averaged ratios over operator insertion time for where Λ is a t s -independent function of momenta and δ is an energy difference which depends on ∆ and p ′ .To improve the signal, we average within each flavor and irrep all ratios for choices (ℓ, p ′ , ∆) that correspond to identical linear combinations of the GFFs up to an overall sign, which we call "c-bins".The coefficients of the GFFs in these linear combinations are defined in Eq. ( 2), rescaled to account for the factor of 2 » E π p E π p ′ in the denominator of Eq. (18).The resulting summed averaged ratios are denoted as ΣiRc .We further partition these into 25 discrete groups, denoted as "t-bins", based on the proximity of their t values, using k-means clustering [55].We form ΣiRct for all contributions to the bare EMT three-point function of Eq. ( 16), and fit them to obtain the corresponding bare matrix elements, ME iRct .The bare GFFs for each irrep R and t-bin are constrained by the system of linear equations where (K A Rct , K D Rct ) are the (unique) coefficients of the c-bin corresponding to matrix element ME iRct , and bold symbols are vectors in the space of c-bins.We discuss the extraction for each contribution i ∈ {l conn , l disco , s, g} individually in the rest of this section.

D. Connected quark contribution
The connected light-quark contribution (i = l conn ) to the three-point function of Eq. ( 16) can be constructed as where Rℓ represents the linear combination of γ µ ← → D ν components corresponding to irrep R and basis element ℓ, defined in Eqs. ( 12) and ( 13), and S l (a; b) are light-quark propagators from lattice coordinate a to b.The symmetric discretized covariant derivative in Eq. ( 21) acts on the (y, τ + t 0 ) argument of the quark propagators directly on its left and right, and shifts them as defined in Eq. ( 8), which is left implicit in the equation for simplicity.Quark smearing at the source and sink is also left implicit.We compute these on 1381 configurations, separated by at least ten trajectories, using the sequential source method, inverting through the sink.On each configuration, we compute correlation functions for 7 distinct temporal source-sink separations t s .The number of source positions computed varies with t s as tabulated in Table II.For every source position and sink time, we compute correlation functions with three different sink momenta, (L/2π)p ′ ∈ {(1, 0, −1), (−2, −1, 0), (−1, −1, −1)}.We use the full set of 9 operators of Eqs. ( 12) and ( 13) and all operator insertion momenta with |∆| 2 < 25(2π/L) 2 , and form 500 bootstrap ensembles from the 1381 source-averaged correlation functions.
In order to extract the bare ground-state matrix elements, we first form the summed ratios of Eqs. ( 17) and ( 19), using source-averaged three-point correlators and two-point correlators computed on the same set of 1381 configurations.Ratios are binned as described in Sec.II C, yielding Σconn lRct with 708 c-bins for τ (3) 1 and 671 for τ (6) 3 .These are fit to the functional form of Eq. ( 19), including the exponential term (cf.Ref. [56]).To enforce a gap over the ground state, the log of ∆E is fit with a prior log ∆E = log[0.2± 0.5], where 0.2 ≈ 3E π − E π is chosen because the first excited state in the spectrum is expected to be a three-pion state.The other parameters are fit with wide, noninformative priors.Model averaging [57][58][59] with Akaike information criterion [60] (AIC) weights over different data cuts is employed to obtain the final estimated matrix element for each cbin.The fit ranges averaged over are τ cut ∈ {2, 3, 4} and t s,min ≤ t s ≤ t s,max , where t s,max ∈ {16, 18}, and t s,min ∈ {6, 8, 10}.Redundant fits where t s,min < 2τ cut are not double-counted.Fits with less than five distinct t s values are not included in the pool of values for averaging.More details are provided in Appendix A 1.

E. Disconnected quark contribution
The disconnected quark three-point function receives a contribution from the light quark and strange quark terms of the EMT.These can be written as: where S s (a; b) is the strange quark propagator.Like in Eq. ( 21), the shifting of the quark propagator on the third line due to the action of the covariant derivative and the quark smearing at the source and sink are left implicit.We measure both light and strange contributions on the same 1381 configurations used for the measurement of the connected contribution discussed in Sec.II D.
We stochastically estimate the trace in the third line of the expression above using one independent sample of Z 4 noise [61] per configuration.We dilute in spacetime using hierarchical probing [62,63] with a basis of 512 Hadamard vectors, and compute the spin-color trace exactly.The three-point function for each configuration is then computed by convolving the trace estimate with the grid of 1024 two-point correlation functions described in Sec.II A to obtain a source-averaged estimate.
To extract the matrix elements, we fit summed averaged ratios Σdisco s/lRct using the functional form of Eq. ( 19) with linear terms only, as the data has insufficient precision to model excited states.We perform fits to summed ratios constructed with all τ cut ≥ 2, for all ranges of 4 or more consecutive sink times in the range [t s,min , t s,max ].The minimum sink time is set to t s,min = 11 for τ 1 , and the maximum sink time to t s,max = 24.These choices are motivated in Appendix A 2. Model averaging with AIC weights [57] then yields the bare matrix elements.
Quark matrix elements that are well-defined under renormalization can be constructed by forming the nonsinglet and isosinglet combinations defined in Eq. ( 11), combining the appropriate connected and disconnected  components.Unlike the connected three-point function dataset, the disconnected data includes all values of t s , which, with the momenta ∆ and p ′ defined above, yield a total of 2179 c-bins for irrep τ 3 .Combining the connected and disconnected data at the level of c-bins would thus involve discarding a considerable fraction of the disconnected data for which no corresponding connected results were computed.Instead, the bare disconnected matrix elements are first fit using Eq. ( 20) to yield the bare GFFs 4 A π,disco,B u+d,Rt , D π,disco,B u+d,Rt , A π,B sRt , and D π,B sRt , with results shown in Fig. 2.These are then used to reconstruct improved estimates of the disconnected matrix elements ME disco u+d,Rct and ME sRct matching the subset of c-bins that are available for the connected quarks.The resulting matrix elements are combined at the bootstrap level with the connected ones, to yield the nonsinglet and isosinglet operator matrix elements, which are renormalized as discussed in Sec.IV.

F. Gluon contribution
The gluon contribution to the three-point function of Eq. ( 16) can be computed as , and similarly for all quantities with a u + d subscript.
where TgRℓ is the gluon EMT projected to the basis defined in Eqs. ( 12) and (13).Quark smearing at the source and sink is left implicit.Gluon operators are computed on 2511 configurations, with the gauge fields first subjected to 200 steps of Wilson flow [64][65][66] up to flow time t flow /a 2 = 2 in order to reduce noise due UV fluctuations.The three-point functions are formed as with the disconnected quark contributions, averaging over 1024 sources on each configuration, vacuum subtracting, and using all momenta with |∆| 2 ≤ 25(2π/L) 2 and |p ′ | 2 ≤ 10(2π/L) 2 , and all sink times t s and operator insertion times τ .
We fit the resulting summed averaged ratios ΣgRct , fitting the linear term only as in the case of the disconnected bare matrix elements.We perform fits to summed ratios constructed with τ cut ≥ 4, avoiding contact terms due to the extension of the flowed operator definition, and to all ranges of t s of length 4 or more in the range [t s,min , t s,max ].As motivated in Appendix A 3, we take t s,max = 25 and choose different values of t s,min for different t-bins due to differing extents of excited state contamination.Model averaging with AIC weights [57] yields the bare matrix elements used to produce the renormalized GFFs of Sec.IV.

III. RENORMALIZATION
In order to obtain renormalized GFFs A π i (t), D π i (t) for i ∈ {g, u, d, s}, we need to compute the renormalization coefficients of Tg,µν , Tq,µν , and Tv,µν .We consider a renormalization scheme defined in the chiral limit where SU(3) isospin symmetry is exact.In this case, the nonsinglet operator Tv,µν is protected from mixing with the gluon operator, and therefore can be multiplicatively renormalized as where T R v,µν and T B v,µν denote the renormalized and bare operators.Z R v (µ 2 , µ 2 R ) denotes a factor that renormalizes the operator in a renormalization scheme R, and runs the result from scale µ R to scale µ.
In contrast, the isosinglet quark and gluon EMTs mix under renormalization, and renormalized operators in scheme R can be defined as (25) To determine the renormalization factor in Eq. ( 24) and matrix in Eq. ( 25), we first compute the renormalization coefficients of each irrep in the RI-MOM scheme [67,68], match them to MS, and run to the scale µ 2 = (2 GeV) 2 .This is the scheme and scale in which we present the renormalized quark and gluon GFFs in Sec.IV.The matching and running is equivalent to multiplication by the conversion factors [69].The inverse nonsinglet renormalization factor can be expressed in terms of the RI-MOM nonsinglet renormalization coefficient R RI vR with R ∈ {τ while the inverse of the isosinglet matrix can be written as 5 In the remainder of this section, we discuss the details of this calculation, which proceeds through the following steps: • For each irrep R, we compute the five different RI-MOM renormalization coefficients that appear in Eqs. ( 26) and ( 27), R RI v and R RI ij with i, j ∈ {q, g}, using renormalization conditions defined in Sec.III A. 5 In our convention, the matrix C RI/MS ij of Eq. ( 27) is the inverse of C RI ′ ij defined in Ref. [69], setting to zero mixing with operators that vanish in matrix elements of physical states [70].In Eq. ( 26),  27) and the nonsinglet renormalization coefficent of Eq. ( 26) for each irrep R.
The RI-MOM renormalization coefficients are computed on ensemble B with parameters as in Table I.While it would be preferable to compute these on ensemble A to match the bare matrix elements presented in Sec.II, far greater statistics than can be practically obtained would be required to make their extraction possible on that ensemble.However, comparison of other renormalization factors on these ensembles [71] suggests that the resulting systematic uncertainty is no more significant than the other unquantified uncertainties in this calculation using a single ensemble of gauge fields.

A. RI-MOM
The RI-MOM scheme [67,68] defines renormalization coefficients by imposing, in a fixed gauge, that amputated forward-limit three-point functions with incoming fields of four-momentum p are equal to their tree-level values at some energy scale µ 2 R = p 2 .For each irrep, four such renormalization conditions are needed to solve for the four coefficients of the mixing matrix in Eq. (25), and an additional one is needed to fix the renormalization of the nonsinglet operator.
Bare three-point functions for the flavor-singlet quark and gluon EMT operators can be constructed with either light quarks or gluons as the external states as where the subscript i ∈ {g, q} denotes the type of operator and the superscript denotes the type of external fields.The variables in the equations above, and in the rest of this section, are written in Euclidean space.The gluon fields are computed from link fields U µ (x) fixed to Landau gauge and can be expressed up to discretization effects as where N c = 3, and μ is a unit vector in the µ-direction.
The bare light quark and gluon propagators needed to amputate the three-point functions are expressed as These are renormalized by The gluon field renormalization is defined as with lattice momenta while the quark field renormalization is defined as We consider the amputated three-point functions, where summation over repeated indices is implied.The tree-level value of the amputated three-point functions with a gluon operator and quark external states or a quark operator with gluon external states are equal to zero at p2 = µ 2 R , while C q,amp q,µν is set equal to at p2 = µ 2 R , and C g,amp g,µναβ is set equal to at p2 = µ 2 R .To isolate the gauge-invariant part of the gluon three-point function, which is proportional to the first term of Eq. ( 38) [72], the momenta and indices must be chosen such that pα = 0, α = β, α ̸ = µ, and α ̸ = ν.The Landau-gauge gluon propagator is not invertible, but for this choice of four-momenta p it takes a block diagonal form of an invertible and a noninvertible piece; for this choice of indices, only the invertible part is needed to amputate the three-point function.
The RI-MOM renormalization coefficients can be obtained by solving the following system of renormalization conditions: The leading-order O(a) correction to Λ q µν (p) calculated in lattice perturbation theory [73,74] is It can be removed by projecting the three-point functions C q,amp i for i ∈ {q, g} to the Λ q µν (p) space.To accomplish this, an inner product ⟨•, •⟩ is introduced on the space of tensors with two Lorentz and two Dirac indices [75], and the renormalization conditions for R RI qq and R RI gq are recast as the matrix equations which are solved for R RI iq with i ∈ {q, g}.To reduce statistical fluctuation in the computation of the renormalization factors involving gluon external states, one can substitute [29,76] ) for one of the gluon propagators in the denominator of the amputated three-point functions.This cancels the dependence of Eq. ( 39) on Z g .
The renormalization condition for the nonsinglet contribution is where C q,amp v,µν is the amputated three-point function of the nonsinglet operator defined in Eq. (11).We project out the leading order lattice artifact of Eq. ( 40) for R RI v using the same approach as in the case of R RI qq and R RI gq described above.

B. Fitting of renormalization coefficients
The For the three-point functions with external quarks, the cut dem ≤ 0.3 both results in acceptable fit quality and retains sufficient data to allow several thousand fits to be performed to different subsets to assess the systematic uncertainties, as discussed in the following subsections.For the three-point functions with external gluons, there is the additional constraint that pα = 0, where α is the Euclidean vector index of the external gluon fields.This reduces the number of available momenta with low democracy, so the cut for these varies between 0.4 and 0.5 instead, as is discussed case-by-case in the following subsections.
The remaining (ap) 2 dependence of the renormalization factors (i.e., due to lattice artifacts and nonperturbative effects) can be modeled as with P n being the polynomial function P log n the polynomial of a logarithm and P inv n the inverse polynomial The renormalization component R RI ij C RI/MS jk can be extracted from Eq. ( 45) as Different R parameters for each (i, j, k) in Eq. ( 45) are implied, and an identical equation can be used to model . In practice, we find that simultaneously including all P n1 , P log n2 , and P inv n3 terms leads to overfitting.The different renormalization terms are instead modeled by different combinations of P n1 , P log n2 , and P inv n3 , according to which momentum dependence(s) is (are) dominant.Specifically: • Logarithmic terms (P log n2 ) are included only for the renormalization terms that are multiplied by the flavor off-diagonal C RI/MS ij for i ̸ = j.In all such cases, we find that a polynomial with n 2 = 2 is sufficient to describe the data, and that including higher-order terms leads to overfitting.Logarithmic terms are not included for the renormalization coefficients multiplied by the flavor-diagonal C RI/MS jj , because the logarithmic dependence for those terms is a subleading correction to the matching coefficient.
• Polynomial terms (P n1 ) are included for renormalization coefficients multiplied by the flavor-diagonal C

RI/MS jj
(including the single isovector contribution), excluding those with unflowed external gluon fields, as discussed below.Including polynomial terms for coefficients multiplied by the flavor offdiagonal C RI/MS ij does not alter the final results, since the logarithmic dependence dominates for these terms.We find that setting n 1 = 1 in the polynomial is sufficient, and including higherorder terms yields orders-of-magnitude smaller AIC weights.
• The inverse polynomial term (P inv n3 ) is included in the fit ansatz for all renormalization coefficients with unflowed external gluon fields.This is necessary because of the strong discretization artifacts of the gluon propagator in the denominator of these coefficients, when it is constructed from unflowed link fields.The parameters of the inverse polynomial term are always chosen such that the term is monotonic in the (ap) 2 region in which the data is fit.The degree of the polynomial n 3 used is different between the different coefficients, and is set to be the integer that yields the highest p-value for each fit.Consistent results are obtained when including several different choices of n 3 in the model averaging.
Further details are provided in the following discussions for each renormalization coefficient.

C. Rv and Rqq
The renormalization condition for R RI vR , Eq. ( 43), depends on the three-point function of the nonsinglet quark EMT, projected to the Euclidean version of the corresponding irrep R basis of Eqs. ( 12) and ( 13), with external u quark fields.Before amputation, this can be expressed for operator ℓ as where The notation is the same as in Fig. 3.We note that R RI gq C RI/MS qg (b) does not correspond to the fit band continued to (ap) 2 = 0, due to the presence of logarithmic terms in the model ansatz, as described in Secs.III B and III D. and x, y, z, p are four-vectors.The u three-point function results in a connected and a disconnected contribution where l = u/d, while the d and s three-point functions are disconnected The connected contribution is computed using the sequential source method, inverting through the operator, on 240 configurations of ensemble B. The light and strange disconnected contributions are computed on 20000 configurations, using hierarchical probing [62]  , as shown in Fig. 3a.The results for each operator ℓ within each irrep are averaged to improve statistics.Fits are performed to all possible sets of 4 or more points with 3 ≤ a 2 p2 ≤ 9, and all fits with p-value > 0.01 are modelaveraged based on their AIC weights.
R RI Rqq is constrained by solving the corresponding renormalization condition of Eq. ( 39) using the threepoint function C q qRℓ (p 2 ) = C q uRℓ (p 2 ) + C q dRℓ (p 2 ) + C q sRℓ (p 2 ) , (54) which also receives a connected and a disconnected contribution.Its computation proceeds as that of R RI v discussed above.It contributes to the MS renormalization matrix of Eq. (27)  are used, with the data and fits presented in Fig. 3c.

D. Rgq
We calculate R RI gqR on 20, 000 configurations of ensemble B by forming vacuum-subtracted three-point functions of the gluon EMT projected to irrep R with external quark propagators.The operator is constructed using gluon link fields flowed to t flow /a 2 = 1.2 to match the physical scale of the flow radius used for the bare gluon operator on the finer ensemble A of Table I, as detailed in Sec.II F. The quark propagator calculations are

R RI
qg is calculated on 20, 000 configurations, using threepoint functions formed by the disconnected singlet quark operator computed as described in Sec.III C. The gluon propagators are formed as defined in Eq. ( 31) for the diagonal indices α = β using gluon fields computed from Eq. ( 29), and are projected to all four-momenta p = 2πn/L that are subject to the constraints µ n 2 µ ≤ 36, n τ = 0, and dem ≤ 0.5.We form R RI qg in two ways, using gluon propagators computed on gluon link fields that are either flowed (to t flow /a 2 = 1.2) or unflowed before gauge fixing.The two approaches are expected to have different discretization artifacts, but to give consistent results once the (ap) 2 dependence has been fit [29,78].In both cases, within each irrep we average all operators and gluon field polarizations α that yield the same R RI qg (a 2 p2 ) as defined in the renormalization condition of Eq. (39).RRI qg C RI/MS gg formed with unflowed external gluon fields is modeled by P inv 2 of Eq. ( 45), as discussed in Sec.III B.

For RRI qg C
RI/MS gg formed with flowed external gluon fields, we find that the discretization effects become milder with increasing flow time, and are well-described by a linear fit [P 1 in Eq. ( 45)] when t flow /a 2 = 1.2.The fitting of both the flowed and unflowed versions yield consistent results for R RI qg C

RI/MS gg
, and the final result is obtained by combined fits to both, sharing only the constant parameter.Fits are performed to all p with pmin ≤ p ≤ pmax , varying the boundaries to all points between 1.5 ≤ a 2 p2 min ≤ 3.5 and a 2 p2 max ≤ 7. RRI qg C RI/MS gq with unflowed external gluon fields similarly demonstrates beyond-linear discretization effects along with logarithmic dependence, and is found to admit high-quality fits when the ansatz P log 2 + P inv 1 of Eq. ( 45) is used, as discussed in Sec.III B. The equivalent quantity formed with external gluon fields flowed to t flow /a 2 = 1.2 demonstrates primarily logarithmic dependence, and only P log

F. Rgg
We calculate R RI gg on 20, 000 configurations by forming vacuum subtracted three-point functions of the flowed gluon operator with external gluon propagators.The computation of the gluon operator is as described in  Sec.III D and of the gluon propagator as in Sec.III F. As in Sec.III E, we form R RI gg using external gluon fields that are either unflowed or flowed to t flow /a 2 = 1.2.In both cases, we average within each irrep all operators and gluon field polarizations α that yield the same R RI gg (a 2 p2 ) as defined in the renormalization condition of Eq. ( 39). with unflowed gluon fields as the external states is fit using P inv 3 , using the same democracy cut as the results above, but with 1 < a 2 p2 min < 2. Restricting the fit to start at small momenta is found to be necessary due to the strongly decaying behavior.The flowed version of RRI gg C RI/MS gg is well-described by a linear model [P 1 of Eq. ( 45)] when restricted to the dataset with momenta of dem ≤ 0.5 and

Similarly to the case of RRI
RI/MS gg is extracted from a combined fit to both versions, with the results shown in Fig. 6b.

G. Renormalization factors
The final results for the isosinglet and gluon inverse MS renormalization matrices of the two irreps are τ R MS qq = 1.056 (27), R MS qg = 0.067(71) , R MS gq = −0.169(22), R MS gg = 1.68 (18) , τ We note that the large values for R MS gg , which contribute to small values around 0.6 for the corresponding component of the inverse matrix Z MS gg , are due to the flowing of the gluon operator; the bare gluon GFFs also increase with increasing flow time, which is compensated by the decreasing renormalization coefficient.

IV. RENORMALIZED RESULTS
The procedure described in Sec.II yields a set of measurements which constrain the bare 2-dimensional vector of GFFs G π,B iRt for each irrep R ∈ {τ 3 } and flavor i ∈ {q, g, v} separately as where ME iRt is the c-dimensional vector of matrix element fits extracted from the lattice data using the summation method, and K t is the kinematic coefficient matrix written as a (c × 2)-dimensional matrix in the space of c-bins and pion GFFs.For i ∈ {q, g}, the bare GFFs are defined in terms of the renormalized GFFs, G π it , as where R MS ijR are the components of the renormalization matrix with values listed in Eq. (55).Since both irreps are expected to yield the same G π it after renormalization, they can be fit simultaneously, with the individual irrep renormalization coefficients used as inputs in the fit.We therefore recast Eq. ( 57) and ( 58) into a single linear system for both irreps where We fit Eq. ( 59) for G π t by minimizing where Cov(ME t ) is the covariance matrix of the matrix element fits for momentum bin t, computed from the bootstrap samples, and RK MS t is assumed to be Gaussian with the mean and covariance determined by the RI-MOM fitting procedure of Sec.III.The solution takes the analytic form and the errors are determined using Gaussian error propagation.This circumvents the d'Agostini bias [79] by neglecting additional correlations between renormalized constraints induced by common factors of R MS ijR .
For the renormalized nonsinglet GFF for which mixing is not present, we follow a combined irrep fitting procedure that is a Bayesian version of the "penalty trick" [79], and was introduced in Ref. [30] for the renormalization of gluon GFFs when mixing with quarks is assumed to be negligible.For a detailed discussion of this procedure, see Appendix 6 of Ref. [30].

A. Pion gravitational form factors
The results for the combined-irrep isosinglet, nonsinglet, and gluon pion GFFs renormalized in the MS scheme at µ = 2 GeV are shown in Fig. 7, along with correlated fits to the data points using the n-pole model where α and Λ are free parameters.We set n = 1 (monopole), which we find to be the choice that yields the highest p-value for all GFF fits.We test model dependence by also using the z-expansion [80] F where t cut = 4m 2 π , α k are free parameters, and we set k max = 2.The nonsinglet GFFs are more precise than the isosinglet ones, due to the cancellation of correlated noise between the light-quark and strange disconnected contributions, and due to the fact that they do not mix with the less precise gluon contribution.The analogous results for each quark flavor, i.e., the light-quark (u + d) and strange (s) contributions, are presented in Fig. 8, while Fig. 9 includes a comparison of the gluon and different flavor quark monopole fits.The gluon GFFs are fit using the monopole and z-expansion models, while the light quark and strange GFFs and their fits are determined by taking linear combinations of the isosinglet and non-singlet components as The parameters of the monopole and z-expansion fits are shown in Table III for A π i (t) and Table IV for D π i (t).In Fig. 10, we present the renormalization scheme and scale independent total GFFs A π (t) and D π (t), obtained by sums of the corresponding isosinglet and gluon GFFs.IV.Fit parameters of the monopole and z-expansion parametrizations of the t-dependence of the pion GFFs D π i (t) renormalized in the MS scheme at scale µ = 2 GeV.We note the high χ 2 /d.o.f of the D π v (t) fits, which are due to three data points fluctuating away from the trend of the curve, as seen in the center right panel of Fig. 7.
While the joint fits discussed above have good fit quality, when alternatively solving for the flavor-singlet GFFs using the two irreps individually instead of performing a combined-irrep fit as shown in this section, we find some tension between the two irreps for the results of D π q (t).The individual-irrep results, along with a discussion of this observation, are included in Appendix B.

B. Forward limits
Our results for the flavor decomposition of the momentum fraction and D-term at scale µ = 2 GeV in the MS scheme using the monopole and z-expansion fit results of Tables III and IV are shown in Table V.A π g (0) is consistent with results obtained in a previous study of the gluon GFFs on an ensemble with m π = 450 MeV [29,30], and with a lattice extraction [33] in the forward limit with N f = 2 + 1 + 1 flavors at quark masses corresponding to the physical pion mass.Our result for A π q (0) is however smaller than the result found in Ref. [33].In contrast, we find a slightly larger contribution from gluons than  from quarks to the momentum fraction of the pion at µ = 2 GeV.The separate quark-flavor contributions, A π u+d (0) and A π s (0), are also smaller than those found in Ref. [33].A possible explanation could be that the latter were computed on an N f = 2 + 1 + 1 ensemble at the physical quark mass, while our results were obtained on a single ensemble at m π ≈ 170 MeV and could not be extrapolated to the physical point.Our results for A π q (0) are also larger than what was found in Ref. [31] after extrapolation to the continuum limit.Our results for the total momentum fraction are slightly larger than the sum rule prediction, A π (0) = 1.
The total D term obtained is consistent with χPT when the first chiral-symmetry breaking correction [16], which for the quark masses of this ensemble is estimated to result in D π (0) ≈ −0.96, is taken into consideration, and smaller in magnitude than the leadingorder chiral limit prediction of D π (0) ≈ −1.Our result for D π q (0) is statistically consistent with the result found in Refs.[27,28], which was computed from several ensembles at heavier pion masses and extrapolated to the physical point, neglecting quark disconnected contributions and mixing with the gluon EMT.We find D π g (0) to be smaller in magnitude than the result at m π ≈ 450 MeV neglecting mixing with quarks found in Ref. [30], [−0.793 (84)].
Interesting physical comparisons can be made by considering the t-dependence of the GFFs, which we find to be different between A π and D π as seen in Fig. 11, and consistent with the NLO χPT prediction [16] in the small |t| region, as seen in Fig. 10.Hadron radii associated with the spatial distributions of their physical properties have historically been defined in relation to the derivative of their form factors at t = 0.For example, the magnitude of the derivative of the pion vector form factor F π v is related to its charge radius, while that of A π (t) is related to its mass radius.We compare the relative sizes of the

NLO PT
FIG. 10.The scale-and scheme-independent total GFFs of the pion, obtained by summing the gluon and quark contributions.The red bands show the next-to-leading order (NLO) χPT prediction for the low |t| region, using a range of estimates for the low-energy constants, as presented in Ref. [16].
Using the PDG averaged value for dF π v /dt| t=0 [81] and the monopole fit of Sec.IV A, we obtain that the charge radius of the pion is approximately 1.6 times larger than its mass radius.From equivalent relations for the mass radii of the individual constituents, we find r π g,mass /r π q,mass ≈ 1.1, and r q,mass in agreement with a phenomenological extraction of the pion quark GFFs from experimental measurements [21].The equivalent quantity for D π 1 D π i (0) is, however, approximately 2.5 times smaller for i = q than what was found in Ref. [21].

V. SUMMARY AND CONCLUSION
We present the first flavor decomposition of the tdependence of the pion GFFs, calculated using a lattice QCD ensemble with quark masses yielding a close-tophysical pion mass.We constrain the gluon, isosinglet, and nonsinglet GFFs renormalized in the MS scheme at scale µ = 2 GeV, accounting for mixing between the gluon and isosinglet contributions.From our results, we perform the first extraction of renormalization scaleand scheme-independent hadron gravitational form factors from lattice QCD.From fitting the t-dependence of A π , we find indications of a smaller mass radius than charge radius for the pion.Our results are consistent with the momentum fraction sum rule, and with the NLO χPT prediction for the pion D term.We observe some tension between D π q (t) obtained by fitting each lattice operator irrep individually, as discussed in Appendix B.Even though this effect is not significant for this calculation, it is important to understand it for future extractions with higher precision.This result also highlights the value of considering different irreducible representations when extracting GFFs from lattice QCD.Another interesting feature in the analysis of the data is that excited state contamination affects the extraction of bare D π g (t) significantly more than any of the other GFFs, as discussed in Appendix A. No such effect was noted in a previous extraction of A π g (t) and D π g (t) at m π ≈ 450 MeV [30].Future improvements, besides the repetition of the calculation on additional lattice ensembles in order to enable the continuum, physical quark mass, and infinite-volume limits to be taken, could include using a variational basis of hadron interpolators [82][83][84] to better control excited state contamination.Another improvement would be the use of gauge-invariant renormalization schemes [85,86] that allow the renormalization mixing matrix of the EMT to be computed nonperturbatively on large-volume ensembles.
As few experimental constraints on the pion GFFs exist to date, the results in this work, and future calculations with further improved systematics, provide particularly valuable information on hadron structure.The structure of the pion is of particular interest as the pseudo-Goldstone boson of dynamical chiral symmetry breaking.Taken together with the first constraints on proton GFFs from experimental measurements in recent years [17,18], and the agreement with lattice QCD results for the gluon contribution [30,87], these developments mark a new milestone in our understanding of the gravitational properties of hadrons.

3 (
FIG. 3. The computed contributions RRI v C RI/MS v (a), RRI qq C RI/MS qq (b), and RRI qq C RI/MS qg (c) are shown as empty markers for irreps τ (3) 1 (blue) and τ (6) 3 (orange).Fits to various sets of points are performed as described in Sec.III C. The model-averaged fits to the data are shown as bands, and the resulting estimates of the renormalization contributions R RI v C RI/MS v , R RI qq C RI/MS qq , and

FIG. 4 .
FIG. 4. The computed RRI gq C RI/MS qq (a) and RRI gq C RI/MS qg (b), and their corresponding fits, performed as discussed in Sec.III D.

FIG. 5 .
FIG. 5. Data for RRI qg C RI/MS gg (a) and RRI qg C RI/MS gq (b) with unflowed external gluon fields are shown with empty markers, while the flowed ones are shown as empty diamonds.The two sets of data are fit simultaneously as described in Sec.III E, and the extracted values for R RI qg C RI/MS gg and R RI qg C RI/MS gqfor each of the two irreps are shown as filled markers at a 2 p2 = 0. We note

2
is included in the model, as discussed in Sec.III B. We extract R RI qg C RI/MS gq from model averages over combined fits to the flowed and unflowed RRI qg C RI/MS gq using the same fit ranges as for R RI qg C RI/MS gg .The data and averaged fits are shown in Fig. 5.

FIG. 6 .
FIG. 6.Data for RRI gg C RI/MS gq (a) and RRI gg C RI/MS gg (b), and corresponding fits performed as described in Sec.III F. The notation E, the flowed version of RRI gg C RI/MS gq demonstrates primarily logarithmic dependence and is fit using P log 2 of Eq. (45), while the unflowed version has strong discretization effects and is modeled using P log 2 +P inv 4 , as discussed in Sec.III B. R RI gg C RI/MS gq is obtained by combined fits to both versions, with the same momentum democracy cut and fit ranges as in Sec.III E. The results are shown in Fig. 6a.As discussed in Sec.III B, RRI gg C RI/MS gg qq = 1.039(28),R MS qg = 0.081(19) , R MS gq = −0.180(23),R MS gg = 1.625(48) .

FIG. 7 .FIG. 8 .
FIG.7.The isosinglet (top), nonsinglet (center), and gluon (bottom) GFFs of the pion renormalized in the MS scheme at scale µ = 2 GeV.The A π i (t) form factors are shown on the left and the D π i (t) on the right.Fits using the monopole model of Eq. (64) and the z-expansion of Eq. (65) are shown, with fit parameters collected in Tables III and IV.
mass radius and the charge radius by taking the ratior π mass r π EM ∼ dA π /dt| t=0 dF π v /dt| t=0.

FIG. 11 .
FIG.11.The pion GFFs rescaled by their central value at t = 0, as obtained from the monopole model with fit parameters collected in Tables III and IV.

FIG. 15 .
FIG.15.Examples of gluon effective GFFs, computed from Eq. (A2) using summed ratios with τcut = 4.The effective D π,eff g , shown in the bottom panels of (a) and (b), does not reach a plateau until later ts values compared to A, a feature that is particularly prominent for smaller t-bins like those shown in (a).The bands are not fit to the data shown but correspond to bare GFFs obtained by fitting the bare matrix elements used to compute the results from the main text.

TABLE I .
Specifics of the lattice ensembles used in this work.Ensemble A, generated by the JLab/LANL/MIT/WM groups

TABLE II .
Number of sources Ns for which the connected quark contribution to the three-point function, Eq. (21), is computed for each sink time ts.
in two ways: through RRI qq C

TABLE V .
The flavor decomposition of the momentum fraction and the D term of the pion, obtained from the monopole and z-expansion fits to the pion GFFs, renormalized at µ = 2 GeV in the MS scheme, with parameters shown in Tables III and IV.