Gluon gravitational form factors of the nucleon and the pion from lattice QCD

A future Electron-Ion Collider will enable the gluon contributions to the gravitational form factors of the proton to be constrained experimentally for the first time. Here, the first calculation of these form factors from lattice Quantum Chromodynamics is presented. The calculations use a larger-than-physical value of the light quark mass corresponding to $m_\pi \sim 450$ MeV. All three form factors, which encode the momentum-dependence of the lowest moment of the spin independent gluon generalised parton distributions and are related to different components of the energy-momentum tensor, are resolved. In particular, the gluon $D$-term form factor, related to the pressure distribution inside the nucleon, is determined for the first time. The gluon contributions to the two gravitational form factors of the pion are also determined, and are compared to existing lattice determinations of the quark contributions to the gravitational form factors and to phenomenology.


I. INTRODUCTION
A defining challenge for hadronic physics research is to achieve a quantitative understanding of the structure of the proton and other hadrons in terms of their fundamental quark and gluon constituents. Generalised parton distributions (GPDs) [1][2][3] provide a framework for a three-dimensional encoding of this structure. These distributions (see Refs. [4][5][6] for reviews) combine and generalise the features of elastic form factors, which describe the charge and magnetisation distributions of the hadron as seen by a photon of a given virtuality, and parton densities, which describe the longitudinal partonic composition of a fast moving hadron. Moreover, the nucleon GPDs encode, for example, the nucleon mass and spin, the quark and gluon contributions to the nucleon angular momentum, and various inter-and multi-parton correlations. They are also directly related to the 'mechanical properties' of the nucleon system, such as the pressure distributions and shear forces [7][8][9]. Given their importance in hadron structure, there are significant experimental and theory efforts targeted at determining GPDs, especially for the proton.
In particular, proton and neutron quark GPDs have been constrained in limited kinematic regions by deeply-virtual Compton scattering (DVCS) and deeply-virtual meson production (DVMP) experiments at Thomas Jefferson National Accelerator Facility (JLab) [10], HERA [11], and COMPASS [12] (Refs. [13][14][15] summarise the world datasets). Ongoing studies at COMPASS and within the 12 GeV program at JLab will significantly improve extractions of these quantities. The quark GPDs have also been studied theoretically in a number of frameworks, including determinations of their connection to different experiments and evolution with factorisation and renormalisation scales [1][2][3][4][5][6], estimates of their forms in hadronic models (see e.g. Refs. [15,16] for reviews), and calculations of their lowest few Mellin moments for the pion [17] and nucleon [18][19][20][21][22][23][24][25] in lattice QCD (LQCD) (see e.g. Ref. [26] for a review). Gluon GPDs, on the other hand, are as-yet unknown from experiment or theory [27]. Performing first measurements of these quantities is a key goal of the planned Electron-Ion Collider (EIC) [28,29], and theory constraints on the gluon GPDs will provide important information as the physics case for an EIC is refined.
This manuscript presents the first LQCD determination of the complete set of gluon gravitational form factors (GFFs) of the nucleon, which are defined as the lowest moments of the spin-independent gluon GPDs. The calculations are undertaken with a larger-than-physical value of the light quark mass that corresponds to a pion mass m π ∼ 450 MeV. All three gluon GFFs of the nucleon are determined at discrete values of the squared momentum transfer t up to |t| ∼ 2 GeV 2 , as are the two gluon GFFs of the pion. The results are presented in a modified minimal subtraction (MS) scheme at a renormalisation scale µ = 2 GeV, by performing a non-perturbative renormalisation using a RI-MOM scheme [30] and a perturbative matching to MS. Mixing of the gluon GFFs with the corresponding quark distributions is neglected; lattice perturbation theory calculations [31] indicate that such effects are likely to be small compared with the statistical and other systematic uncertainties of this study. For both the pion and the nucleon, the gluon momentum fraction is found to be approximately 0.5-0.6, somewhat larger than the phenomenological value in both cases, while the fractional contributions of gluons to the nucleon lightcone momentum and angular momentum are found to be consistent within uncertainties. The |t|-dependences of two of the three nucleon gluon GFFs are consistent with dipole forms, while the third GFF shows no |t|-dependence and is consistent with zero. Comparing the results with those of previous calculations of the quark GFFs using similar lattice discretisations and at similar values of the quark masses reveals that the gluon radii defined as the slopes of the nucleon gluon GFFs at t = 0 are larger than the corresponding quark radii for each form factor. The pion gluon GFFs also have dipole-like dependences on |t|, but are consistent with the corresponding quark GFFs within uncertainties, revealing no clear ordering of pion quark and gluon radii. Compared with the nucleon gluon GFFs, the pion gluon GFFs define consistently smaller radii, consistent with the ordering of the nucleon and pion charge radii (defined from the electric form factors) determined experimentally.
In the following section, the quark and gluon GPDs and GFFs of the nucleon and pion are defined. Section III details the LQCD calculations that are performed to extract the gluon GFFs, while the results of those calculations are presented in Section IV. The extracted gluon GFFs are compared with the corresponding quark GFFs, which have been previously calculated using similar lattice discretisations at quark masses corresponding to a similar value of the pion mass. Earlier LQCD results for the gluon momentum fractions of the nucleon and pion, defined as the forward limits of the appropriate GFFs, are also collated for comparison. Finally, Sec. V highlights the conclusions that can be drawn from this study.

A. Nucleon
GPDs encode the three-dimensional distribution of quarks and gluons in the nucleon [1][2][3][4][5][6]. In the deep inelastic regime, the leading contributions arise from the lowest-twist operators. For the nucleon, the leading spin-independent quark and gluon distributions are twist-two [32][33][34], and, following the conventions of Ref. [5], can be expressed in terms of matrix elements of non-local light-ray operators as: where ψ q is a quark field of flavour q, G a µν = (∂ µ A a ν − ∂ ν A a µ + gf abc A b µ A c ν ) is the gluon field strength tensor built from the gluon field A a µ , and the ellipses denote structures with twist greater than two. Here, n µ is a light-like vector with n 2 = 0, the momenta and spins of the initial and final nucleons are (p, s) and (p , s ) respectively, and it is convenient to define P = 1 2 (p +p), ∆ = p −p, t = ∆ 2 , Bjorken x = 1 2 ∆ 2 /p·∆, and skewness ξ = − 1 2 n·∆/n·P . The path-ordered gauge links in the fundamental and adjoint representations are where t a are SU(3) generators in the fundamental representation and f abc are the structure constants defining the adjoint representation. The inclusion of the gauge links in Eqs. (1) and (2) ensures the gauge-invariance of these expressions (in the case of the gluon operator, alternate gauge link choices are also possible [35]). Braces denote symmetrisation and trace-subtraction in the free indices, i.e., a {µ b ν} = 1 2 (a µ b ν + a ν b µ ) − 1 4 g µν a α b α , and the covariant normalisation of states p , s | p, s = 2p 0 (2π) 3 δ s s δ (3) (p − p) is used along with the spinor normalisation U (p, s) U (p, s) = 2M N . In the forward limit, the distributions H a (x, 0, 0), for a = {u, d, . . . , g}, define the familiar unpolarised quark and gluon PDFs, i.e., H q (x, 0, 0) = q(x) and H g (x, 0, 0) = xg(x).
The operator product expansion (OPE) relates the Bjorken-x (Mellin) moments of the GPDs H a (x, ξ, t) and E a (x, ξ, t) to matrix elements of local twist-two operators. The focus of this work is on the lowest moments of the spin-independent gluon GPDs, which are related to the nucleon matrix element of the gluon contribution to the (traceless, symmetric) energy-momentum tensor (EMT) 1 , and are encoded in three scalar GFFs that are functions of t [9]: where (again, following the conventions of Ref. [5]) An exactly analogous decomposition exists for matrix elements of the quark contribution of flavour q to the traceless part of the EMT: For each q = {u, d, . . . }, the GFFs are related to the lowest Mellin moments of the relevant unpolarised GPDs defined in Eq. (1): and similarly the gluon GFFs are related to the GPDs defined in Eq. (2): Since the quark and gluon pieces of the EMT are not separately conserved, the individual form factors A a (t), B a (t) and D a (t) are scale-and scheme-dependent, although the total form factors A(t), B(t), D(t), where X(t) ≡ a X a (t) with a = {u, d, . . . , g}, are renormalisation-scale invariant. The GFFs A a (t) encode the distribution of the nucleon's momentum among its constituents (and momentum conservation implies A(0) = 1), while the angular momentum distributions are described by J a (t) = 1 2 (A a (t) + B a (t)) (and total spin constrains J(0) = 1 2 ). The D a (t) terms encode the shear forces acting on the quarks and gluons in the nucleon while their sum D(t) determines the pressure distribution [7][8][9].

B. Pion
The spin-independent pion GPDs are defined by pion matrix elements of the lowest-twist light-ray quark and gluon operators: for q = {u, d, . . .}, and where the notation is as in Eqs. (1) and (2). A covariant normalisation of pion states has been used: . The lowest moments of these GPDs are related to the pion matrix elements of the quark and gluon pieces of the traceless EMT, which are described by two scalar GFFs for each flavour a, labelled A and similarly for the quark operators, Just as for the nucleon, the GFFs which describe pion matrix elements of the EMT correspond to the quark and gluon gravitational form factors of the pion, and can be expressed as Mellin moments of the pion GPDs: The forward limit A

III. LATTICE QCD CALCULATION
In this work, a single ensemble of isotropic gauge-field configurations is used to determine the matrix elements corresponding to the gravitational form factors of the nucleon and pion, Eqs. (4) and (11), respectively. Simulations are performed with N f = 2 + 1 flavours of quarks, with quark masses chosen such that m π ∼ 450(5) MeV. A cloverimproved quark action [37] and Lüscher-Weisz gauge action [38] are used, with the clover coefficient set equal to its tree-level tadpole-improved value. The configurations have dimensions L 3 × T = 32 3 × 96, with lattice spacing a = 0.1167(16) fm [39]. Details of this ensemble are given in Table I and in Ref. [40]. Subsections III A, III B and III C define the Euclidean-space gluon operators studied here, detail the renormalisation prescription, and outline the extraction of the gluon GFFs from Euclidean correlation functions, respectively.

A. Operators
To determine the spin-independent gluon GFFs, matrix elements of the gluon operators 2 are constructed, where the brackets denote symmetrisation and tracelessness in the µ and ν indices by In Euclidean space, the unrenormalised gluon operators are defined using the clover definition of the discretised Euclidean-space field-strength tensor derived from the combination of plaquettes which are in turn built from gauge link fields that have been subject to Wilson flow to flow-time t = 1.0 [41] in order to increase the signal-to-noise ratio of the calculation. In a previous study of these operators in a φ meson [42,43], the effects of different flow times and different choices of smearing prescription on the bare matrix elements have been found to be mild. Since a non-perturbative renormalisation procedure is used here (discussed in the next section), the differences between bare matrix elements calculated with different smearing prescriptions will be compensated for by differences in the renormalisation. Because of the reduced symmetry of the lattice geometry, the discretised operators transform in particular representations of the hypercubic group H(4). Specifically, the operators in Eq. (14) subduce into traceless, symmetric representations of H(4), two of which do not mix with same or lower-dimension operators (labelled τ and τ (6) 3 in the notation of Refs. [44,45]). In Minkowski space 3 , a basis of operators in the three-dimensional τ 2 Since the Lüscher-Weisz gauge action is used in this work, the Belinfante procedure [36] produces a gluon EMT that has an additional contribution that is higher-order in a. This term is neglected in the present work. 3 The Euclidean operators are related to these by G while a basis for the six-dimensional τ representation is: All basis operators in each of these two representations are studied here. Within each representation, the renormalisation of the different operators are related by symmetries, while the renormalisations of operators in the two different representations are only constrained to be the same in the continuum limit. Studying both representations thus permits a test of the discretisation artefacts in this calculation.

B. Renormalisation
The unrenormalised operators in Eq. (14) mix with the flavour-singlet quark operators Q µν = q∈{u,d,s} ψ q γ {µ i ↔ D ν} ψ q such that the renormalised gluon operator O ren. µν (in any particular scheme) is described by O ren. µν = Z gg O µν + Z gq Q µν . It was shown in Ref. [31] that the mixing of the quark operator into the gluon operator is a few-percent effect, using a one-loop perturbative renormalisation procedure and a similar action to the one used here. Consequently, this mixing is assumed to be negligible relative to the statistical uncertainties of this calculation and is neglected here.
The bare lattice operators described in the previous section are renormalised via a non-perturbative RI-MOM prescription [30,46], similar to that recently investigated for gluon operators in Refs. [47,48]. A perturbative matching is used to relate the renormalised operators to the MS scheme. A bare lattice operator O latt is thus renormalised as 4 : The conversion factor R MS (µ 2 , µ 2 R ) from the RI-MOM scheme to MS is calculated in continuum perturbation theory [49], while the RI-MOM renormalisation constant Z RI-MOM which relates the bare and tree-level amputated Green's functions Λ bare/tree O (p) for the operator O in a Landau-gaugefixed gluon state of momentum p 2 = µ 2 R . Here, Z 1/2 g (p 2 ) denotes the gluon field renormalisation. For the particular operator of interest, O µν , the tree-level amputated Green's function can be expressed as [50]: As discussed in Ref. [50], and also noted in Ref. [47], only the first structure in this expression is protected from mixing with the gauge-variant parts of the energy-momentum tensor. Consequently, choosing renormalisation conditions that only involve this term allows a purely multiplicative renormalisation procedure even for gauge-fixed states. The operators in Eqs. (17) and (18) can be arranged into the forms 5 1 , 4 In the general case this is a matrix equation that accounts for mixing among a set of bare operators O latt i . 5 For the operators in representation R = τ , and the combination (1/ L/a T /a β am l ams N cfg 12 24 6.1 -0.2800 -0.2450 24600 TABLE II: Details of the gauge ensemble used to determine the non-perturbative operator renormalisation. Other than the lattice volume L 3 × T , which is smaller, the parameters are the same as those of the ensemble (detailed in Table I) used for the primary calculation. A total of N cfg configurations are used, generated in 50 streams of configurations, with samples separated by 10 hybrid Monte-Carlo trajectories in each stream.
with no summation over repeated indices implied. For these operators, the constraint σ = τ = α = β is sufficient to isolate the desired term in Eq. (21), leading to In general, the construction of the forward, amputated bare Green's function Λ bare O (p) will depend on the operator O. For the operatorsÔ αβ considered here, with the same conditions on the external states as applied in Eq. (23), the additional condition p τ = 0 (where τ is the Lorentz index of the external gluon fields) permits a simple form: where no summation over repeated indices is implied, and where the second line follows by substitution of the trace of the gluon propagator for one of the gluon terms in the denominator. The RI-MOM renormalisation constant Z RI-MOM O (p 2 ) can thus be determined by combining Eqs. (23) and (24) as prescribed by Eq. (20), and taking N c = 3: The gluon three-point and two-point functions that appear in this expression are computed on the ensemble detailed in Table II, which has the same bare parameters, but smaller lattice volume and an order of magnitude more configurations, as the ensemble used for the main calculation. The gluon fields are computed from Landau-gauge-fixed links (gauge-fixed using an iterative procedure with tolerance 10 −5 ) U µ (x): which holds up to O(a 2 ) corrections. Momentum-space lattice gluon fields are defined by the discrete Fourier transform: where L µ denotes the number of lattices sites in dimension µ and where the discretised momenta accessible on the finite lattice volume arep Gluon two-point functions D τ τ (p) = Tr[A τ (p)A τ (−p)] are constructed for all four-momentap µ corresponding to p µ with µ n 2 µ ≤ 36 (Eq. (28) 17) and (18), calculated on the ensemble detailed in Table II, with cuts on four-momenta such that (µ 2 R ) using flowed or unflowed fields in the propagators will agree up to discretisation artefacts, and comparing the two determinations provides a measure of such effects. Gluon three-point functions Ô αβ Tr[A τ (p)A τ (−p)] are constructed on each configuration by correlating the gluon two-point functions with the operatorsÔ αβ , computed as described in Section III A and projected to zero four-momentum, and subtracting the vacuum contribution.
At each unique squared four-momentump 2 , the right-hand side of Eq. (26) is computed for each corresponding p, for all operators in a given representation R ∈ {τ 3 }, and for all allowed choices of the Lorentz index τ of the external gluon states. Fits to these results are performed in a correlated manner to determine the RI-MOM renormalisation factor Z RI-MOM R (p 2 ) for that scale and representation. The correlations are propagated using the bootstrap resampling procedure described in Sec. III C. Choices of the number of bootstraps N boot from 200 to 1000 are tested and found to give consistent results and uncertainties. As discussed in Ref. [51], combining data from all operators in a given irreducible representation of the hypercubic group, as is done here, in general reduces the amount of O(4) violation and produces a smoother dependence of the common renormalisation factor on the scalep 2 than choosing a single operator.
In addition to the RI-MOM factors Z RI-MOM R (p 2 ) for the two representations, the complete multiplicative renor- includes a perturbative matching factor which converts from the RI-MOM renormalisation at scalep 2 to the MS scheme at µ = 2 GeV. In this work, the 1-loop expression for this matching, derived in Ref. [49], is used: For these calculations in the Landau gauge, ξ = 0, N c = 3 = N f , and g 2 is defined by α(µ MS ) evaluated to three loops [52][53][54].
The extracted renormalisation constants Z MS R (µ = 2 GeV), determined from the RI-MOM factors Z RI-MOM R (p 2 ) at a range of scales (ap) 2 , are displayed in Fig. 1. In the absence of discretisation artefacts and in the perturbative regime, each renormalisation constant would be independent of the intermediate scale (ap) 2 . It is apparent that the results obtained using Wilson-flowed fields in the gluon two-point functions have smaller discretisation artefacts than those with unflowed fields: while the former are consistent with a linear form in (ap) 2 at large scales, the latter display significant quadratic effects. Nevertheless, in the limit a → 0, the renormalisation constants constructed using flowed and unflowed gauge fields in the gluon propagators agree for each representation R. The results for the τ 8 representation display larger discretisation artefacts and statistical fluctuations than those for the τ (6) 3 representation. Values for the constants Z MS R (µ = 2 GeV) that are used to renormalise the bare lattice results for the GFFs are taken from the a = 0 extrapolations of linear fits in (ap) 2 to the renormalisation constants constructed using flowed gauge fields with different intermediate scales (ap) 2 . The significance of discretisation artefacts is assessed by taking various cuts on the four-momenta included in the fits, such that µp There is insufficient data to constrain the extrapolation for smaller X, while the fit quality decreases for larger X, indicating significant contamination from discretisation effects when more momentum components are included. For each cut, linear fits in (ap) 2 over all fit ranges with a lower bound of (ap) 2 ≥ 1 are performed. The standard deviation of the variation of the central values over all fits with acceptable χ 2 /d.o.f for all cuts is included in quadrature with the statistical uncertainty of the best-fit extrapolation to a = 0. The resulting values of the renormalisation constants for the respective representations are:

C. Matrix Elements
Bare matrix elements of the operators in Sec. III A are extracted from ratios of two-and three-point correlation functions built from quark propagators originating from an APE-smeared [55] source and having either an APEsmeared or point sink 6 . The two sets of resulting correlation functions are labelled as smeared-smeared (SS) and smeared-point (SP), respectively. For the nucleon, two-point correlation functions are defined as Here ( x 0 , t 0 ) denotes the source position, p is the chosen momentum projection, and χ α ( N + | p | 2 is the nucleon energy for a given momentum p, and Z(p) (Z(p)) controls the overlap factor of the source (sink) interpolating operator onto the nucleon state (the source and sink overlaps are distinct for the SP correlation functions). The ellipsis denotes contributions from higher excitations, which are exponentially suppressed for t f t 0 . Similarly, the pion two-point correlation function is defined by where E (π) p = m 2 π + | p | 2 and Z π (p) (Z π (p)) controls the overlap factor of the source (sink) interpolating operator onto the pion states. The interpolating operator is χ (π) ( x, t) = ψ u ( x, t)γ 5 ψ d ( x, t), constructed both with and without APE smearing as described for the nucleon.
Nucleon and pion two-point functions are evaluated for all three-momenta such that | p | 2 ≤ 5(2π/L) 2 , and for both spin components of the nucleon. Effective mass functions defined as constructed from the two-point functions of the nucleon (averaged over spins) and pion, are shown in Fig. 2 for the SP correlators. As the momentum increases, the signal quality degrades for both the nucleon and pion. Energies are extracted from constant fits to the effective masses over the longest time region with χ 2 /d.o.f ≤ 1, accounting for the correlations in the data. These time windows define the range of sink times t f where excited state contamination is small in comparison with the statistical uncertainties of the data. The energies extracted as a function of momentum using both SS and SP two-point correlators are used to construct the effective speed of light (in units of c) shown for both hadrons in Fig. 3; comparison of these quantities to unity provides a measure of discretisation errors in this calculation. On this ensemble, discretisation effects on the speed of light are at the percent level for all momenta considered. Consistent values for c 2 (N/π) were found on this ensemble in Ref. [40]. For the gluon operators O R i defined in Eqs. (17) and (18), nucleon three-point correlation functions are defined by where F R i for i = {1, 2, 3} denotes the linear combination of F µν (defined in Eq. (5)) with indices matching the structure of the corresponding operator O R i . Similarly, the three-point correlation functions of the pion are defined by where again the representation and subscript labels {R, i} correspond to the operators defined in Eqs. (17) and (18), as discussed for the nucleon. For both the nucleon and the pion, three-point functions are constructed with all possible sink three-momenta that satisfy | p | 2 ≤ 5(2π/L) 2 , with operator three-momenta | ∆| 2 ≤ 18(2π/L) 2 .
The two-and three-point correlation functions are evaluated on an average of N src = 203 randomly placed sources on each of the N cfg = 2821 configurations of the ensemble detailed in Table I. At the first stage of analysis, results on each configuration are averaged (after translation such that all sources coincide at the origin). In the discussion of further analysis, the x 0 and t 0 labels are thus omitted. A bootstrap resampling procedure over the N cfg independent samples is used to propagate the statistical uncertainties of the two-and three-point functions. In this procedure, N boot = 200 bootstrap ensembles each with N cfg elements are randomly drawn (allowing replacement). Repeating the analysis with N boot = 100 or 1000 yields consistent values and uncertainties. To test the assumption that the configurations, each separated by 10 hybrid Monte-Carlo trajectories, are independent, the analysis is also undertaken with correlation functions calculated on sets of N block = 10 successive configurations (still spaced by 10 hybrid Monte-Carlo trajectories) averaged before bootstrap resampling. This blocking process does not modify the results at a statistically significant level.
To extract the GFFs, ratios of the nucleon and pion three-point and two-point functions are formed for each of the N boot bootstrap resamplings: where the nucleon and pion energies are constructed as E (N/π) p = M 2 N/π + | p | 2 , rather than determined from the two-point functions at each three-momentum. In these ratios, the exponential time-dependence and overlap factors cancel for the ground state contribution for 0 τ t f . The decompositions of the nucleon and pion matrix elements in Eqs. (4) and (11) thus allow the ratios in this limit to be constructed in terms of only the unknown GFFs, which are functions of the squared momentum transfer t = (p − p) 2 , and known kinematic factors. At each value of the squared momentum transfer t, the various consistent choices of p, p , and operator index i for a given representation R (and spin s for the nucleon) thus provide a system of equations that can be solved to isolate the GFFs at that t. Of the large number of t values accessible using the three-momenta considered here, many are very close together (in comparison with the overall scale of momenta). To provide the best determination of the GFFs, nearby t values are thus binned as illustrated graphically in Fig. 4, and their constraints are treated as a single system of equations at each average t. The bins are defined such that no adjacent accessible t values that differ by 0.03 GeV 2 or more will be in the same bin.
To reduce the dimensionality of the highly overdetermined systems of equations that must be solved to isolate the GFFs for each momentum transfer bin, ratios that give the same linear combination of GFFs (up to a sign) in the limit 0 τ t f are averaged (including the appropriate signs). The averaged ratios are denoted R (π/N ) where the subscript k now enumerates the different averaged ratios, rather than indexing specific operators. This averaging is performed separately for ratios constructed with SS and SP correlation functions, and separately for each representation R. For the nucleon there are between 4 and 101 averaged ratios for each smearing at different momenta, and for the pion there are between 2 and 48. Subsequently, constant fits to the t f and τ dependence of the averaged ratios are performed to extract the ground state contribution. For each averaged ratio, fits are performed over windows of timeslices in the two-dimensional {t f , τ } plane with t min f < t f < t max f and ∆τ < τ < (t f − ∆ τ + 1). Fit windows are constrained to have the earliest sink time t min f no earlier than the time at which the hadron two-point functions are consistent with a single state. This constraint is imposed despite the fact that the ratios are typically noisier than the two-point functions and are consistent with a single state considerably earlier than the two-point correlators. For a given t min f , the maximum sink time t max f in the window is constrained to be at least 4 timeslices larger, and no later than than the latest time at which the hadron two-point functions are consistent with a single state (under correlated fits). The gap ∆τ is also chosen to be at least 4. Subject to these constraints, all possible constant fits are performed to a given average ratio by minimizing the total χ 2 function including both SS and SP ratios for {t f , τ } pairs in each window. The mean of the bootstrap results for the best fit is taken as the central value, while the uncertainty is formed by taking half of the variation in central values of all fits with χ 2 /d.o.f. ≤ 1 in quadrature with the standard deviation of the bootstrap results for the best fit. For the pion ratios, the variation over acceptable fits is typically small compared with the statistical uncertainty, while for the nucleon, whose signal degenerates more quickly with sink time t f , this systematic uncertainty dominates. An alternative approach to the analysis using the summation method [56] yields fits for each averaged ratio that are consistent within uncertainties.
The GFFs A g (t), B g (t), and D g (t) for the nucleon, and A

IV. RESULTS
The three gluon GFFs of the nucleon and two gluon GFFs of the pion that are extracted from the LQCD calculations detailed in Sec. III are shown in Figs. 5 and 6, respectively. Results are shown in the MS scheme at µ = 2 GeV, where the renormalisation procedure is as described in Sec. III B. The A(t) and D(t) GFFs for the nucleon and pion fall off as |t| increases and are well-described by dipole forms as well as the more general z-expansion [57] parametrisation: where k max = 2 and for both pion and nucleon t cut = 4m 2 π [57]. Higher-order z-expansion fits yield a χ 2 /d.o.f. 1, overfitting statistical fluctuations in the data. Fit parameters for each parametrisation are tabulated in Table III. Notably, the dipole masses are a factor of two larger for the pion than for the nucleon for both A(t) and D(t). Correspondingly, the gluon radii defined from either form factor are smaller for the pion than the nucleon, as is found in experiment for the respective charge radii [58]. The nucleon B g (t) GFF is consistent with zero over the entire range of t of this study.
Ratios of the GFFs are renormalisation-independent if mixing with the corresponding quark GFFs is negligible. Fig. 7 displays these ratios for both the nucleon and the pion, along with linear and quadratic fits to their t-dependence. Several notable features are apparent. First, the ratios of D (π) g /A (π) g and D g /A g are approximately linear over a wide range of t, and the nucleon ratio is larger than the pion ratio at t = 0. Second, since B q (t) ∼ 0 , the ratio (A(t) g + B(t) g )/A(t) g is approximately unity, and, correspondingly, the fractional contributions of gluons to the nucleon angular momentum and momentum are similar.
While there are no previous QCD calculations of the t-dependence of the three gluon GFFs of the nucleon, a quenched QCD calculation of A g (t) and B g (t) was presented in Ref. [59] at larger than physical quark masses. The A g (t) form factor determined in that study is considerably smaller in magnitude although with a similar tdependence to that calculated here, while the B g (t) form factor is consistent with zero. The gluon momentum fraction of the nucleon, which corresponds to the forward limit A g (0) = x g , is found in the present study to be x g (µ = 2 GeV) = 0.54 (8) at m π = 450 MeV in the MS scheme. This quantity has previously been determined in a number of LQCD calculations [31,48,[59][60][61]. These results are collated in Fig. 8(a), which also includes the phenomenological result for the momentum fraction from the CT14 PDF parametrisation [62]. It is notable that while there is some scatter in the LQCD results, likely a result of systematic uncertainties that are as-yet uncontrolled, the gluon momentum fraction is approximately constant with changes in quark masses within each study (which one could expect to have correlated systematic effects at different masses). It is expected that the GFFs at nonzero t will also be approximately independent of quark mass. The gluon momentum fraction of the pion from this work is x (π) g (µ = 2 GeV) = 0.61 (9) at m π = 450 MeV in the MS scheme. A phenomenological estimate of x (π) g ∼ 0.3 in reported in Ref. [64] but without an uncertainty. A quenched QCD calculation of this quantity has been presented in Ref. [63], and a comparison with the results of this study is shown in Fig. 8(b). As was found for the nucleon, no significant quark-mass dependence is evident in this quantity.
The quark GFFs of the nucleon and pion have been previously computed at various quark masses in LQCD using a variety of lattice actions and quark masses. Only the connected quark-line contributions have been determined in most cases. Where they have been computed at similar quark masses to those in this study, disconnected contributions are found to be at the percent level [59]. Figure 9 shows a comparison of the three gluon GFFs of the nucleon to the corresponding connected isoscalar quark GFFs computed at a similar quark mass (corresponding to m π ∼ 496 MeV) to that used in this work using domain wall valence quarks on configurations generated with an ASqTad staggered quark action at a similar lattice spacing (a = 0.125 fm) [23] to that used here. Fits using the z-parametrisation of Eq. (40) are shown to both quark and gluon of GFFs. The A g (t) and D g (t) GFFs fall off more quickly with t than the corresponding quark GFFs, and thus the generalised nucleon radii defined from the gluon GFFs are larger than those from the corresponding quark quantities. In the forward limit, A g (0) ≈ A u+d,conn. (0), indicating that quarks and gluons each contribute approximately half of the nucleon momentum at this unphysically heavy value of the quark mass (the sum of the gluon and connected quark momentum fraction is ∼ 1.07 (9), indicating that undetermined systematic effects from disconnected contributions and lattice artefacts are likely small). The quark and gluon Dterms are both negative, with D g (t) ∼ 2D u+d,(conn.) (t). The B g and B u+d,(conn.) GFFs are both consistent with zero.
The connected quark GFFs of the pion have previously been calculated using a different formulation of the clover quark action at somewhat heavier quark masses corresponding to a pion mass of 842 MeV and a lattice spacing of a = 0.073 fm [65]. These results are shown alongside the gluon GFFs of the pion computed here in Fig. 10. Fits using the z-parametrisation of Eq. (40) are shown for both quark and gluon GFFs. As for the nucleon, the gluon and quark contributions to the lightcone momentum of the pion, defined by the forward limits A g (0) and A u+d,conn. (0) respectively, are similar. Unlike for the nucleon, the corresponding pion quark and gluon GFFs A a (t) remain similar over the entire range of t that is investigated. The quark and gluon D-term GFFs in the pion are also similar in magnitude, relative to the uncertainties of each calculation.

V. CONCLUSION
In this article, the first determination of the complete set of gluon generalised gravitational form factors of the nucleon and pion from lattice QCD is presented. All GFFs are found to have dipole-like dependence on the squared momentum transfer t, with the exception of the B g (t) GFF of the nucleon that is consistent with zero over the Results from quenched QCD are also shown: the purple inverted triangles show the results of Ref. [59] (χQCD collaboration) determined using quenched QCD, the orange triangles show those from Ref. [60] (QCDSF collaboration), and the yellow filled triangles denote those from Ref. [61] (QCDSF collaboration). The experimental value for the proton is shown as the red star and is taken from the CT14 PDF parametrisation [62]. In subfigure (b), the blue squares show data from the quenched QCD calculation reported in Ref. [63].
entire rage of t that is investigated. For the nucleon, the gluon GFFs fall off faster in |t| and can be parametrised with larger dipole masses than the corresponding quark GFFs computed using similar lattice discretisations and at a similar value of the quark masses, indicating the gluon distributions have a smaller spatial size than those of the quarks. In contrast, the quark and gluon GFFs of the pion have very similar t-dependences. For both the pion and the nucleon, the gluon momentum fraction, corresponding to the forward limit of one of the GFFs, is found to be approximately 0.5-0.6, somewhat larger the phenomenological value in both cases. The gluon contributions to the nucleon momentum and angular momentum are of similar relative size.
All calculations presented here have been performed at a single lattice spacing and volume and at a single unphysical value of the light quark masses, and mixing of the isoscalar quark GFFs with the gluon GFFs has been neglected based on expectations from lattice perturbation theory [31] that these effects are small. The as-yet-unquantified systematic uncertainties that result from the lattice spacing and finite volume effects are expected to be considerably smaller than the uncertainties reported on the renormalised GFFs. Since the gluon GFFs are determined from purely gluonic operators (up to effects of mixing), the quark-mass-dependence is also expected to be mild, and extrapolation to the physical quark masses will likely not shift the GFFs outside their uncertainties. Future calculations will control these remaining systematic uncertainties and thereby allow more precise comparisons with phenomenology and also controlled predictions for the gluon contributions to the shear and pressure distributions of the nucleon and pion that are determined by the D-term GFFs [9,66].  This Appendix displays examples of raw LQCD data for the averaged ratios of three-and two-point functions R (π/N ) R;k (t f , τ ) (defined in Section III), and illustrates the results of plateau fits to the t f and τ dependence of these ratios. Figures 11-14 show data for the nucleon at a selection of values of the squared momentum transfer t. In each case, ratios are plotted vs τ at two different sink times for both SP and SS correlation functions, and vs sink time t f for two choices of operator insertion time. Also shown on each figure are both the fit band resulting from a simultaneous fit to the t f and τ dependence of the ratios formed with both SS and SP three-point functions within the plateau region (discussed in Section III), and the central value and uncertainty of the final fitted result for the GFFs at that momentum, projected back to the linear combination of that particular ratio. Analogous figures for the pion are shown in Figs. 15-18 (a later sink time is shown for the pion than for the nucleon since the signal in the statistically cleaner pion data continues to later times). Figures 19 and 20 show the constraints from the fits to each averaged ratio on the GFFs graphically at a selection of values of the squared momentum transfer t.