Gluon gravitational structure of hadrons of different spin

The gravitational form factors (GFFs) of hadrons encode the matrix elements of the energy momentum tensor of QCD. These quantities describe how energy, spin, and various mechanical properties of hadrons are carried by their quark and gluon constituents. We present the gluon GFFs of the pion, nucleon, $\rho$ meson, and $\Delta$ baryon as functions of the squared momentum transfer $t$ in the region $0 \leq -t<2 \; \text{GeV}^2$, as determined in a lattice QCD study with pion mass $m_{\pi} = 450(5) \; \text{MeV}$. By fitting the extracted GFFs using multipole and z-parameter expansion functional forms, we extract various gluon contributions to the energy, pressure, and shear force distributions of the hadrons in the 3D and 2D Breit frames as well as in the infinite momentum frame. We also obtain estimates for the corresponding gluon mechanical and mass radii, as well as the forward-limit gluon contributions to the momentum fraction and angular momentum of the hadrons.


I. INTRODUCTION
Understanding the internal dynamics of hadrons in terms of their fundamental quark and gluon constituents has been a goal of particle and nuclear physics since the first experimental probe of proton substructure at SLAC [1] and the subsequent development of the theory of quantum chromodynamics (QCD) [2][3][4]. However, many aspects of hadron structure have not yet been fully constrained from theory or experiment, including the gravitational form factors (GFFs) [5] of hadrons, defined from matrix elements of the QCD energy-momentum tensor (EMT). These form factors describe how energy, spin, pressure, and shear forces are distributed within hadrons [6]; therefore, their determination is of fundamental significance.
The off-forward hadron matrix elements of the symmetric 1  (1) where |h(p, s) denotes a hadronic state with fourmomentum p and polarization s, and K {µν},h,j ss are kinematic coefficients symmetrized over their Lorentz indices as a {µ b ν} = (a µ b ν + b µ a ν )/2, written in terms of P = (p + p )/2 and ∆ = p − p. The GFFs G h,j i (t) are functions of the Mandelstam variable t = ∆ 2 , and j indexes the different GFFs in the decomposition for hadron h. An analogous decomposition of the total conserved EMT T µν = q T µν q + T µν g yields the total GFFs G h,j (t) = q G h,j q (t) + G h,j g (t).
The GFFs associated with the symmetric traceless part of the EMT correspond to the second Mellin moments of the corresponding generalized parton distributions (GPDs) [8][9][10], which allows them to be constrained by experimental data from deeply virtual Compton scattering (DVCS) [11][12][13] and meson production [14,15]. For example, data from the Belle experiment at KEKB [16,17] has been used to constrain the pion quark GFFs [18], while the nucleon quark GFFs have been studied from DVCS with the CLAS detector [19][20][21][22] at the Thomas Jefferson National Accelerator Facility (JLab). The PANDA experiment at the Facility for Antiproton and Ion Research (FAIR) [23], as well as future experiments at SuperKEKB, the International Linear Collider (ILC), the Japan proton accelerator complex (J-PARC) [24] and the nuclotron-based ion collider facility (NICA) [25] will further constrain the quark GPDs and thus GFFs of various hadrons.
The gluon contributions to the GFFs, on the other hand, are far less well constrained and have so far only been studied in an extended holographic light-front QCD framework for the pion and nucleon [59], in lattice QCD calculations for the pion [60], nucleon [52,56,[60][61][62][63], and φ meson [64], in almost all cases at larger-thanphysical values of the quark masses. While no exper-  [58].
imental constraints on the gluon GFFs of any hadron have been achieved to date, the gluon GFFs of the nucleon are accessible via photo-or leptoproduction of J/ψ and Υ [65][66][67]; J/ψ production is studied in experiments that are ongoing at JLab [68], while Υ production studies are planned at the electron-ion collider (EIC) [69]. Improved QCD constraints on the gluon GFFs of the nucleon and other hadrons are particularly valuable at the current time as they can inform the target kinematics for these experiments and provide theory predictions to test against future experimental results.
In this work, we present a lattice QCD calculation of gluon GFFs of the pion and nucleon at unphysically heavy quark masses, incorporating additional data corresponding to spin-nonconserving channels and an improved statistical analysis compared with the previous study of Ref. [60]. We further undertake a first study of the complete set of gluon GFFs of the ρ meson and ∆ baryon, which are stable at these quark masses, to investigate the gluon radii and gluon energy, pressure, and shear force distributions of hadrons of higher spin. In Sec. II we outline the lattice QCD calculation and analysis, discussed more extensively in Appendix A, and show the extracted renormalized GFFs for all hadrons considered. In Sec. III, we present our results for the radii and densities in different frames. In Sec. IV, we provide a summary and outlook.

II. GRAVITATIONAL FORM FACTORS FROM LATTICE QCD
In this section we discuss the decompositions of hadronic matrix elements of the gluon EMT into gluon GFFs for the pion, nucleon, ρ meson, and ∆ baryon, and present the results of our lattice QCD extraction of these quantities.
We use a single ensemble of 2820 configurations of lattice volume 32 3 × 96 with N f = 2 + 1 quark flavors, with a heavier-thanphysical pion mass of m π = 450(5) MeV and lattice spacing a = 0.1167 (16) fm [70]. The ensemble was generated using the Lüscher-Weisz gauge action [71] and cloverimproved Wilson quarks [72] with clover coefficient set to the tree-level tadpole-improved value and constructed using stout-smeared links [73]. The specifics of the ensemble are summarized in Table I [58,74]. Our results for the nucleon and pion GFFs are consistent with but more precise than those of Ref. [60], which studied those states on a subset of the data used in this work, including only spin-conserving channels.
Our methods are similar to those of Ref. [60], but with an improved statistical analysis and with the necessary extensions to treat hadrons of higher spin. The calculation proceeds independently but analogously for each hadron, in several stages detailed in Appendix A: (i) Compute (hadron-independent) measurements of the symmetric traceless gluon EMT, discretized and projected to irreducible representations (irreps) of the hypercubic group that are protected from mixing with lower-dimensional operators (Appendix A 1).
(ii) Compute hadronic two-point correlation functions (Appendix A 2) and three-point correlation functions (Appendix A 3) including insertions of the gluon EMT.
(iii) Extract hadronic matrix elements of the gluon EMT by fitting ratios of three-and two-point correlation functions (Appendix A 4).
(iv) Extract the renormalized gluon GFFs by fitting the constraints defined by the measured matrix elements and Eq. (1) (Appendixes A 5, A 6), incorporating the renormalization factors for the different irreps computed using a nonperturbative RI-MOM prescription and a one-loop perturbative matching to MS at µ = 2 GeV. We neglect effects due to mixing with the quark GFFs under renormalization, which are expected to be O(10%) [56,75] and thus small compared with the statistical uncertainty of the calculation.
Our lattice calculation yields the GFFs at a set of discrete values of the squared momentum transfer t but, as discussed in Sec. III, subsequent extrapolations to the forward limit and derivations of densities and radii require models of the t dependence of the GFFs. We consider two different ansätze, a multipole and a modification of the z-parameter expansion or "z-expansion" [76,77]. The multipole form is defined as where α and Λ are free parameters and we set n = 3 (tripole) in order for all the integrals that define the energy, pressure, and shear force densities discussed in Sec. III to converge. As introduced in Ref. [77], the modified z-expansion we consider is FIG. 1. A π g (t) and D π g (t) renormalized at µ = 2 GeV in the MS scheme. The bands correspond to the multipole form (Eq. (2)) with n = 3 and the modified z-expansion (Eq. (3)) with kmax = 2, with fit parameters shown in Tab. II. where α k are free parameters, and Λ is constrained as described below. This functional form can be interpreted as a series of corrections to the envelope defined by the multipole form, which coincides with Eq. (3) when k max = 0. The multipole envelope is necessary for convergence of the density integrals discussed in Sec. III. Following Ref. [77], we set k max = 2, t 0 = t cut (1 − 1 + (2 GeV) 2 /t cut ), and t cut = 4m 2 π , using m π = 450 MeV. For each GFF, we use Λ obtained from the multipole fit to the same GFF as a prior for the parameter Λ in the z-expansion, retaining correlations between the prior and the data, 2 thereby reducing the number of free parameters to prevent overfitting and explicitly enforcing the notion of the modified z-expansion as a correction to the multipole envelope. As detailed in Appendix A, we fit the models to bare GFFs and renormalize afterwards to circumvent the d'Agostini bias [79]; in this section we present the resulting renormalized parameters α and α k .

A. Pion
The pion matrix element of the symmetric gluon or quark contribution to the energy-momentum tensor can be decomposed as 2 Compare with the discussion of "chained fitting" in Ref. [78]. (3)] with kmax = 2 and n = 3 models for the t dependence of the renormalized pion GFFs. For all GFFs, the parameter Λ of the z-expansion fit is consistent with the prior provided by the tripole fit and is thus not shown. The parameters α and α k are renormalized at µ = 2 GeV after fitting the bare GFFs as described in Appendix A 5.
where i ∈ {q, g}, g µν is the Minkowski space-time metric, and we have defined the traceless piece O µν (π) for later convenience. A π i (0) is the traceless contribution to the momentum fraction carried by the quarks or gluons and must satisfy A π (0) = i A π i (0) = 1 because of Poincaré symmetry. D π i (t) is related to the mechanical properties of the pion. In the forward and chiral limits, the total D π (0), also called the D-term or Druck term, is predicted to be −1 up to chiral-symmetry breaking effects [26,80,81].c π i (t) appears due to the nonconservation of the separate quark and gluon contributions and vanishes for the total EMT, i.e.c π (t) =c π g (t) + qc π q (t) = 0. Our results for the two renormalized traceless gluon GFFs of the pion, A π g (t) and D π g (t), are shown in Fig. 1. The fit parameters for the two ansätze, Eqs. (2) and (3), are shown in Table II, and the predicted forward-limit gluon momentum fraction and D-term are shown in Table III. We note that the sum of our gluon D-term with tripole z-expansion A π g (0) 0.537 (45) 0.544(46) D π g (0) −0.793(84) −0.74 (21) TABLE III. The forward-limit values of the gluon momentum fraction and the gluon D-term, obtained from the tripole and modified z-expansion fits to the pion GFFs, renormalized at µ = 2 GeV in the MS scheme, with parameters shown in Table II. the value D π u+d (0) = −0.264(32) (extrapolated to the physical pion mass) from Ref. [51] is statistically consistent with the chiral prediction, although this may be a coincidence that does not survive a chiral and continuum limit extrapolation.
Our results for the three renormalized traceless gluon GFFs of the nucleon are shown in Fig. 2, the tripole and modified z-expansion fit parameters are shown in Table IV, and the predictions for the forward-limit gluon contributions to the momentum fraction, D-term, and angular momentum are shown in Table V. Note that the increase order-by-order of the parameters in the modified z-expansion fit for the D-term are not of concern, as we consider here a modification of the standard z-expansion  (3)] with kmax = 2 and n = 3 models for the t dependence of the renormalized nucleon GFFs. For all GFFs, the parameter Λ of the z-expansion fit is consistent with the prior provided by the tripole fit and is thus not shown. The parameters α and α k are renormalized at µ = 2 GeV after fitting the bare GFFs as described in (57) TABLE V. The forward-limit values of the gluon contributions to the momentum fraction, D-term, and angular momentum of the nucleon, obtained from the tripole and modified z-expansion fits of the nucleon GFFs renormalized at µ = 2 GeV in the MS scheme, with parameters shown in Table IV. for which the guarantees of convergence for the standard form do not apply. These results are based on a superset of the data presented in Ref. [60], including a larger number of sources per configuration and all four nucleon polarization channels, rather than only the two spin-conserving ones. The additional data allows a nonzero functional fit for B N g (t) to be resolved. Additionally, the behavior of A N g (t) and D N g (t) in the lower half of the −t region studied is shifted slightly compared to what was found in Ref. [60]. These updated results show a qualitative difference between the low-|t| behavior of the model obtained for D N g (t) with the tripole functional form (that is by definition monotonic), and the modified z-expansion fit (that is allowed to be nonmonotonic). In general, monotonically increasing behavior is universally expected for the total D-term of any stable mechanical system [90]; however, no such prediction exists for the individual quark and gluon contributions. It is not clear whether the suppression of the lowest-|t| point of D N g (t) is physical or due to a statistical fluctuation or unquantified systematic uncertainty. If physical, it would imply a qualitatively different t dependence of the gluon GFFs compared with that typi-0.00 0.25 0.50 0. 75 Table IV. cally assumed for the quark GFFs, and suggest that the multipole functional form often used by default for these quantities is not a good model at low |t|. If the first data point is considered an outlier and excluded from the z-expansion fit, the resulting error bands do not exclude nonmonotonicity but encompass both monotonic and nonmonotonic forms.
Our predictions for A N g (0) and J N g (0) are statistically consistent with the equivalent quantities found in a lattice QCD calculation at quark masses corresponding to the physical value of the pion mass [56].

C. ρ Meson
Following the conventions of Ref. [91], the matrix elements of the gluon or quark contribution to the symmetric EMT for the ρ meson can be decomposed as where m ρ denotes the mass of the ρ meson, and α (p, λ) is the polarization 4-vector for a massive spin-1 particle, which satisfies αβ (p) (9) with λ ∈ {1, 0, −1}. Note that the subscript ρ in m ρ and the superscript ρ in the GFFs such as A ρ 0,i (t) is a label for the ρ meson and not a Lorentz index.
The momentum sum rule constrains the total momentum fraction to be A ρ (0) ≡ A ρ 0 (0) = 1, and the forwardlimit angular momentum must be equal to the spin of the hadron, i.e., J ρ (0) = 1. The interpretation of the D-term of the ρ meson is more complicated than those of the pion and the nucleon, since there are three such terms [48,91], one of monopole and two of quadrupole order, corresponding to frame-dependent linear combinations of D ρ 0 (0), D ρ 1 (0), and E ρ (0). Here we focus on the forward limit of the form factor that is the coefficient of the same Lorentz structure corresponding to the D-term GFFs of the nucleon and pion, namely D ρ (0) ≡ −D ρ 0 (0), which is unconstrained from theory. There are two GFFs that arise from the trace of the EMT,c 1,ρ i (t) andc 2,ρ i (t), and their contribution to the total EMT must be equal to zero. In contrast to the pion and the nucleon GFF decompositions, the traceless piece of the EMT matrix element for hadrons of spin-1 gives rise to a nonconserved GFF,f ρ i (t), which we can access in our calculation and vanishes when summed over the quark and gluon contributions.
Our results for the seven renormalized traceless gluon GFFs are shown in Fig. 3. Model fit parameters are tabulated in Table VI, excluding for the GFF D ρ 1,g (t), which is not well described by either model ansatz. The conserved gluon predictions of these quantities are presented in Table VII.
We find that approximately half of the angular momentum of the ρ meson is carried by gluons. Interestingly, the NJL model [33] predicts that half is carried by the quark spin. Just as in the nucleon case, we find a significant difference between the forward limit of the D form factor resulting from the tripole fit and the z-expansion (see Table VI) which can be traced to the difference of the two fits in the low momentum region of D ρ 0,g (t) as seen in Fig.  3(c). Here the first data point again suggests nonmonotonic behavior that cannot be captured by a multipole form, but it is unclear whether this is a physical effect or due to a statistical fluctuation or underestimated systematic uncertainty. In Ref. [48], the total D ρ 0 (t) from the light-front constituent quark model was found to be nonmonotonic, in contrast with the prediction from the NJL model [33].
The total momentum fraction A ∆ (0) ≡ F ∆ 10 (0) is constrained to equal 1. As with the ρ, there are three total D-terms of different order and we again focus our discussion on the form factor that is the coefficient of the same Lorentz structure as the nucleon and pion D-terms, namely D ∆ (0) ≡ F ∆ 20 (0). The total forward-limit angular momentum, J ∆ (0) = F ∆ 40 (0), is constrained to be equal to 3/2 [37]. Just as in the case of the ρ meson, there are nonconserved GFFs related to the trace piece, F ∆,i 30 (t) and F ∆,i 31 (t), that we do not have access to in this calculation, and one nonconserved GFF F ∆,i 60 (t), that arises from the traceless EMT and that we are able to constrain.    Table VI. No fit is shown for D g 1 (t) because the corresponding data is not well described by either functional form.       Table VIII. No fit is shown for the four GFFs that are not resolved from zero.  (3)] with kmax = 2 and n = 3 models for the t dependence of the renormalized ρ GFFs. The signal for D ρ 1,g (t) is not well described by either functional form and therefore no fit is shown. For all GFFs that have been fit, the parameter Λ of the z-expansion fit is consistent with the prior provided by the tripole fit and is thus not shown. The parameters α and α k are renormalized at µ = 2 GeV after fitting the bare GFFs as described in Appendix A 5. Our results for the eight renormalized traceless gluon GFFs of the ∆ baryon are shown in Fig. 4. Only four of them are resolved from zero, and their fit parameters are shown in Table VIII. The conserved gluon contributions to the forward limit quantities obtained from the tripole and z-expansion fits are shown in Table IX.

III. DENSITIES AND RADII FROM GFFS
In the decomposition of the matrix element h(p , s )| T µν i |h(p, s) , the GFFs are Lorentz scalars but their coefficients depend on the reference frame. The spatial energy, pressure, and shear force densities of hadrons are related to Fourier transforms of the   (3)] with kmax = 2 and n = 3 models for the t dependence of the renormalized ∆ GFFs. Only four of the eight GFFs are fit, as the others are not resolved from zero as seen in Fig. 4. For all GFFs that have been fit, the parameter Λ of the z-expansion fit is consistent with the prior provided by the tripole fit and is thus not shown. The parameters α and α k are renormalized at µ = 2 GeV after fitting the bare GFFs as described in (18) TABLE IX. The forward-limit values of the conserved gluon contribution to the momentum fraction, D form factor, and angular momentum of the ∆ baryon, obtained from the tripole and modified z-expansion fits of the ∆ GFFs, renormalized at µ = 2 GeV in the MS scheme, with parameters shown in Table VIII. momentum space matrix elements, and therefore are also frame dependent. In this work we consider these densities in two frames, namely the Breit frame and the infinite momentum frame.
Breit frame (BF).-The "brick-wall" frame in which there is no energy transfer to the system, ∆ 0 = 0, and additionally P = 0. This is the frame traditionally used to define spatial distributions, such as the charge distribution in terms of the electromagnetic form factors [92]. The equivalent 3D density for the EMT in the BF (the BF3 density) is It is known that the Fourier transform of a form factor in the BF is not a relativistically correct way to define the corresponding spatial distributions [93,94] since one is free to multiply the distribution by a boost factor that cannot be uniquely defined in relativistic quantum field theory. However, following the phase-space approach introduced in Refs. [95,96], Eqs. (12) and (13) can be interpreted as quasidistributions instead of densities, and there is no ambiguity with respect to the boost factor.
Infinite momentum frame (IMF).-The elastic frame in which ∆ · P = 0 and P z → ∞. In this frame there is Galilean symmetry in the transverse plane and thus 2D Fourier transforms of the momentum tensor matrix elements can be interpreted as spatial densities. 3 The expression for the EMT density in this frame is We note that for the case of the pressure and shear forces of spherically symmetric hadrons, it was recently shown that the densities in the two frames are related by Abel transformations [101].

A. Pion
The expressions for the various EMT distributions of the pion in terms of the corresponding GFFs are listed in Appendix B 1. In Fig. 5, we present our results for the gluon contribution to the energy density ε(r), the pressure p(r), and the shear forces s(r) in the 3D and 2D BF, and in the IMF. The definitions of the energy and pressure densities of the individual constituents include the nonconserved GFFc(t), which cancels between the quark and gluon contributions. Since we cannot constrain this term from our calculations, the results shown in Fig. 5, and for the rest of the hadrons in the sections below, include only the traceless gluon contribution to the densities. 3 Another frame that can be considered is the front-form Drell-Yan frame, in which ∆ − = 0 and ∆ + = 0. 2D Fourier transforms in the Drell-Yan frame can be correctly interpreted as spatial distributions [33,95,[97][98][99], since in the light-cone transverse boosts are Galilean [100]. We choose not to discuss this frame here, since the energy density corresponds to a different component of the energy momentum tensor and is thus not directly comparable with the instant-form energy density. The pressure and shear forces for the pion and the nucleon are identical in the infinite momentum frame and the Drell-Yan frame.
From the pressure density, it is interesting to test whether the 3D and 2D von Laue stability conditions [102] for the total pressure of a composite particle ∞ 0 dr r 2 p BF3 (r) = 0, ∞ 0 dr ⊥ r ⊥ p BF2/IMF (r ⊥ ) = 0 (15) hold for the traceless gluon piece alone. Indeed, by numerical integration we find that the pressures are consistent with the von Laue condition in all frames and using both multipole and z-expansion functional forms to model the t dependence of the GFFs. Another important stability condition first shown in Ref. [88] and recently extended in Ref. [98] is that for the total D-term which is satisfied by the gluon contribution to the pion D-term in Table II. We also find that the hadron stability conditions [88,99] hold for the traceless gluon piece of the pion pressure and shear force, which allows us to define a partial gluon mechanical mean square radius for the pion r 2 π,g [see Eq. (B5)]. Our results for the mechanical radii as well as the mass radii, defined in Appendix B 1 as appropriate averages of the distance from the center of the hadron weighted by the energy density, are presented in Table X.

B. Nucleon
The BF and IMF densities of the nucleon EMT have been studied previously [19,77,88,95] and are defined in Appendix B 2. We note that the 2D Breit frame pressure and shear force coincide with their IMF equivalents. Our results for the symmetric traceless gluon contributions to the densities are shown in Fig. 6. The difference between the results based on tripole and z-expansion fits to the GFFs is due to the difference between the two fits of D N g (t) in the low −t region, as discussed in Sec. II B. A nonmonotonic gluon D N g (t) causes the traceless gluon pressure to have two nodes, which is different than the form of the quark pressure distribution as found in Ref. [19]. Future lattice, experimental, and phenomenological extractions of D N g (t) and D N q (t) will help to clarify the picture.
From the nucleon density results, we find that the pressures in all frames and models are consistent with the von Laue condition (15); however, Eq. (17) is only satisfied for the tripole fit. Moreover, the D-term fit by the zexpansion is not strictly positive within uncertainty, and thus we only present the mechanical radius of Eq. (B5) for the tripole fit. The mass radii definitions are shown in Appendix B 2, and the corresponding numerical results of all radii are presented in Table X Table IV. The expressions for the pressure and shear force are identical in the 2D BF and the IMF.

C. ρ Meson
Beyond the lowest-order energy, pressure, and shear force densities of the pion and the nucleon, the structure of hadrons of spin = 1 or higher depends on additional quadrupole densities. The BF3 densities and mass radii of the ρ meson were derived in Refs. [48,91] and are listed in Appendix B 3. We also derive expressions for the lowest-order BF2 and IMF distributions.
Our numerical results for the lowest-order (Fig. 7) and higher-order densities (Fig. 8) are partial and exclude terms depending on D ρ,g 1 (t), since its signal is not well modeled by the ansätze considered here, as discussed in Sec. II C. We note that almost all of the quadrupole densities are poorly constrained, and therefore we only consider the mechanical stability conditions for the lowestorder densities. In particular, Eq. (17) only holds for the BF3 and BF2 from the tripole GFFs fits. The corresponding mechanical radii of Eq. (B5) and the mass radii are shown in Table X.

D. ∆ Baryon
The BF3 distributions for the ∆ baryon were first presented in Ref. [103] and extended in Ref. [36], and are listed in Appendix B 4, along with our derived lowest-order BF2 and IMF contributions to the densities. The higher-order densities are not well constrained, and we limit our mechanical stability check to the lowest-order densities only. We find that Eq. (17) only holds for the tripole fits to the GFFs. Our estimates for the accessible gluon mechanical and mass radii are shown in Table X. As discussed in Sec. II D, we were able to model only four of the ∆ gluon GFFs with the tripole and zexpansion fits, and therefore our results for the gluon densities shown in Figs. 9 and 10 and the radii in Table X are partial.

IV. SUMMARY AND CONCLUSION
In this work, we present a lattice QCD extraction of the gluon GFFs of the pion, nucleon, ρ meson, and ∆ baryon at quark masses corresponding to a pion mass m π = 450(5) MeV in the range 0 ≤ −t < 2 GeV 2 . All of the pion and nucleon GFFs, along with six of seven of the ρ and four of eight of the ∆ GFFs, are fit using the multipole and modified z-expansion models of Eqs. (2) and (3). Most of the GFF fits are consistent between the two models considered, with the most significant exception being the nucleon D N g (t) which, while consistent within error, shows a qualitative difference in behavior between the (monotonic) tripole fit and (nonmonotonic) z-expansion fit. Further calculations with different ensembles are needed in order to determine whether this nonmonotonicity is physical, or the result of poorly quantified systematic uncertainties or a statistical fluctuation in the data. This inconsistency, however, brings to attention the importance of considering a variety of models when studying the GFFs of hadrons at different values of the energy transfer.
In Figs. 12, 13, and 14 we provide a summary of the gluon contributions to the momentum fraction form factor A h (t), angular momentum form factor J h (t), and the D h (t) form factor. In all cases, the meson GFFs fall off more slowly as functions of −t than the baryon GFFs. It will be important for future work to study them at higher magnitudes of energy transfer in order to fully quantify their behavior. In Fig. 11, we summarize the forward limit results for the gluon momentum and spin fractions, with the gluon spin fraction defined as the ratio J h g (0)/J h (0) of the gluon contribution to the forward-limit angular momentum J h g (0) and the total spins J h (0) of the corresponding hadrons. The gluon momentum fraction is larger for the mesons and smaller for the baryons, decreasing with increasing hadron mass.
Additionally, for each model of the t dependence of the GFFs we compute the gluon contributions to the hadron energy, pressure, and shear force spatial densities and mass radii, as well as their mechanical radii given stability conditions are satisfied. We consider both the Breit frame and the infinite momentum frame of the densities and the radii. The estimates made using each model are consistent with each other, with the exception of a dis- Tripole Fit BF3 r 2 p ,g 0 (r) BF2 r p ,g 0 (r ) IMF r p ,g 0 (r ) Z-Expansion Fit BF3 r 2 p ,g 0 (r) BF2 r p ,g 0 (r ) IMF r p ,g 0 (r ) Tripole Fit BF3 r 2 s ,g 0 (r) BF2 r s ,g 0 (r ) IMF r s ,g 0 (r )  Z-Expansion Fit BF3 r 2 ,g 0 (r) BF2 r ,g 0 (r ) IMF r ,g 0 (r )    Tables III, V, VII, and IX crepancy in the nucleon gluon D-term and energy, pressure, and shear force densities between the two ansätze, that can be traced to the nonmonotonicity of the D N g (t) z-expansion fit. 4 For the aforementioned hadronic quantities that additionally depend on trace or antisymmetric GFFs, or on the ρ and ∆ GFFs that are not consistent with our fit models, we provide partial contributions from only the GFFs that we have constrained. In order to better understand in what ways the gluons and quarks separately contribute to the gravitational structure of hadrons, it will be important to constrain all of the GFFs in future studies. Moreover, our renormalization procedure for the gluon EMT does not include mixing with quark operators. Due to the small magnitude of the mixing renormalization coefficient [56,75], the effect is expected to contribute at the level of a few percent, which is negligible compared to the current statistical uncertainties of 4 The IMF energy density does not depend on D N g (t) and is consistent between the tripole and the modified z-expansion fit. This study uses heavier-than-physical quark masses at a single lattice spacing and volume, and therefore our results are subject to unquantified systematic uncertainties that need to be addressed in future studies. For the pion and nucleon, repeating the calculation of the gluon GFFs using different ensembles is critical in order to control the effect of systematic uncertainties for comparisons with future experimental data from J/ψ and Υ production processes. For unstable hadrons like ρ and ∆, lattice QCD methods are the only known way to access their gluon GFFs; studying them at lighter quark masses, where they are not stable, will require more computationally involved Lüscher method analyses [104][105][106][107].   Tables II, IV, VI and VIII. sources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357, and the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI-1548562. Computations were carried out in part on facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy. The authors thank Robert Edwards, Balint Joo, Kostas Orginos, and the NPLQCD Collaboration for generating the ensembles used in this study. The Chroma [108], QLua [109], QUDA [110][111][112], QDP-JIT [113], and QPhiX [114] software libraries were used in this work. Data analysis used NumPy [115], SciPy [116], pandas [117,118], lsqfit [119], and gvar [120]. Figures were produced using matplotlib [121], seaborn [122], and Mathematica [123].
Appendix A: NUMERICAL AND ANALYSIS DETAILS

Symmetric traceless gluon EMT and renormalization
The symmetric traceless part of the EMT,T µν , can be obtained via the Belinfante-Rosenfeld process 5 [125] and can be decomposed into gluon and quark terms,T µν =T µν g + qT µν q , wherê where F µν is the gluon field-strength tensor of QCD, D µ is the covariant derivative, , γ µ are the Dirac matrices, the repeated indices are contracted with the Minkowski space-time metric g µν , and the trace is over the color indices.
The gluon field-strength tensor F µν can be defined up to finite lattice spacing corrections on a Euclidean spacetime lattice as where g 0 is the bare gauge coupling, the label E denotes Euclidean spacetime, and P µν (x) is the clover term defined in terms of gauge links U µ (x) as The momentum-projected traceless symmetric piece of the gluon EMT in Euclidean space can be defined aŝ where the repeated indices are contracted with the Euclidean metric δ µν (i.e. the Kronecker delta in four dimensions) and the trace is over the color indices. In continuous spacetime, a traceless symmetric tensor transforms in the (1, 1) representation of the Lorentz group. However, Lorentz symmetry is reduced to hypercubic symmetry on the lattice and therefore lattice operators transform in irreps of the symmetry group H(4). There are two choices of irreps that are safe from power-divergent mixing with lower-dimensional operators, namely τ and τ (6) 3 [126]. We thus project all operator measurements ofT gE µν to particular bases for these two irreps, We compute all operators in the τ irreps for all lattice momenta satisfying |∆| 2 ≤ 18(2π/L) 2 . To suppress gauge noise, we improve the operators by constructing them from gauge links that have been subjected to Wilson flow [127][128][129] to flow time t/a 2 = 1 (with integrator step size = 0.01). Different choices of flow time, as well as use of hypercubic smearing instead of Wilson flow, have been shown to give consistent results [64,130].
At finite lattice spacing, the two irreps τ renormalize differently and only coincide in the continuum limit, but all operators within each irrep share the same renormalization factor by symmetry. Ref. [60] carried out a nonperturbative RI-MOM calculation [131,132] of the renormalization factors of these operators on a smaller-volume ensemble for the same parameters as those used in this work. With a one-loop perturbative matching to the MS scheme [133], this yielded the renormalization coefficients  [60]. The uncertainties on these quantities are dominantly systematic and, because they were computed on a different ensemble, uncorrelated with the rest of the data. We choose to model their distribution as uncorrelated Gaussians. As in Ref. [60], we neglect mixing with the quark operators under renormalization, which is expected to contribute at the few-percent level [56,75].

Two-point correlation functions
As detailed in the main text, we use a single lattice ensemble in this work; parameters for this ensemble are listed in Table I. Using matching valence and sea quark actions, we compute two-point functions for varying numbers of light-quark sources on each configuration, using an average of 235 randomly chosen locations (240 for 80% of the configurations, 200 for 90%). As described below, our analysis accounts for differing numbers of sources by weighting configurations proportionately when drawing bootstrap ensembles. For each source position, we invert from a smeared source (S) and construct propagators for both a point (P) and smeared sink (S), with matching source and sink smearing for the SS propagators, using APE smearing [134] with 35 steps of gauge-invariant Gaussian smearing with width ρ = 4.7. From each propagator we construct two-point correlation functions for each hadron using the interpolating operators where all gamma matrices are Euclidean, C is the charge conjugation matrix, and color and spinor indices are left implicit.
The momentum projected two-point correlation function of the pion can be expressed as where E π p is the energy of the lowest-lying state with momentum p, andZ p = Z p when the source and sink are smeared differently. The two-point correlation function of the nucleon for spin channel σ → σ is traces are over Dirac indices, and Γ σ σ is a 2 × 2 block matrix that projects the four different spin channels of the nucleon [135], i.e.
where P + ≡ 1 2 (1 + γ t ) is a positive-energy projector. We use all four possible channels σ, σ ∈ {+1/2, −1/2}, adding significant additional data over the analysis in Ref. [60] where only the two spin-conserving channels were used. The two-point correlation function of the ρ meson, in the spherical basis of one of the 9 spin channels a → a , can be expressed as (9)]. Finally, we compute the two-point correlator of the ∆ baryon for the 10 spin channels ξ → ξ where ξ ≥ ξ , where repeated indices are summed over, Λ and D ξ σ,a = 0 for all other choices of {ξ, σ, a}. 7 We average over sources to obtain per-configuration measurements C h,2pt ss (p, t f = t − t 0 ) for each hadron h, weighting this average by the number of sources on each configuration when forming bootstrap ensembles as discussed below. The effective mass for each hadron h is defined as and constructed from the spin-averaged (over diagonal spin channels for states with spin = 0) two-point functions.
The results for each hadron are shown in Fig. 15, along with the numerical values that we use for the hadron masses m h throughout this work, which are obtained via single-state correlated fits to Eq. (A16), in regions in which the excited-state contamination is smaller than the statistical uncertainties of the effective mass function. The numerical values are given in Table XI.  Table XI.

Three-point correlation functions and ratios
We construct hadronic three-point functions of the gluon EMT operator, which are defined as where repeated indices are summed over, (x 0 , t 0 ) and (x, t ) are the source and sink positions, (y, τ ) is the operator insertion position, p and p are the three-momenta of the hadron at the source and sink, ∆ = p − p is the threemomentum injected by the operator, andT gE R,i is a gluon EMT operator projected to irrep R and basis element i as defined in Eq. (A5). These three-point functions are entirely disconnected and so may be computed by correlating two-point functions with measurements of the gluon EMT, i.e. by computing whereT gE R,i (∆, τ ) is the gluon EMT operator projected to momentum ∆ as in Eq. (A4). We average measurements of the three-point correlation functions over sources (translating appropriately) to obtain per-configuration measurements, denoted by C h,3pt ss ;R,i (p, p ; t f = t − t 0 , τ ). Given per-configuration measurements of the two-and three-point functions, we draw 1000 bootstrap ensembles, weighting the probability of drawing each configuration by the number of sources measured on that configuration. To improve signal-to-noise as discussed in Ref. [64], we perform a vacuum subtraction of each three-point correlation function within each bootstrap ensemble, where . . . indicates an ensemble average and the explicit sum is an average over sources. We form ratios of two-and three-point functions to isolate the matrix elements of interest. For all hadrons the appropriate ratio is the same, and is constructed as within each bootstrap ensemble. We have suppressed dependence on the source and sink smearing, but we carry out this computation separately using correlation functions constructed from SS-and SP-smeared propagators, yielding separate SS and SP measurements of each ratio.

Coefficients, binning, and ratio fits
The ratio in Eq. (A23) is chosen such that the leading-order t f , τ dependence and the overlap factors between the hadronic ground state and the interpolating operator cancel. Thus, for sufficiently large separation between the source, sink, and operator insertion times, the ratio asymptotically approaches a value proportional to the matrix element h(p , s )|T g µν |h(p, s) with exponentially suppressed excited state contamination. Specifically, for the four states of interest h ∈ {π, N, ρ, ∆}, the computed ratios are related to the matrix elements of interest as are Euclideanized using the Euclidean-to-Minkowski matching relation where i δµt generates a factor of i on the temporal component. It follows directly from F µν ∝ [D µ , D ν ] that the Euclidean and Minkowski matrix elements of the gluon EMT are related as Each ratio is associated with a different set of momenta ∆ µ and P µ , operator basis element R, i, and spin channel s → s , all of which define a set of kinematic coefficients K h,j ss ;R,i (P µ , ∆ µ ) for the bare GFFs for irrep R in the decomposition The GFFs are real, but the kinematic coefficients and ratio measurements are generically complex, so the real and imaginary parts of each ratio measurement provide independent constraints on the GFFs; we thus treat each part as a separate real-valued ratio associated with real coefficients. We discard any ratio for which all kinematic coefficients are zero.     16. For each of the states π, N , ρ, and ∆, the figure shows how different combinations of the discretized momenta p, p , and ∆ = p − p used in this calculation are associated with discrete t bins as described in the text. Each marker represents a collection of different momentum combinations that result in the same value of −t. Colors correspond to different values of |∆| 2 . The area of each marker is proportional to the number of associated momentum combinations. Gray bands indicate each t bin, and markers are associated with the band that contains their central point. Table XI lists the number of bins for each state.

using the dispersion relation
Although the kinematic coefficients and values of t = ∆ 2 associated with each ratio are functions of the hadron mass m h and lattice spacing a and so in principle are only known up to some uncertainty (correlated with the ratios), these errors are subdominant, so we neglect them and evaluate the coefficients using a = 0.1167 fm and the numerical values of m h listed in Table XI, obtained from single-state fits to the effective mass as shown in Fig. 15.
We associate each ratio with a "t bin" so that we can estimate model-independent values of the GFFs at discrete values of t. Bins are defined by grouping together any two ratios associated with values of t that differ by less than 0.03 GeV 2 , with no additional restriction on the maximum width of each bin. 8 We define the value of t for each bin as the average over t for all ratios in the bin. Figure 16 illustrates the resulting associations for each hadron. There is a one-to-one correspondence between t-bins and values of |∆| 2 in the case of the baryons, but not in that of the mesons, which is due to the smaller masses of ρ and π compared to N and ∆. Within each t-bin, we average any ratios  -bins), and the number of ratios before and after combining ratios with kinematic coefficients related by an overall sign.
associated with kinematic coefficients related by an overall sign within each bootstrap draw, with each ratio multiplied by the appropriate sign. We do not combine ratios from different irreps, as they are renormalized differently, and we continue to keep SS and SP ratios separate. This additional averaging helps to compensate for gauge noise, providing clearer signals for subsequent fitting. The resulting averaged ratios R Rtc (t f , τ ) in momentum bin t for irrep R are no longer associated with specific momenta, irrep basis elements, spin channels, or real/imaginary parts, and are instead associated simply with some particular set of kinematic coefficients indexed by c. The resulting reduction in data volume is significant, as tabulated in Table XI.
To extract the t f τ 0 asymptotic values of the ratios R Rtc (t f , τ ), which we denote by R Rtc with no argument, we perform correlated χ 2 fits of a constant to each ratio for every triangular connected region in the (t f , τ ) plane that satisfies t f < 25, τ > 4, and t f − τ > 4. The minimum cuts on τ and t f − τ guarantee a transfer matrix exists between the source and operator insertion, and the insertion and sink. t f ≈ 8 is the approximate time after which the effective masses are consistent with a single state for all hadrons, momenta, and smearings, and the upper bound t f ≈ 25 removes the bulk of the noise-dominated region. To combine the separate SS and SP ratios, we simultaneously fit the same region in each to a single value of R Rtc . We combine the results of fits to different regions of (t f , τ ) using a scheme inspired by Bayesian model averaging [136]. Denoting by r m the values of R found by fits to each region m, we associate each fit with a weight [137] where for the fit to region m, δr stat m is the statistical error found by the fit and p m = Prob( is the p-value of the fit. Normalizing the weights such that m w m = 1, we obtain the mean value asR = m w m r m and the total variance (δR) 2 as the sum of statistical and systematic contributions defined as [136] (δR stat ) 2 = We find that typically δR stat ≈ δR syst . In practice, we perform "central-value fits" to the median of R(t f , τ ) over bootstraps, from which we compute a set of weights w * m and averaged error δR * . For subsequent error propagation, we compute bootstrapped fit results by averaging over fits within bootstraps b using the central-value weights w * m , then rescaling to obtain a set of results whose spread reproduces δR * . In detail: we fit all R b (t f , τ ) to obtain r bm for only the subset of highest-weight regions making up 99% of the total weight, which reduces the computational cost by excluding the bulk of fit regions. We then average to obtainR b = m w * m r bm , with w * m suitably re-normalized to account for the exclusion of low-weight fits. The spread inR b obtained in this way only reproduces δR * stat , so we rescale each set ofR b around their mean by δR * /δR * stat . Note that we use the same covariance matrix for both the central-value fits and bootstrap fits, computed over R b (t f , τ ) using an outlier-robust ±1σ percentile definition of the error. 9 As shown in Figs. 17 and 18, the ratios typically exhibit plateaus in t f , suggesting that excited-state contamination will not significantly affect the results. To check this, we perform a simplified version of this analysis for the pion, nucleon, and rho using a two-state ansatz and only bootstraps from the highest-weight fits; the resulting GFFs are consistent within uncertainties in all cases. The precision of the ratio data for the delta baryon does not admit two-state fits.
The results of this fitting and averaging procedure are generically robust against varying the lower bounds on fit regions, but increasing the upper bound on t f results in sudden catastrophic increases in error and destabilization of  central values. This effect can be traced back to fits to pure-noise regions which are excluded by the t f cut. These fits are apparently good, as measured by their χ 2 /d.o.f. or p-values, but the loss of Gaussianity in noise regions (the onset of which occurs at t f ≈ 25 in the nucleon two-point correlator as diagnosed using both cumulant expansions [140] and Shapiro-Wilk testing [141]) renders these metrics of fit quality meaningless. Noisy regions must thus be excluded using a t f cut to prevent them from dominating the averages. We chose to use the ad hoc weight definition described above because we found it to be practically more robust against this effect (due to the inverse variance factor) than the better-motivated AIC weighting of Ref. [136].
In the analysis described above, the choice to rescale the bootstraps around their means amounts to an assumption that systematic errors due to the choice of fit range have the same correlation structure as the statistical errors. This is different from the typical assumption of uncorrelated systematics [60], but both are strong assumptions. To check that this choice does not bias our results, we applied the subsequent analysis to the nucleon data with all correlations between ratios either artificially scaled down by overall factors or completely neglected, as well as using best fits or fits to a fiducial (t f , τ ) region rather than averaging, and found no systematic shift in the results. Further work is needed to more gracefully reconcile frequentist resampling techniques with Bayesian model averaging methods and avoid the need for such ad hoc constructions. For further analysis, we take the median over (rescaled) bootstraps for the central value of each R Rtc and construct their covariance matrix using the outlier-robust estimator noted above. Parametrizing the fit results as central values and a covariance matrix amounts to modeling their distribution as a multivariate Gaussian. We check this assumption by examining the bootstrap distribution of fit results, and find that histograms of marginal distributions are either consistent with or contained in their Gaussian approximations. We have also checked that bootstrapping through the further analysis detailed below produces marginal distributions consistent with or narrower than the ones presented in the main text, which are obtained using linear error propagation from this Gaussian model.

Constraint fitting
To compactify notation, throughout this section we use 1 for τ whenever an irrep label appears in a subscript, and switch to vector notation for the kinematic coefficients and GFFs, i.e. K j , G j ⇔ K, G.
The procedure described in the previous section yields a set of measurements which constrain the bare GFFs of each irrep R ∈ {τ 3 } separately as where K and G are N h -element vectors over the set of different GFFs, t indexes the discrete t-bin, and c indexes the different combined ratios with shared kinematic factors as described in Sec. A 4. Extracting the renormalized GFFs from these constraints, as well as subsequent model fitting of the GFFs, requires careful treatment to avoid the d'Agostini bias [79]. This bias is an effect caused by violation of implicit Gaussianity assumptions in correlated χ 2 fitting by non-Gaussianity arising from multiplication by the renormalization factors. To circumvent it, we use a Bayesian version of the "penalty trick" [79], performing combined fits of data from both irreps to estimate the bare GFFs G 1t and update the renormalization factors Z 1 , Z 2 → Z 1 , Z 2 . The updated renormalization may be applied immediately to obtain the renormalized GFFs G t = Z 1 G 1t , or deferred until after subsequent model fitting to again circumvent the bias as discussed below. We defer detailed discussion of the bias and the derivation of the fitting procedure presented here to Sec. A 6. Our procedure estimates a Gaussian approximation of the posterior distribution where G 1t are the t-bin-dependent bare GFFs for irrep τ 1 , Z 1 and Z 2 are the updated renormalization factors which are shared across all t bins, R 1 and R 2 represent the full set of ratio fit results R Rtc , the factor p(R 1 , R 2 ) is the usual uninteresting data normalization factor in Bayes's theorem, the prior p (Z 1 , Z 2 ) is the multivariate Gaussian defined by Eq. (A6), and the likelihood L is multivariate Gaussian, defined in terms of the measured meansR Rtc and covariance matrix Σ −1 Rtc,R t c of the ratio fit results R Rtc . Equation (A34) should be read as one overall distribution for all t bins and not a set of separate equations for each bin. The data only constrain the ratio of the Z factors and not their overall magnitude, which corresponds to a flat direction in the likelihood function that is only regulated in the posterior by p(Z 1 , Z 2 ). Note that we have left implicit the uniform prior over G 1t to emphasize that, although our analysis is phrased in Bayesian language, it involves no informative priors.
We estimate the parameters of the posterior distribution using two stages of fitting. In the first stage, we introduce a separate ratio (Z 1 /Z 2 ) t for each t bin, defining an extended version of the likelihood which we approximate as a Gaussian distribution around the maximum likelihood parameters G * 1t and (Z 1 /Z 2 ) * t . We obtain these parameters by fitting each t bin separately, using linear error propagation to obtain covariances between the parameters (both within and between t bins). The posterior of interest can then be written in terms of this extended likelihood function as which, after evaluating the δ functions, provides a new merit function which we can re-fit (i.e. minimize and expand about) to estimate the parameters of the Gaussian posterior. This second stage of fitting incorporates the measured distribution of Z factors [Eq. (A6)] and the constraint that the ratio Z 1 /Z 2 is the same for all t bins. We again estimate the covariances of this distribution using linear error propagation. For the pion and nucleon, χ 2 /d.o.f. ≈ 1 and p > 0.1 for all first-stage fits to individual t bins; for most t bins, p ≈ 1. The second-stage fits are of similarly high quality. However, for the ρ and ∆, we observe that p 1 in fits to t bins with more than ≈ 600 constraints. We trace the source of this effect to finite-statistics limitations, which we circumvent by combining constraints. When more than 600 constraints are present in a t bin, we apply a "pair binning" procedure to that bin to reduce the number of constraints before fitting. To choose which pairs of constraints are binned together in a way that heuristically minimizes loss of orthogonality in the set of constraints, we use a greedy algorithm which repeatedly associates the two unpaired constraints K and K with the least angle cos −1 (K · K /|K||K |) between them until all constraints are paired (with possibly one left unpaired, which is retained). Paired constraints are combined by taking weighted averages at the per-bootstrap level, using weights proportional to the number of ratio measurements averaged into each constraint. For the ρ, no t bin requires more than one application of this procedure, while for the ∆, some bins require two applications. After applying this procedure, first-stage fits to the pair-binned constraints for the ρ and ∆ satisfy p > 0.1 for all bins, with p ≈ 1 for most; second-stage fits are also of high quality. This procedure could have instead been applied to combine the (t f , τ )-dependent ratios before fitting, and less naive clustering algorithms than the one used here may allow more effective use of the data; we did not explore either direction in this study, but they are interesting topics for future work. To check that pair binning does not bias the results, we instead discard random subsets of the data to equivalently reduce the number of constraints and find consistent but noisier results. The statistical limitations addressed by pair binning may be artificial and due to the limited number (B = 1000) of bootstraps, as we observe similar failures in the fits for the pion and nucleon when using B = 200 that are resolved when using more; 10 however, we find that our results do not depend significantly on the number of bootstraps B, after pair binning or discarding constraints to ensure all first-stage fits are of good quality.
The renormalized GFFs are distributed as the product under the posterior distribution of the bare GFFs G 1t and renormalization factor Z 1 for irrep τ 1 . Starting from the Gaussian approximation of the posterior computed using the procedure described above, we obtain the uncertainties of the renormalized GFFs presented in the main text using linear error propagation; we find this approximation to be consistent with the spreads in Z 1 G 1t computed over samples drawn from the Gaussian posterior. However, as discussed in Appendix A 6, the renormalized GFFs are sufficiently non-Gaussian that subsequent fits of the models of Eqs. (2) and (3) to them would again be victim to the d'Agostini bias. We instead fit the models to the bare GFFs G 1t , then renormalize afterwards by multiplying with Z 1 . Due to the structure of the model functions, the factor of Z 1 may be absorbed into the parameters α and α k , defining the renormalized fit parameters presented throughout this work.
Reference [60] instead circumvented the d'Agostini bias by neglecting additional correlations between renormalized constraints induced by common factors of Z. The results obtained using the method presented here are consistent with the ones from that study, but with narrower and more correlated uncertainties on the GFF estimates and wider ones on the model fits and densities. For all GFFs, the results of this sampling procedure are consistent within error with the results of the procedure used here. Employing the sampling procedure while accounting for correlations induced by shared Z factors would require performing the entire analysis for each sample from the Z distribution, including the expensive density estimations, which would require significant additional computational effort.
As mentioned throughout the discussion above, starting from the model of the ratio distribution as Gaussian, we use linear error propagation to propagate uncertainty through the rest of the analysis and obtain the presented results, amounting to repeatedly approximating intermediate distributions as Gaussian. Other than the checks of these approximations described above, we have also checked that bootstrapping through the entire analysis, as well as just through the first stage of fitting and using the bootstrap results to construct a covariance matrix before the second stage, produces marginal distributions of GFFs consistent with or contained within the marginal distributions obtained with linear error propagation.

Gaussianity and the d'Agostini bias
In this section we discuss the d'Agostini bias, identify where the non-Gaussianities that trigger it arise in our analysis, introduce and discuss the penalty trick fitting procedure in a Bayesian framework, and motivate and derive the modified version described in Sec. A 5.
In its simplest form, the d'Agostini bias occurs when performing a correlated χ 2 fit of some linear model [or a nonlinear model whose form accommodates arbitrary rescaling, like the model ansatzë Eqs. (2) and (3)] to some data which has been multiplied by an overall normalization factor with a large relative uncertainty; the result is different than what is obtained by first fitting then normalizing after, and thus obviously incorrect. This occurs because the χ 2 fitting procedure takes the covariance matrix of the data as input, and thus implicitly truncates the data distribution to Gaussian; products of Gaussian-distributed variables are not Gaussian distributed, and the bias occurs when this truncation yields a poor approximation of the true product distribution. While resampling through a fit allows for treatment of non-Gaussianity in distributions of fit parameters due to nonlinear model functions, it cannot correct for the d'Agostini bias, which occurs because the fit assumes an inaccurate representation of the data.
Given our multivariate Gaussian models of the distributions of the bare ratios and renormalization factors, the bare GFFs are Gaussian but the renormalized ratios and GFFs are not. The bare ratios are Gaussian by assumption and constrain the bare GFFs linearly per Eq. (A33), so the bare GFFs inherit the Gaussianity of the ratios. However, the renormalized ratios Z R R Rtc are non-Gaussian, as shown in Fig. 19a and discussed in the caption. It follows that the renormalized GFFs, which are linearly constrained by the renormalized ratios, are also non-Gaussian, intrinsically and independently of how we extract them, as shown in Figs. 19b and 19c. These non-Gaussianities trigger the d'Agostini bias both when fitting ratios to extract GFFs, as demonstrated in Fig 19a, as well as subsequently when fitting the GFFs to model functions, as shown in Figs. 19b and 19c. The fitting procedure described in Sec. A 5 circumvents the bias in the former case using the penalty trick, and in the latter case by extracting the Gaussian-distributed bare GFFs for one irrep and allowing the problematic multiplication by a Z factor to be deferred until after fitting models to the bare GFFs. Note that while the histograms of marginal distributions shown in Fig. 19 naively appear close enough to Gaussian to justify approximation as Gaussian, inspection of the joint histograms reveals the asymmetry of the distribution that leads to the bias.
The penalty trick is a common prescription for circumventing the d'Agostini bias [79]. Our choice to phrase the fitting problem as an estimation of a posterior distribution (as described in Sec. A 5), with the measured distribution of the renormalization factors entering as a prior to be updated, amounts to a Bayesian reframing of this technique. Generally, for a fit of a model function f (θ) to some data y times a normalization factor Z, where y ∼ N (ŷ, Σ y ) and Z ∼ N (Ẑ, σ 2 Z ) are Gaussian, but a Gaussian N ( " yZ, Σ yZ ) is a poor approximation of the distribution of the product yZ, the penalty trick prescribes the replacement allowing a fit using the original covariance matrix Σ y , assumed to be a good description of the data. This comes at the cost of replacing the fixed normalization Z with an additional nuisance parameter z which is constrained to be consistent with the provided Z and discarded after fitting. In the limit σ Z → 0 the two fit procedures are equivalent. While usually motivated as an ad hoc frequentist procedure, the penalty trick can be more naturally understood in a Bayesian context, wherein it is structurally equivalent to updating a prior for Z with the data then marginalizing over it, assuming a Gaussian posterior. The right-hand side of Eq. (A37) can be interpreted as a log-likelihood and log-prior for the data and z, defining a posterior distribution via Bayes's theorem as p(θ, z|y) = L(y|θ, z) p(z) p(θ) / p(y) where p(y) is the data normalization and p(θ) is a trivial factor of the uniform distribution added as a prior for the fit parameters. Fitting the penalty trick χ 2 to obtain the best-fit θ * and z * and fit parameter covariance matrix Σ * θ,z corresponds to approximating the posterior distribution as Gaussian, i.e.
suppressing normalization factors. Discarding z after fitting corresponds to marginalizing over z in the posterior, as marginalizing over a dimension of a multivariate Gaussian is equivalent to dropping it. The generalization to the case of multiple different normalization factors for different subsets of the data is straightforward: the prior p(z) becomes multidimensional, and now  . 19. Examples of non-Gaussianities in various distributions in the analysis of the nucleon data (left panels), and the resulting effects of d'Agostini bias (right panels). Note: these plots are an illustration of the bias and are not the final results of our calculation. In left panels, the orange features show the Gaussian approximations to these distributions obtained with linear error propagation. In the joint histograms, the ellipses denote the 3σ contour and the dot denotes the mean. In right panels, blue bands show results obtained by fitting a bare quantity first then renormalizing afterwards to circumvent the bias, whereas orange bands are biased fits to data renormalized before fitting. (a) For irrep τ only to obtain A1 then renormalizing afterwards, example of a joint distribution of the renormalized GFF A(t) = Z1A1(t) in two different t bins, and fits of a tripole model to these GFFs. (c) For renormalized GFFs obtained using the fitting procedure described in Sec. A 5 incorporating data from both irreps, joint distribution of the same renormalized GFFs as in (c), and fits of a tripole model to these GFFs.
where i indexes different subsets of the data.
We modify the penalty trick procedure to estimate the bare GFFs G 1t (corresponding to θ/z 1 ) instead of the non-Gaussian renormalized GFFs G t (corresponding to θ). The modification singles out one particular normalization as special, multiplying it onto the model function f so that the data are modeled as If the model function f is linear in the parameters (e.g. K · G is linear in the GFFs G), then this procedure extracts θ = θ/z 1 (corresponding to G 1t ) rather than θ (corresponding to G t ). In this modified form one still (trivially) marginalizes over all z i for i = 1, but z 1 must be retained to examine θ = z 1 θ (corresponding to renormalizing the bare GFFs as G R = Z 1 G 1 ). While fitting procedures exist for treating the d'Agostini bias other than the penalty trick [143], a model function and a data distribution define a distribution of model parameters (e.g. GFFs) independent of the choice of biascircumventing fitting procedure. The Bayesian framework makes clear that the renormalized GFFs extracted by this procedure may themselves be non-Gaussian, such that subsequent fits are also vulnerable to the bias. This will hold independent of the fitting procedure used.

Appendix B: DENSITY DEFINITIONS
This section lists the expressions for the energy, pressure, and shear force distributions in the 3D Breit frame (BF3), 2D Breit frame (BF2), and infinite momentum frame (IMF) used to generate the results of Sec. III. To simplify the expressions below, we define bracket notation for the relevant integrals, where I is a generic integrand and J 0 is a Bessel function of the first kind. We compute the presented densities, defined by Eq. (B1) and the expressions below, using numerical integration. The analysis of Ref. [60] propagated uncertainty on model parameters into the densities by sampling from the multivariate Gaussian distribution of the model parameters, evaluating the integrals for each draw. The large number of densities considered in the present study make this approach impractical. We instead used linearized error propagation: by differentiating under the integral sign with respect to model parameters θ, we obtain the Jacobian J(r; θ) i = ∂I(r; θ)/∂θ i , where I(r; θ) is an integral evaluated to obtain a density at radius r, as a matrix of integrals that can each be evaluated numerically. The covariance matrix for the r-dependent density is then obtained as Cov[ρ(r), ρ(r )] = ij J(r; θ) i Cov[θ i , θ j ]J(r ; θ) j where Cov[θ i , θ j ] is the covariance matrix of the parameters of the model integrated to obtain the density.
The model functions are linear in some parameters (α for the multipole and α k for the modified z-expansion) but not others (multipole masses), so this approach is approximate. However, for all densities for the nucleon and pion, as well as for the monopole densities for the ρ meson, we found consistent results for all r by computing integrals for samples from the distribution of renormalized model parameters. For tripole models of the nucleon GFFs in the 3D Breit frame, we also checked our numerically integrated density results against ones derived from the closed-form solution using linear error propagation from the tripole model parameters α and Λ, and found indistinguishable results. The uncertainties on the densities presented in the main text are derived from the distribution of renormalized model parameters, obtained by combining the uncertainties of the bare α and α k parameters and fitted values of Z 1 using linear error propagation as described in Sec. A 5. We found consistent results by computing densities from bare model parameters, using linear error propagation to obtain correlations between Z 1 and the resulting bare densities, then applying the renormalization factor Z 1 and propagating uncertainities either linearly or by drawing correlated samples of Z 1 and the bare densities and multiplying within samples.
The mass mean square radii are defined identically for all hadrons as .

(B6)
The mechanical radius results presented throughout this work are computed by numerically approximating the integrals with the trapezoidal rule, evaluated at 500 values of the integrands evenly spaced in 0 ≤ r ≤ 2 fm. We obtain error estimates using linear error propagation from the values of p(r) and s(r) at each r, computed as described above, and corresponding to the results presented in = 12/Λ 2 from Eq. (B2) and find it yields results consistent within uncertainty. The shear and pressure densities for other models, frames, and hadrons are comparably small by r = 2 fm, so we expect this quality of approximation to hold generally.
The subsection below lists the various densities computed for each hadron. In all expressions for the densities and radii, we use the definitions ∂ 2 = 1 The BF3 densities of the ∆ baryon were derived in Ref. [37] and are