Lattice evidence that scalar glueballs are small

TL;DR

This study uses lattice field theory to calculate the gravitational form factors of the scalar glueball, revealing it to be significantly smaller than other hadrons, which could aid experimental identification.

Contribution

First lattice calculation of scalar glueball GFFs, showing its distinct gluonic structure and small size compared to other hadrons.

Findings

Mass radius of 0.263(31) fm for the scalar glueball

Glueballs have a different gluonic structure than typical hadrons

Supports the idea that scalar glueballs are smaller than other hadrons

Abstract

This work reports the first calculation of the gravitational form factors (GFFs) of the scalar glueball, performed via lattice field theory in Yang-Mills theory at a single lattice spacing. The glueball GFFs are compared with those of other hadrons as determined in previous lattice calculations, providing strong indications that glueballs have a different gluonic structure than typical hadronic states. A mass radius of 0.263(31) fm is predicted, supporting previous suggestions that the scalar glueball is significantly smaller than other hadrons. These results point towards a potential smoking-gun characteristic to target by experimental glueball searches.

Equations

$$\begin{array}{rl}⟨G[0^{++}](p^{\prime })|T^{\mu \nu }|G[0^{++}](p)⟩& =\\ 2P^{\mu }{P}^{\nu }A(t)+& \frac{{\Delta }^{\mu }{\Delta }^{\nu }-g^{\mu \nu }{\Delta }^2}{2}D(t),\end{array}$$

$$\begin{array}{r}{\chi }_1(x)=\frac{1}{4}\sum _{\mu \ne \nu }\text{Re Tr}{U}_{\mu \nu }^2(x),\\ {\chi }_2(x)=\frac{1}{4}\sum _{\mu \ne \nu }\text{Re Tr}{U}_{\mu \nu }^7(x),\end{array}$$

$$\chi_{i}({\mathbf{p}},t)=\sum_{{\mathbf{x}}}e^{-i{\mathbf{p}}\cdot{\mathbf{x}}}\bigg{[}\chi_{i}({\mathbf{x}},t)-\braket{\chi_{i}({\mathbf{x}},t)}\bigg{]}.$$

$$C_{ij}^{2\text{pt}}(\mathbf{p},t)=\frac{1}{T}\sum _{t_0}⟨{\chi }_i(\mathbf{p},t+t_0){\chi }_j(\mathbf{p},t_0{)}^{†}⟩$$

$${\chi }_0(\mathbf{p},t)=\sum _i{w}_{0i}(|\mathbf{p}{|}^2){\chi }_i(\mathbf{p},t).$$

$$\begin{array}{rl}{C}_{0^{++}}^{2\text{pt}}(\mathbf{p},t)& =\frac{1}{T}\sum _{t_0}⟨{\chi }_0(\mathbf{p},t+t_0){\chi }_0(\mathbf{p},t_0{)}^{†}⟩\\ & =\sum _{i,j}{w}_{0i}(|\mathbf{p}{|}^2)C_{ij}^{2\mathrm{pt}}(\mathbf{p},t)w_{0j}^{\ast }(|\mathbf{p}{|}^2).\end{array}$$

$$\begin{array}{rl}{C}_{0^{++},\mathcal{R}\ell }^{3\mathrm{pt}}({\mathbf{p}}^{\prime }& ,\mathbf{\Delta },t_s,\tau )=\\ \frac{1}{T}\sum _{t_0}& ⟨{\chi }_0({\mathbf{p}}^{\prime },t_s+t_0)T_{\mathcal{R}\ell }(\mathbf{\Delta },\tau +t_0){\chi }_0({\mathbf{p}}^{\prime },t_0{)}^{†}⟩\end{array}$$

$$\begin{array}{rl}{R}_{0^{++},\mathcal{R}\ell }(& {\mathbf{p}}^{\prime },\mathbf{\Delta },t_s,\tau )=\frac{C_{0^{++},\mathcal{R}\ell }^{3\mathrm{pt}}({\mathbf{p}}^{\prime },\mathbf{\Delta },t_s,\tau )}{C_{0^{++}}^{2\mathrm{pt}}({\mathbf{p}}^{\prime },t_s)}\\ \times & \sqrt{\frac{C_{0^{++}}^{2\mathrm{pt}}(\mathbf{p},t_s-\tau )}{C_{0^{++}}^{2\mathrm{pt}}({\mathbf{p}}^{\prime },t_s-\tau )}\frac{C_{0^{++}}^{2\mathrm{pt}}({\mathbf{p}}^{\prime },t_s)}{C_{0^{++}}^{2\mathrm{pt}}(\mathbf{p},t_s)}\frac{C_{0^{++}}^{2\mathrm{pt}}({\mathbf{p}}^{\prime },\tau )}{C_{0^{++}}^{2\mathrm{pt}}(\mathbf{p},\tau )}}\\ & \underset{\tau ,t_s-\tau \gg 0}{\to }\frac{⟨G[0^{++}](p^{\prime })|T^{\mu \nu }|G[0^{++}](p)⟩}{2\sqrt{E_{\mathbf{p}}{E}_{{\mathbf{p}}^{\prime }}}}.\end{array}$$

$$r_{\mathrm{mass},\mathrm{BF}}^2=\frac{1}{A(0)}{\left[6\frac{dA(t)}{dt}\right|}_{t=0}-\frac{3}{4m^2}(A(0)+2D(0))],$$

Full text

††preprint: MIT-CTP/5907, FERMILAB-PUB-25-0621-T, INT-PUB-25-021

Introduction: Since the inception of quantum chromodynamics (QCD) [1], glueballs have been postulated as hadronic states with purely gluonic degrees of freedom. Decades of experiments in hadron spectroscopy have collected experimental glueball candidates with a variety of allowed quantum numbers, e.g., ${J^{PC}=0^{++},0^{-+},1^{++},1^{+-},1^{--},2^{++}}$, etc. [2, 3, 4, 5, 6, 7]. Lattice QCD calculations of hadronic resonances provide valuable guidance towards the identification of putative glueball states [8, 9, 10, 11, 12]—see Refs. [13, 14] for a review. However, identifying observed hadrons as glueballs or glueball-like states based on spectroscopy remains challenging due to the mixing of hadrons with the same quantum numbers [15, 16].

Beyond the spectrum, information about the internal structure of hadrons may provide an additional pathway for classification of observed states as glueball-like or non-glueball objects. For example, features such as the momentum fraction carried by gluons may be qualitatively different for glueball and non-glueball states, and could serve as evidence of a hadron having predominantly gluonic degrees of freedom. The size of glueballs, which can be defined in various ways, can serve as another way to distinguish them from typical hadrons; model studies [17, 18, 19] have speculated that the scalar glueball might be more compact than the $\Lambda_{\text{QCD}}$ scale that is thought to define the size of more typical hadrons. This was supported by early lattice studies which, however, relied on a model-dependent Bethe-Salpeter approach to define the radius [20, 21, 22], or were performed in $SU(2)$ Yang-Mills theory [23]. The typical scales found, around 0.2 fm, are also significantly smaller than the radius of the $\sigma$ meson, which is the lightest resonance with the same quantum numbers, as predicted from a dispersive analysis [24].

The gravitational form factors (GFFs) of hadrons encode the gluon momentum fraction and distribution of quantities like the energy inside hadrons [25], offering a robust way to define their size through, e.g., their mass radius, defined from the root mean square of the energy density. The GFFs are defined from the matrix elements of the energy-momentum tensor (EMT) $T^{\mu\nu}$ of QCD [26, 25, 27, 28]. Lattice QCD provides a tool to determine these matrix elements, and thus the glueball GFFs. While previous pioneering work [29, 23] has investigated glueball matrix elements using lattice field theory, the GFFs of glueballs have not been previously constrained.

This work presents a first step towards quantitative studies of glueball structure through a calculation of the GFFs of the lowest-lying scalar glueball state111See Ref. [30] for preliminary results of this work., denoted as $G[0^{++}]$, in $SU(3)$ Yang-Mills theory. In a pure gauge theory, the EMT is purely gluonic, defined as ${T^{\mu\nu}=2\;\text{Tr}[-F^{\mu}_{\alpha}F^{\alpha\nu}+\frac{1}{4}g^{\mu\nu}F^{\alpha\beta}F_{\alpha\beta}]}$, where $F^{\mu\nu}$ is the gluon field strength tensor. The EMT matrix element for the scalar glueball can be decomposed in terms of two GFFs, $A(t)$ and $D(t)$ [31, 32, 33]:

[TABLE]

where the four-momenta of the incoming and outgoing states are denoted $p$ and $p^{\prime}$, $P=(p+p^{\prime})/2$, $\Delta=p^{\prime}-p$, and $t=\Delta^{2}$. Note that $A(0)=1$ follows from the momentum sum rule, while the $D$-term, $D(0)$, is an unconstrained and previously unknown quantity.

Methodology: Calculations are performed using a single pure-gauge ensemble defined by the $SU(3)$ Wilson gauge action with $\beta=5.95$. This results in $a=0.098$ fm, using the Sommer parameter to set the scale [34, 35]. The lattice geometry is $L^{3}\times T=24^{3}\times 48$. Calculations are performed using $\mathcal{O}(10^{7})$ configurations generated using $\mathcal{O}(10^{5})$ independent streams of heatbath with overrelaxation [36, 37, 38, 39, 40], which were saved every 25 heatbath hits. Two-point correlation functions are constructed on each configuration using two interpolating operators:

[TABLE]

where $\mu,\nu\in\{x,y,z\}$. Here, $U^{n}_{\mu\nu}$ is an $n\times n$ Wilson loop which extends in the $\mu$ and $\nu$ directions constructed from links stout-smeared [41] by 3 steps in the spatial directions only. It is convenient to define vacuum-subtracted operators of definite three-momentum ${\mathbf{p}}$ as:

[TABLE]

The positive-parity $0^{++}$ glueball is the lowest-energy state excited by these operators; the summations over $\mu,\nu$ in Eq. (2) project to the $A_{1}^{+}$ (rest frame) or $A_{1}$ (moving frames) irreducible representations (irreps) of the finite-volume symmetry group, while taking the real part projects to positive charge conjugation quantum numbers. These irreps have maximal overlap with the scalar glueball; however, depending on the momentum frame, the spectrum may also include heavier glueballs with other quantum numbers, e.g. tensor or pseudoscalar glueballs, or multi-glueball or higher-lying ditorelon states [42].

The spectrum of glueball states is constrained via the generalized eigenvalue problem (GEVP) method [43, 44, 45, 46, 47, 48, 49]. The GEVP is applied to the $2\times 2$ matrix of momentum-projected two-point functions averaged over all timeslices:

[TABLE]

for all $\left|{\mathbf{p}}\right|^{2}\leq 6(2\pi/L)^{2}$ on 200 bootstrap ensembles after binning the $\mathcal{O}(10^{7})$ configurations into groups of $1000$ consecutive configurations. Averaging over equivalent momenta yields two-point functions for the 7 distinct $|{\mathbf{p}}|^{2}$; analyzing them in a “fixed pivot” [43, 44, 45] mode with $t_{0}=1$ and diagonalization time $t_{d}=3$ gives 7 sets of weights $w_{ij}{(|{\mathbf{p}}|^{2})}$ used to construct optimized ground-state interpolating operators:

[TABLE]

These are used to construct two-point correlation functions as:

[TABLE]

For each $\left|{\mathbf{p}}\right|^{2}$, the generalized eigenvalue problem is then solved to extract the ground state, which is identified as the scalar glueball. Figures showing examples of the two-point correlation functions, and the corresponding analysis, are included in the Supplementary Material (SM).

The matrix elements defined in Eq. (1) are obtained by using GEVP-optimized interpolating operators $\chi_{0}$ to compute vacuum-subtracted three-point functions [46]

[TABLE]

for $\left|{\mathbf{\Delta}}\right|^{2}\leq 10(2\pi/L)^{2}$ and $\left|{\mathbf{p}}^{\prime}\right|^{2}\leq 6(2\pi/L)^{2}$. $T_{\mathcal{R}\ell}$ denotes the vacuum-subtracted gluon EMT projected to the $\ell$ row of the hypercubic group irreps [53] $\mathcal{R}\in\{\tau_{1}^{(3)},\tau_{3}^{(6)}\}$. The individual EMT components are computed using the clover definition of the gluon field strength tensor, constructed using gluon links that have been acted on with 2 or 3 steps of stout smearing [41]. The two- and three-point correlation functions are then combined into ratios:

[TABLE]

Ratios with choices of momenta, $\mathcal{R}$, and $\ell$ that result in the same linear combination of GFFs as defined in Eq. (1) are averaged. These averaged ratios are fit to extract the ground-state contribution. The fitting procedure is cross-checked against recent generalizations of the Lanczos method [54, 55, 56, 57], which provides statistically compatible results. Additional details, which closely follow previous lattice QCD studies of GFFs, e.g., Refs. [58, 59, 60, 61, 51, 52], and figures showing the results of each analysis step, are given in the SM.

Finally, the GFFs $A(t)$ and $D(t)$ are obtained by first grouping the fit results for the averaged ratios $\bar{R}_{0^{++},\mathcal{R}\ell}$ for each irrep $\mathcal{R}$ separately into 14 bins using k-means clustering [62] on the momentum transfer squared $t=\Delta^{2}$, then solving the overconstrained systems of linear equations dictated by Eq. (1) to obtain the bare GFFs $A_{\mathcal{R}}(t)$ and $D_{\mathcal{R}}(t)$ for each irrep and bin. They can be renormalized by imposing the sum rule $A(0)=1$. The renormalization factors $1/A^{\mathrm{bare}}_{\mathcal{R}}(0)$ are obtained from a fit of a $n$-pole model $\alpha/(1-t/\Lambda^{2})^{n}$ to the bare GFF $A_{\mathcal{R}}(t)$; $\alpha$ and $\Lambda$ are fitted parameters, and $\alpha=A^{\mathrm{bare}}_{\mathcal{R}}(0)$. The bare GFFs in each momentum bin for each irrep are then multiplied by the renormalization factors, and then averaged together. Subsequently, the final results are obtained from a fit to the $n$-pole form, with $\alpha$ set to $1$ to enforce $A(0)=1$. Various analysis choices do not affect the results, as shown in the SM. Specifically, the choices $n\in\{1,2,3\}$ yield consistent values for the GFFs; in the main text $n=3$. Similarly, separate analyses using the bare three-point functions with 2 and 3 steps of stout smearing are consistent; in the main text 2 steps of stout smearing are used, as the $n$-pole fits to this data have a significantly larger p-value.

Results: Figure 1 shows $A(t)$ and $D(t)$ of the $0^{++}$ glueball determined in this work. Also shown for comparison are the gluon GFFs, up to a factor222Mesons and baryons with valence quarks and antiquarks also receive a quark contribution to their GFFs, which mixes with the gluonic component. The work of Ref. [50] constrained only the gluon contribution of the GFFs, also neglecting the mixing. The comparison of the overall normalization between those results and the gluon GFFs in this work—which coincide with the total GFFs in a theory with only gluonic degrees of freedom, as investigated here—is not meaningful. of $A(0)$, of four hadrons with quantum numbers $J^{P}=0^{-}$, $1^{-}$, $1/2^{+}$, and $3/2^{+}$, corresponding to the pion, $\rho$ meson, nucleon, and $\Delta$ baryon, computed using a single lattice QCD ensemble with $N_{f}=2+1$ clover-improved dynamical quark flavors, $a\simeq 0.12~\text{fm}$, and $m_{\pi}\simeq 450~\text{MeV}$ [50]. While the present calculation is undertaken in Yang-Mills theory, this qualitative comparison may nevertheless yield interesting conclusions; a quenched calculation of the gluon GFFs of a pseudoscalar meson state with valence quark masses tuned such that the hadron mass matches that of Ref. [50] (shown in the SM) is consistent with the QCD results for the pion shown in Fig. 1, indicating that quenching has a negligible effect on these quantities at the level of uncertainties of this calculation. The glueball $A(t)$ form factor decays more slowly than that of the pion, corresponding to a smaller mass radius contribution. The same is observed in the respective $D$ form factor, although the statistical uncertainty of $D(t)$ is larger for the glueball333The flat $t$-dependence of $D(t)$ is difficult to interpret at this level of statistics, especially when taking into consideration the 3-stout smeared result shown in the SM, where the behavior is more consistent with monotonically increasing..

Fig. 2 shows a comparison of the results of this work with the total GFFs of the pion [52] and the nucleon [51], computed using an ensemble with the same action as in Ref. [50] but smaller lattice spacing $a\simeq 0.091~\mathrm{fm}$ and closer to the physical point with $m_{\pi}\simeq 170~\text{MeV}$. Since this comparison shows total GFFs, no division by $A(0)$ is necessary. The same qualitative conclusions are drawn; the glueball GFFs decay more slowly than those of the other hadrons, implying a smaller radius.

To carefully assess the radii comparison, we consider the root mean square radius of the energy density of the glueball in the Breit frame, also known as the mass radius. It is [27, 63]

[TABLE]

and is found to be $0.263(31)~\text{fm}$. This is remarkably small, and close to older lattice QCD and Bethe-Salpeter amplitude findings, which estimated the size of the scalar glueball to be in the range $0.1-0.3~\text{fm}$ [20, 21, 22, 17, 18] by analyzing its (Bethe-Salpeter) wavefunction. This supports the view that the scalar glueball size is set by short-range gluon interactions (e.g. instanton-induced forces [17, 18]) and not the confinement scale of QCD. It is also smaller than a dispersive prediction of the scalar radius of the $\sigma$ meson [24], suggesting that the scalar glueball is the most compact state in that channel.

Figure 3 shows a comparison of this quantity with the results of the $n$-pole fits of Refs. [50, 51] for the other hadron states444The numbers in Fig. 3 taken from Ref. [50] differ from those currently in Table X of the published article. This is due to an error in normalization in the results of that table only; an erratum is forthcoming.555We do not compare with the radius of the $m_{\pi}\simeq 170$ MeV pion study of Ref. [52], because for light hadrons the 3D BF density interpretation is expected to break down [64, 65, 66, 67, 68, 69, 26, 70].. Note that other than for the pion, the BF mass radius of the rest of the hadrons depends on additional GFFs besides $A(t)$ and $D(t)$ [27, 71, 72, 73], and their contribution is included in the radii values of Fig. 3666This excludes the $D_{1,g}^{\rho}$ form factor for $\rho$ and the $F_{11}^{\Delta,g}$, $F_{21}^{\Delta,g}$, $F_{41}^{\Delta,g}$, and $F_{51}^{\Delta,g}$ form factors for $\Delta$ that were consistent with zero and were not fit with a tripole form in Ref. [50].. The scalar glueball radius is considerably smaller than that of the other hadrons.

Conclusion: These results constitute the first investigation of the internal structure of glueballs in an SU(3) lattice gauge theory, marking a promising advance in the understanding of the internal structure of potential glueball-like hadrons in nature.

While the calculation is undertaken in Yang-Mills theory instead of QCD, comparison with (quenched and full) QCD calculations of the (gluon and total) GFFs of other hadrons provides the first indication that glueballs may have a different gluonic structure than more typical hadronic states; specifically, the $A(t)$ GFF decays more slowly, corresponding to a smaller gluonic radius. This observation, if it is maintained in future QCD calculations of glueball GFFs, would open the potential for future smoking gun identification of glueball states directly via their gluon structure.

It would also be interesting to repeat this study with other glueball states in Yang-Mills theory, for instance with pseudoscalar or tensor quantum numbers. Notably, early lattice QCD studies based on wavefunction analyses suggested that the tensor glueball has a significantly larger radius than the scalar [20, 21, 22]. Future calculations could benefit from recent multi-level algorithm development [74, 75].

Acknowledgments: The authors acknowledge the MIT SuperCloud [76] and Lincoln Laboratory Supercomputing Center for providing HPC resources that have contributed to the research results reported within this paper. This work made use of resources provided by subMIT at MIT Physics [77].

The authors thank Julian Urban for contributions at early stages of the project. We also thank Lorenzo Barca, Max Hansen, and André Walker-Loud for useful discussions. PES is supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under grant Contract Number DE-SC0011090, by Early Career Award DE-SC0021006, by Simons Foundation grant 994314 (Simons Collaboration on Confinement and QCD Strings), by the U.S. Department of Energy SciDAC5 award DE-SC0023116, and has benefited from the QGT Topical Collaboration DE-SC0023646. DAP is supported from the Office of Nuclear Physics, Department of Energy, under contract DE-SC0004658. This document was prepared using the resources of the Fermi National Accelerator Laboratory (Fermilab), a U.S. Department of Energy, Office of Science, Office of High Energy Physics HEP User Facility. Fermilab is managed by Fermi Forward Discovery Group, LLC, acting under Contract No. 89243024CSC000002. DAP thanks the Albert Einstein Center at the University of Bern for its hospitality during a visit that significantly advanced this work. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. RA is supported by the U.S. Department of Energy SciDAC5 award DE-SC0023116 and the High Energy Physics Computing Traineeship for Lattice Gauge Theory (DE-SC0024053). PES thanks the Institute for Nuclear Theory at the University of Washington for its kind hospitality and stimulating research environment. This research was supported in part by the INT’s U.S. Department of Energy grant No. DE-FG02-00ER41132. This manuscript was finalized at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-2210452.

The Chroma [78] and LALIBE [79] and JAX [80] libraries were used in this work. Data analysis used NumPy [81], SciPy [82], lsqfit [83], and gvar [84]. Figures were produced using matplotlib [85].

I Additional results

GEVP analysis: Figure 4 compares the ground-state energies extracted using GEVP-optimized interpolators with the individual interpolators in the basis.

Matrix element fits: The matrix elements are extracted via a constant fit to all possible ($t_{s}$,$\tau$) time ranges of the ratios of Eq. (8), after having averaged over all kinematics and irrep vectors that correspond to identical linear combinations of GFFs. The minimum number of data points per fit is set to 8, and additionally $t_{s_{\text{min}}}\geq 5$, $t_{s_{\text{max}}}\leq 13$, $\tau\geq 2$, $\tau\leq t_{s}-2$ are imposed. The final matrix element results are averages of all fits weighted by their $p$-value weights [86, 87, 88]. Examples of ratios and corresponding fits are shown in Fig 5.

Fig. 6 shows the results of a cross-check of our fitting procedure against a block Lanczos matrix element analysis [54, 55, 56, 57]. The Lanczos approach requires a three-point function and the corresponding initial- and final-state two-point functions as inputs, and so is inapplicable to the binned ratios fit for the primary analysis. Data were not available for a full reanalysis using Lanczos methods. Thus, for the test, we repeat the ratio fitting procedure on a small subset of unbinned ratios with definite kinematics. For the Lanczos analysis of the corresponding two- and three-point functions, we use Hermitian subspace [55] and ZCW [56] filtering with $\epsilon^{\mathrm{ZCW}}=10^{-3}$, and use a nested bootstrap median estimator to compute values and errors. The resulting estimates shown in Fig 6 converge to a stable value after a few iterations. In all cases we have examined, the value provided is statistically consistent with the fit result; the uncertainties on fit and Lanczos estimates are comparable.

Additional GFF results: Figure 7 shows a comparison of the GFFs using 3 versus 2 steps of stout smearing for the EMT. In Fig. 8, we compare different choices of the $n$-pole model ($n=1,2,3$) for the 2-stout smeared operator.

II Pseudoscalar meson analysis

Figure 1 of the main text presents a comparison of the results of this work against a previous extraction of the gluon GFFs of four hadrons of different spin [50]. The older calculation was performed with a Lüscher-Weisz gauge action [89], and $N_{f}=2+1$ flavors of non-perturbatively improved clover Wilson quarks [90], with their bare masses tuned yield a pion mass of $m_{\pi}\simeq 450~\text{MeV}$. Since for the glueball in Yang-Mills the overall normalization of $A_{g}(t)$ is trivial and set as $A_{g}(0)=1$, the glueball GFFs are compared with $A_{g}(t)/A_{g}(0)$ and $D_{g}(t)/A_{g}(0)$ of the other hadrons. However, given the difference in actions used between the two calculations, it is important to assess to what extent the comparison of the rest of the properties—the $t$-dependence of $A(t)$ and $D(t)$, as well as $D(0)/A(0)$—is physically meaningful. For that purpose, a quenched calculation of the pion gluon GFFs is undertaken using the same Yang-Mills action as was used for the glueball calculation, and clover-improved Wilson valence quarks with clover coefficient set to 1, and bare mass tuned to match the pion mass of Ref. [50], $m_{\pi}\simeq 450~\text{MeV}$. Mixing between gluon and quark GFFs is neglected, as was done in Ref. [50].

A single interpolating operator $\overline{u}(x)\gamma_{5}d(x)$ is used, with quark fields that have been smeared by Gaussian gauge-invariant smearing with radius $4.7a$. The two-point correlator is measured at 12 spacetime sources of 2760 configurations, and momenta $|{\mathbf{p}}|^{2}\leq 6(2\pi/L)^{2}$. A GEVP analysis using a $2\times 2$ matrix of two-point correlation functions constructed using both smeared-smeared and point-smeared propagators yields consistent results as the non-GEVP smeared-smeared only analysis, so the latter is taken as the primary result. The two-point correlator is multiplied with the EMT operator on each of the 2760 configurations, using an EMT with 2 and 3 steps of stout smearing and momenta $\left|{\mathbf{\Delta}}\right|^{2}\leq 10(2\pi/L)^{2}$. 200 bootstrap ensembles of vacuum substracted three-point functions are taken to form the ratios of Eq. (8). The rest of the analysis proceeds identically to the glueball GFF extraction.

The final results for $A_{g}(t)/A_{g}(0)$ and $D_{g}(t)/A_{g}(0)$ of the quenched pion are shown in Fig. 9, and compared against the pion GFFs of Ref. [50]. The two sets of results are in good agreement.