Off-forward gluonic structure of vector mesons

TL;DR

This study uses lattice QCD to explore the gluonic structure of the $$ meson, revealing differences in spatial distributions of gluons and quarks, advancing understanding of hadronic internal structure.

Contribution

First lattice QCD calculation of the gluonic GFFs of the $$ meson, providing new insights into gluon distributions and their spatial characteristics.

Findings

Quantitative differences between gluon and quark distributions.

First determination of transversity gluon GFFs for the $$ meson.

Comparison with quark GFFs and electromagnetic form factors.

Abstract

The spin-independent and transversity generalised form factors (GFFs) of the $\phi$ meson are studied using lattice QCD calculations with light quark masses corresponding to a pion mass $m_\pi\sim450(5)$ MeV. One transversity and three spin-independent GFFs related to the lowest moments of leading-twist spin-independent and transversity gluon distributions are obtained at six non-zero values of the momentum transfer up to 1.2 GeV$^2$. These quantities are compared with the analogous spin-independent quark GFFs and the electromagnetic form factors determined on the same lattice ensemble. The results show quantitative distinction between the spatial distribution of transversely polarised gluons, unpolarised gluons, and quarks, and point the way towards further investigations of the gluon structure of nucleons and nuclei.

Figures

Tables

Table 1. Table 1: LQCD simulation details. The gauge configurations have dimensions L 3 × T superscript 𝐿 3 𝑇 L^{3}\times T , lattice spacing a 𝑎 a , and bare quark masses a ​ m q 𝑎 subscript 𝑚 𝑞 am_{q} (in lattice units). A total of N meas subscript 𝑁 meas N_{\textrm{meas}} light-quark sources were used to perform measurements across N cfg subscript 𝑁 cfg N_{\textrm{cfg}} configurations.

$L/a$	$T/a$	$\beta$	$a{m}_l$	$a{m}_s$	$a$ (fm)	$L$ (fm)	$T$ (fm)	$m_{\pi }$ (MeV)	$m_K$ (MeV)	$m_{\varphi }$ (MeV)	$m_{\pi }L$	$m_{\pi }T$	$N_{\text{cfg}}$	$N_{\text{meas}}$
24	64	6.1	-0.2800	-0.2450	0.1167(16)	2.801(29)	7.469(77)	450(5)	596(6)	1040(3)	6.390	17.04	1042	${10}^5$

Table 2. Table 2: One choice of the initial and final three-momenta (in lattice units) and polarisations (in a cartesian basis) { p → , p → ′ , E , E ′ } → 𝑝 superscript → 𝑝 ′ 𝐸 superscript 𝐸 ′ \{\vec{p},\vec{p}\,^{\prime},E,E^{\prime}\} contributing to each reduced set of averaged ratios R ¯ ¯ 𝑅 \overline{R} . The ordering of the rows here is the same as in Eq. ( 35 ).

$\overrightarrow{p}$	${\overrightarrow{p}}^{\prime }$	$E$	$E^{\prime }$
[1, 0, 0]	[1, -1, 0]	3	3
[0, 1, 0]	[0, 1, -1]	3	3
[0, 0, 0]	[0, -1, 0]	3	3
[0, 0, 1]	[0, -1, 1]	3	3
[0, 0, 0]	[0, 0, -1]	3	3
[1, 0, 0]	[1, -1, 0]	2	2
[0, 1, 0]	[0, 1, -1]	1	1
[0, 0, 1]	[0, -1, 1]	1	1
[0, 0, 0]	[0, -1, 0]	1	1
[1, 0, 0]	[1, -1, 0]	1	1
[0, 0, 1]	[0, -1, 1]	2	2
[0, 0, 0]	[0, -1, 0]	2	2
[0, 1, 0]	[0, 1, -1]	2	2
[0, 0, 0]	[0, 0, -1]	2	2
[1, 0, 0]	[1, -1, 0]	2	1
[0, 0, 1]	[0, -1, 1]	2	3
[0, 1, 0]	[0, 1, -1]	2	3
[0, 0, 1]	[0, -1, 1]	3	2
[1, 0, 0]	[1, -1, 0]	1	2
[0, 1, 0]	[0, 1, -1]	3	2

Table 3. Table 3: One choice of ratios of two and three-point functions that allow extraction of the electromagnetic form factors. Here, E = E ​ ( Δ 2 ) 𝐸 𝐸 superscript Δ 2 E=E(\Delta^{2}) and p 𝑝 p denotes one unit of momentum in lattice units. For each R X , a subscript 𝑅 𝑋 𝑎 R_{X,a} , rotationally equivalent contributions are averaged.

$\overrightarrow{\Delta }/p$	$X$	${ℳ}_X$	$N_{X,a}$	$R_{X,a}$	$N_{X,b}$	$R_{X,b}$	$N_{X,c}$	$R_{X,c}$
	C	$1/3$	1	${ℛ}_{11}^4$	1	${ℛ}_{22}^4$	1	${ℛ}_{33}^4$
(0,0,0)	M	–	–	–	–	–	–	–
	Q	–	–	–	–	–	–	–
	C	$\frac{2\sqrt{Em}}{3\left(E+m\right)}$	1	${ℛ}_{11}^4$	2	${ℛ}_{22}^4$	–	–
(1,0,0)	M	$-\frac{2\sqrt{E^3}}{p\sqrt{m}}$	–	–	–	–	1	${ℛ}_{12}^2$
	Q	$\frac{2\sqrt{E{m}^3}}{p^2}$	1	${ℛ}_{11}^4$	-1	${ℛ}_{22}^4$	–	–
	C	$\frac{2\sqrt{m}}{3\sqrt{E}\left(E+m\right)}$	$\sqrt{2}\sqrt{E^2+m^2}$	${ℛ}_{11}^4$	$\left(2E-m\right)$	${ℛ}_{33}^4$	–	–
(1,1,0)	M	$-\frac{2\sqrt{E\left(E^2+m^2\right)}}{p\sqrt{m}}$	–	–	–	–	1	${ℛ}_{13}^3$
	Q	$\frac{\sqrt{m^3}}{p^2\sqrt{E}}$	$\sqrt{2}\sqrt{E^2+m^2}$	${ℛ}_{11}^4$	$-\left(E+m\right)$	${ℛ}_{33}^4$	–	–
	C	$\frac{2\sqrt{E^2+2m^2}}{3\sqrt{Em}\left(E+m\right)}$	$\frac{\sqrt{3}\left(2E+m\right)}{3}$	${ℛ}_{11}^4$	$-\frac{4\sqrt{3}p\left(E+2m\right)}{3\left(E+m\right)}$	${ℛ}_{12}^1$	$-\frac{4\sqrt{3}p\left(2E+m\right)}{3\left(E+m\right)}$	${ℛ}_{12}^2$
(1,1,1)	M	$\frac{2\sqrt{Em}\sqrt{2m^2+E^2}}{\sqrt{3}p\left(E+m\right)}$	–	–	1	${ℛ}_{12}^1$	-1	${ℛ}_{12}^2$
	Q	$\frac{2\sqrt{m}\sqrt{E^2+2m^2}}{3p^2\sqrt{3E}}$	$\left(m-E\right)$	${ℛ}_{11}^4$	$\frac{{\left(E+2m\right)}^2}{3p}$	${ℛ}_{12}^1$	$\frac{\left(E+2m\right)\left(2E+m\right)}{6p}$	${ℛ}_{12}^2$
	C	$\frac{2\sqrt{Em}}{3\left(E+m\right)}$	$1$	${ℛ}_{11}^4$	$2$	${ℛ}_{22}^4$	–	–
(2,0,0)	M	$-\frac{\sqrt{E^3}}{p\sqrt{m}}$	–	–	–	–	1	${ℛ}_{12}^2$
	Q	$\frac{\sqrt{E{m}^3}}{2p^2}$	$1$	${ℛ}_{11}^4$	$-1$	${ℛ}_{22}^4$	–	–
	C	$\frac{2\sqrt{m}}{3\sqrt{E}\left(E+m\right)}$	$\sqrt{5}\sqrt{E^2+4m^2}$	${ℛ}_{11}^4$	$2\left(E-2m\right)$	${ℛ}_{33}^4$	–	–
(1,2,0)	M	$-\frac{2\sqrt{E}\sqrt{E^2+4m^2}}{p\sqrt{5m}}$	–	–	–	–	1	${ℛ}_{13}^3$
	Q	$\frac{2\sqrt{m^3}}{5p^2\sqrt{E}}$	$\sqrt{5}\sqrt{E^2+4m^2}$	${ℛ}_{11}^4$	$-\left(E+4m\right)$	${ℛ}_{33}^4$	–	–
	C	$\frac{2\sqrt{2E^2+m^2}}{9\sqrt{3Em}\left(E+m\right)}$	$\left(5m-2E\right)$	${ℛ}_{11}^4$	$\frac{\sqrt{2}\sqrt{E^2+5m^2}\left(4E-m\right)}{\sqrt{2E^2+m^2}}$	${ℛ}_{22}^4$	–	–
(2,1,1)	M	$\frac{\sqrt{2Em}\sqrt{m^2+2E^2}}{\left(E+m\right)\left(E+5m\right)}$	$\frac{2\sqrt{6}}{3}$	${ℛ}_{11}^4$	$-\frac{2\sqrt{5m^2+E^2}}{\sqrt{m^2+2E^2}}$	${ℛ}_{22}^4$	$-\frac{\sqrt{6}\left(E+m\right)}{p}$	${ℛ}_{12}^1)$
	Q	$\frac{\sqrt{m}\sqrt{2E^2+m^2}}{9\sqrt{3E}{p}^2}$	$\left(E+5m\right)$	${ℛ}_{11}^4$	$-\frac{\sqrt{2}\sqrt{E^2+5m^2}\left(2E+m\right)}{\sqrt{2E^2+m^2}}$	${ℛ}_{22}^4$	–	–

Equations

$${\overline{\mathcal{O}}}_{\mu \nu {\mu }_1\dots {\mu }_n}$$

$${\tilde{\mathcal{O}}}_{\mu \nu {\mu }_1\dots {\mu }_n}$$

$${\mathcal{O}}_{\mu \nu {\mu }_1\dots {\mu }_n}$$

$$E^{\mu }(\vec{p},\lambda )=\left(\frac{\vec{p}\cdot {\vec{e}}_{\lambda }}{m},{\vec{e}}_{\lambda }+\frac{\vec{p}\cdot {\vec{e}}_{\lambda }}{m(m+\mathcal{E})}\vec{p}\right),$$

$${\vec{e}}_{\pm }$$

$${\vec{e}}_0$$

$$⟨p^{\prime }{E}^{\prime }|S[G_{\mu \alpha }i{\overleftrightarrow{D}}_{{\mu }_1}$$

$$=\sum _{\begin{array}{c}m\text{ even}\\ m=0\end{array}}^n\{$$

$$+B_{2,m}^{(n+2)}({\Delta }^2)S\left[(E\cdot E^{\prime \ast })P_{\mu }{P}_{\nu }{\Delta }_{{\mu }_1}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}\right]$$

$$+B_{3,m}^{(n+2)}({\Delta }^2)S\left[(E\cdot E^{\prime \ast }){\Delta }_{\mu }{\Delta }_{\nu }{\Delta }_{{\mu }_1}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}\right]$$

$$+B_{4,m}^{(n+2)}({\Delta }^2)S\left[\left((E^{\prime \ast }\cdot P)E_{\mu }+(E\cdot P)E_{\mu }^{\prime \ast }\right)P_{\nu }{\Delta }_{{\mu }_1}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}\right]$$

$$+B_{5,m}^{(n+2)}({\Delta }^2)S\left[\left((E^{\prime \ast }\cdot P)E_{\mu }-(E\cdot P)E_{\mu }^{\prime \ast }\right){\Delta }_{\nu }{\Delta }_{{\mu }_1}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}\right]$$

$$+\frac{B_{6,m}^{(n+2)}({\Delta }^2)}{M^2}S\left[(E\cdot P)(E^{\prime \ast }\cdot P)P_{\mu }{P}_{\nu }{\Delta }_{{\mu }_1}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}\right]$$

$$+\frac{B_{7,m}^{(n+2)}({\Delta }^2)}{M^2}S\left[(E\cdot P)(E^{\prime \ast }\cdot P){\Delta }_{\mu }{\Delta }_{\nu }{\Delta }_{{\mu }_1}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}\right]\}$$

$$+\mathrm{mod}(n,2)B_{4,n}^{(n+2)}((E^{\prime \ast }\cdot P)E_{\mu }+(E\cdot P)E_{\mu }^{\prime \ast }){\Delta }_{\nu }{\Delta }_{{\mu }_1}\dots {\Delta }_{{\mu }_n}.$$

$$⟨p^{\prime }{E}^{\prime }|S[G_{\mu {\mu }_1}$$

$$=\sum _{\begin{array}{c}m\text{ even}\\ m=2\end{array}}^n\{$$

$$+A_{2,m-2}^{(n)}({\Delta }^2)S\left[({\Delta }_{\mu }{E}_{{\mu }_1}-E_{\mu }{\Delta }_{{\mu }_1})({\Delta }_{\nu }{E}_{{\mu }_2}^{\prime \ast }-E_{\nu }^{\prime \ast }{\Delta }_{{\mu }_2}){\Delta }_{{\mu }_3}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}\right]$$

$$+A_{3,m-2}^{(n)}({\Delta }^2)S[\left(({\Delta }_{\mu }{E}_{{\mu }_1}-E_{\mu }{\Delta }_{{\mu }_1})(P_{\nu }{E}_{{\mu }_2}^{\prime \ast }-E_{\nu }^{\prime \ast }{P}_{{\mu }_2})-({\Delta }_{\mu }{E}_{{\mu }_1}^{\prime \ast }-E_{\mu }^{\prime \ast }{\Delta }_{{\mu }_1})(P_{\nu }{E}_{{\mu }_2}-E_{\nu }{P}_{{\mu }_2})\right)$$

$$\times {\Delta }_{{\mu }_3}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}]$$

$$+A_{4,m-2}^{(n)}({\Delta }^2)S\left[(E_{\mu }{E}_{{\mu }_1}^{\prime \ast }-E_{{\mu }_1}{E}_{\mu }^{\prime \ast })(P_{\nu }{\Delta }_{{\mu }_2}-P_{{\mu }_2}{\Delta }_{\nu }){\Delta }_{{\mu }_3}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}\right]$$

$$+\frac{A_{5,m-2}^{(n)}({\Delta }^2)}{M^2}S[((E\cdot P)(P_{\mu }{\Delta }_{{\mu }_1}-{\Delta }_{\mu }{P}_{{\mu }_1})({\Delta }_{\nu }{E}_{{\mu }_2}^{\prime \ast }-E_{\nu }^{\prime \ast }{\Delta }_{{\mu }_2})$$

$$+(E^{\prime \ast }\cdot P)(P_{\mu }{\Delta }_{{\mu }_1}-{\Delta }_{\mu }{P}_{{\mu }_1})({\Delta }_{\nu }{E}_{{\mu }_2}-E_{\nu }{\Delta }_{{\mu }_2})){\Delta }_{{\mu }_3}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}]$$

$$+\frac{A_{6,m-2}^{(n)}({\Delta }^2)}{M^2}S[((E\cdot P)(P_{\mu }{\Delta }_{{\mu }_1}-{\Delta }_{\mu }{P}_{{\mu }_1})(P_{\nu }{E}_{{\mu }_2}^{\prime \ast }-E_{\nu }^{\prime \ast }{P}_{{\mu }_2})$$

$$-(E^{\prime \ast }\cdot P)(P_{\mu }{\Delta }_{{\mu }_1}-{\Delta }_{\mu }{P}_{{\mu }_1})(P_{\nu }{E}_{{\mu }_2}-E_{\nu }{P}_{{\mu }_2})){\Delta }_{{\mu }_3}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}]$$

$$+\frac{A_{7,m-2}^{(n)}({\Delta }^2)}{M^2}(E^{\prime \ast }\cdot E)S\left[(P_{\mu }{\Delta }_{{\mu }_1}-{\Delta }_{\mu }{P}_{{\mu }_1})(P_{\nu }{\Delta }_{{\mu }_2}-{\Delta }_{\nu }{P}_{{\mu }_2}){\Delta }_{{\mu }_3}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}\right]$$

$$+\frac{A_{8,m-2}^{(n)}({\Delta }^2)}{M^4}(E\cdot P)(E^{\prime \ast }\cdot P)S\left[(P_{\mu }{\Delta }_{{\mu }_1}-{\Delta }_{\mu }{P}_{{\mu }_1})(P_{\nu }{\Delta }_{{\mu }_2}-{\Delta }_{\nu }{P}_{{\mu }_2}){\Delta }_{{\mu }_3}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}\right]\}$$

$$+\sum _{\begin{array}{c}m\text{ odd}\\ m=3\end{array}}^n\{A_{9,m-2}^{(n)}({\Delta }^2)S[(P_{\mu }{\Delta }_{{\mu }_1}-{\Delta }_{\mu }{P}_{{\mu }_1})(E_{\nu }^{\prime \ast }{P}_{{\mu }_2}{E}_{{\mu }_3}+E_{\nu }{P}_{{\mu }_2}{E}_{{\mu }_3}^{\prime \ast }-E_{{\mu }_2}^{\prime \ast }{P}_{\nu }{E}_{{\mu }_3}-E_{{\mu }_2}{P}_{\nu }{E}_{{\mu }_3}^{\prime \ast })$$

$$\times {\Delta }_{{\mu }_4}\dots {\Delta }_{{\mu }_{m+1}}{P}_{{\mu }_{m+2}}\dots P_{{\mu }_n}]$$

$$+\frac{A_{10,m-2}^{(n)}({\Delta }^2)}{M^2}[(P_{\mu }{\Delta }_{{\mu }_1}-{\Delta }_{\mu }{P}_{{\mu }_1})(E_{\nu }^{\prime \ast }{\Delta }_{{\mu }_2}E\cdot P-E_{\nu }{\Delta }_{{\mu }_2}{E}^{\prime \ast }\cdot P-E_{{\mu }_2}^{\prime \ast }{\Delta }_{\nu }E\cdot P+E_{{\mu }_2}{\Delta }_{\nu }{E}^{\prime \ast }\cdot P)$$

$$\times {\Delta }_{{\mu }_3}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}]$$

$$+\frac{A_{11,m-2}^{(n)}({\Delta }^2)}{M^2}[(P_{\mu }{\Delta }_{{\mu }_1}-{\Delta }_{\mu }{P}_{{\mu }_1})(E_{\nu }^{\prime \ast }{P}_{{\mu }_2}E\cdot P+E_{\nu }{P}_{{\mu }_2}{E}^{\prime \ast }\cdot P-E_{{\mu }_2}^{\prime \ast }{P}_{\nu }E\cdot P-E_{{\mu }_2}{P}_{\nu }{E}^{\prime \ast }\cdot P)$$

$$\times {\Delta }_{{\mu }_3}\dots {\Delta }_{{\mu }_m}{P}_{{\mu }_{m+1}}\dots P_{{\mu }_n}]\}.$$

$${\mathcal{O}}_{l,m,n}^{(E)}=Z_{l,m}{O}_{l,m,n}^{\text{latt}},$$

$${\overline{\mathcal{O}}}_{l,m,n}^{(E)}=Z_{l,m}^{gg}{\overline{\mathcal{O}}}_{l,m,n}^{\text{latt}}+Z_{l,m}^{gq}{\overline{\mathcal{Q}}}_{l,m,n}^{\text{latt}}.$$

$$⟨0|{\eta }_j(\vec{p})|\vec{p},\lambda ⟩=Z_{\varphi }(\vec{p}){ϵ}_j(\vec{p},\lambda ),$$

$$C_{jk}^{2\text{pt}}(\vec{p},t)$$

$$=|Z_{\varphi }(\vec{p}){|}^2\left(e^{-\mathcal{E}t}+e^{-\mathcal{E}(T-t)}\right)\sum _{\lambda }{ϵ}_k(\vec{p},\lambda ){ϵ}_j^{\ast }(\vec{p},\lambda ).$$

$$C_{jk}^{3\text{pt}}$$

$$-⟨{\eta }_k(\vec{p},t){\eta }_j^{†}(\vec{p}{}^{\prime },0)⟩⟨\mathcal{O}(\vec{p}{}^{\prime }-\vec{p},\tau )⟩$$

$$=Z_{\varphi }^{†}(\vec{p})Z_{\varphi }(\vec{p}{}^{\prime })e^{-\mathcal{E}t}\sum _{\lambda {\lambda }^{\prime }}{ϵ}_k(\vec{p},\lambda ){ϵ}_j^{\ast }(\vec{p}{}^{\prime },{\lambda }^{\prime })⟨\vec{p},\lambda |\mathcal{O}|\vec{p}{}^{\prime },{\lambda }^{\prime }⟩$$

$$R_{jk}(\vec{p},\vec{p}{}^{\prime },t,\tau ,\mathcal{O})=$$

$$\frac{C_{jk}^{3\text{pt}}(\vec{p},\vec{p}{}^{\prime },t,\tau ,\mathcal{O})}{C_{kk}^{2\text{pt}}(\vec{p}{}^{\prime },t)}$$

$${\chi }_A^2(b,f_{i\in S,j})=\sum _{\overline{R}}\frac{{\left(M_{\overline{R}}(f_{i\in S,j}){|}_{(f_{i\notin S,j}=0)}-\overline{R}(b)\right)}^2}{\Delta {\overline{R}}^2+\hat{\Delta }(\overline{R}{)}^2}.$$

$$\hat{\Delta}(\overline{R})=\sqrt{\sum_{k\notin S}\left(M_{\overline{R}}(f_{i,j})\big{|}_{(f_{k,j}=10)}-M_{\overline{R}}(f_{i,j})\big{|}_{(f_{k,j}=0)}\right)^{2}}.$$

$${\chi }_B^2(b,f_{i,j},{\mu }_i,a_i)=\sum _{i\in S,j}\frac{{\left(\frac{{\mu }_i}{(1+a_i{\Delta }^2{)}^2}-f_{i,j}\right)}^2}{{\left(\frac{1}{2}\frac{{\partial }^2{\chi }_A^2}{\partial f_{i,j}^2}\right)}^2},$$

Full text

††preprint:

MIT-CTP/4894

I Introduction

Understanding the quark and gluon structure of hadrons and nuclei is a fundamental and compelling goal of nuclear physics. Over the last 60 years there has been extraordinary progress in both experimental measurement and theoretical understanding of the distributions of quarks inside hadrons and nuclei. There are now precise pictures of both the electromagnetic form factors and the quark distributions of nucleons and nuclei. Determining the gluonic structure of these objects to a similar level, however, requires a new generation of experiments with higher luminosity and better detectors. In particular, an electron-ion collider (EIC) designed to fulfil these needs is currently in the planning phase Accardi et al. (2012). Such a machine will provide access to a host of information about the gluonic structure of hadrons and nuclei, including transverse-momentum dependent distributions (TMDs) and gluon generalised parton distributions (GPDs). While many aspects of gluonic structure can be investigated, generalised transversity gluon distributions are of particular interest since they are purely gluonic; they do not mix with quark distributions at leading twist. In addition, the forward limit of these quantities are the double helicity flip parton distributions introduced by Jaffe and Manohar Jaffe and Manohar (1989), which provide a clean signature of non-nucleonic degrees of freedom in nuclei of spin $\geq$ 1. Away from the forward limit, these distributions exist for targets of any spin Hoodbhoy and Ji (1998); Belitsky and Mueller (2000).

In this work the generalised gluon distributions (gluon GPDs) of the spin-1 $\phi$ meson are investigated. In particular, the generalised form factors (GFFs) corresponding to the first Mellin moments of the unpolarised and transversity gluon GPDs are determined for the first time using lattice QCD (LQCD) calculations, albeit at unphysical quark masses. The forward limits of these form factors correspond to the gluon momentum fraction and a transverse momentum asymmetry, respectively. Since the transversity gluon GPD is non-zero in the forward limit only in targets of spin$\geq$1, the $\phi$ meson, which is the simplest spin-1 system, is chosen for this exploratory study.

Because of the large number of both transversity and unpolarised GFFs that contribute to the first Mellin moments of the GPDs away from the forward limit, only a subset can be cleanly determined within the technical limitations of the LQCD calculation presented here. Nevertheless, three of the seven unpolarised gluon GFFs, and one of the eight transversity GFFs, are extracted for six non-zero momentum transfers in the range $0<|\Delta^{2}|<1.2$ GeV2. The unpolarised gluon GFFs are compared with the analogous quark GFFs, to which they have a one-to-one correspondence. While there are clear quantitative differences between the unpolarised gluon and quark distributions, and also between these and the gluon transversity distribution, interpreting these differences is a challenging problem. Resolving a full three-dimensional picture of the gluon structure of the $\phi$ meson will require more precise calculations that extend to more than the lowest moment of the GPDs.

This work represents the first probe of detailed aspects of the gluonic structure of a hadron using LQCD. A number of technical aspects of the study are novel, and the improvements presented here set the stage for future studies that will map out the full gluonic structure of the $\phi$ meson and other hadrons. Applied to nucleons and nuclei, these techniques will set QCD benchmarks for an EIC.

II Gluon GFFs for spin-1 particles

GPDs encode the three-dimensional quark and gluon structure of hadrons and nuclei. They encompass the information carried by the parton distribution functions and the elastic electromagnetic form factors, describing the distribution of partons both in the transverse plane and in the longitudinal direction Diehl (2003). Through the operator product expansion, the towers of Bjorken-$x$ (Mellin) moments of the GPDs are related to matrix elements of towers of local twist-two operators. These matrix elements, in turn, are parametrised in terms of the GFFs which are the focus of this work.

There are three towers of moments of twist-2 gluon GPDs, encoding the spin-independent, spin-dependent and transversity distributions. These moments are related to matrix elements of the operators

[TABLE]

respectively, where the gluon field strength tensor is $G_{\mu\nu}$, the dual field-strength tensor is $\tilde{G}_{\mu\nu}=\frac{1}{2}\epsilon_{\mu\nu\alpha\beta}G^{\alpha\beta}$, and $\overleftrightarrow{D}=\frac{1}{2}\left(\overrightarrow{D}-\overleftarrow{D}\right)$. ‘S’ denotes symmetrisation and trace-subtraction in all free indices for Eqs. (1) and (2), and symmetrisation in the $\mu_{i}$ and and trace-subtraction in all indices for Eq. (3). The matrix elements of these operators in spin-1 states, at lowest $n$, are the focus of this work.

The off-forward matrix elements of the twist-2 operators defined above are described by GFFs. The spin-dependent gluon GFFs, which vanish at lowest-$n$ through operator symmetries, are not considered numerically in this work but are enumerated in Appendix B. With the polarisation vectors of massive spin-1 particles defined in Minkowski space as

[TABLE]

where $\lambda=\{+,-,0\}$, $m$ and $\mathcal{E}=\sqrt{|\vec{p}|^{2}+m^{2}}$ are the rest mass and energy of the state, and

[TABLE]

the spin-independent gluon GFFs are defined111This choice of basis is slightly different from that in Ref. Taneja et al. (2012), where the decomposition also includes a trace term. Note also that not all GFFs are independent for all $n$. Taneja et al. (2012) through

[TABLE]

Here, $P=(p+p^{\prime})/2$ is the average momentum and the momentum transfer is defined as $\Delta=p^{\prime}-p$. ‘S’ denotes symmetrisation and trace-subtraction in all free indices. Of these GFFs, only $B^{(n)}_{1,0}(\Delta^{2})$ and $B^{(n)}_{2,0}(\Delta^{2})$ contribute to forward-limit matrix elements. The renormalisation scheme and scale-dependence of the GFFs is suppressed here.

The transversity GFFs are defined through

[TABLE]

where the polarisation vectors $E$ and the momenta $P$ and $\Delta$ are as defined above. Here, ‘S’ denotes symmetrisation in the indices $\mu_{i}$ (the pairs {$\mu$,$\mu_{1}$} and {$\nu$,$\mu_{2}$} are antisymmetric), symmetrisation of $\mu$ and $\nu$, and trace-subtraction in all free indices. The construction of this decomposition and that of Eq. (7) follows from applying discrete symmetries and demanding the correct Lorentz structure. Only $A^{(n)}_{1,0}(\Delta^{2})$ contributes to forward-limit gluon transversity matrix elements.

III Lattice QCD calculation

In this work, a single ensemble of isotropic gauge-field configurations is used to determine the matrix elements discussed above at lowest $n$. The solutions of the systems of equations generated by various choices of polarisations and momenta in Eqs. (7) and (8) allow subsets of the GFFs to be extracted, as will be discussed in detail below. Simulations are performed with $N_{f}=2+1$ flavours of dynamical quarks, with quark masses chosen such that222Throughout this work, the $\phi$ meson is assumed to have a flavour content that is purely $\overline{s}s$ and annihilation contributions are ignored in two and three-point correlation functions. Such terms are suppressed by the Zweig rule. $m_{\pi}\sim 450(5)$ MeV and $m_{\phi}\sim 1040(3)$ MeV. A clover-improved quark action Sheikholeslami and Wohlert (1985) and Lüscher-Weisz gauge action Lüscher and Weisz (1985) are used, with the clover coefficient set equal to its tree-level tadpole-improved value. The lattices have dimensions $L^{3}\times T=24^{3}\times 64$, with lattice spacing $a=0.1167(16)$ fm Meinel . Details of this ensemble are given in Table 1 Orginos et al. (2015).

III.1 Lattice operator construction

The lowest-$n$ operators of the towers given in Eqs. (1) and (3) are considered here. Symmetrised and trace-subtracted, the Minkowski-space gluonic transversity operator for $n=2$ (Eq. (3)) does not mix with quark-bilinear operators of the same or lower dimension under renormalisation. The spin-independent gluonic operator with $n=0$ (Eq. (1)), however, mixes with the flavour singlet quark operator $\sum_{f=\{u,d,s\}}S\left[\overline{\psi}_{f}\gamma_{\mu}\overleftrightarrow{D_{\nu}}\psi_{f}\right]$, as discussed in more detail below. Moreover, the discrete symmetries of a hypercubic lattice reduce the Lorentz group to the hypercubic group H(4), creating the possibility of further mixing. Lattice operators with the appropriate continuum behaviour that do not have additional mixing with lower or same-dimensional operators were constructed, for the cases considered here, in Refs. Göckeler et al. (1996) and Detmold and Shanahan (2016).

For the gluon transversity operator in Eq. (3), operators in two irreducible representations of H(4) that do not mix with operators of same or lower dimension are investigated. Lattice operators which define bases of these representations are given explicitly in Appendix A along with their Minkowski-space analogues. These operators are constructed using the clover definition of the gluon field strength tensor, with gradient flow Lüscher (2010) applied to the links in the lattice gluon operators. The results shown use operators flowed to a total time of 1 in lattice units using a step size of 0.01.

As discussed in Ref. Detmold and Shanahan (2016), the transversity lattice operators are related to continuum Euclidean-space operators through a finite multiplicative renormalisation factor:

[TABLE]

where the subscript $(l,m,n)$ denotes the $n$th vector from the $m$th representation of a lattice operator (where such operators are identified by the subscript $l$), and $Z_{l,m}=1+\mathcal{O}(\alpha_{s})$. In this work the renormalisation factors are not computed, but it is expected based on studies of similar gluonic operators Alexandrou et al. (2016) that they are $\mathcal{O}(1)$.

Lattice operators from two representations of H(4) are considered for the spin-independent gluon operator defined in Eq. (1), with explicit definitions given in Appendix A. As noted above, these operators mix with the quark operator $\sum_{f}S\left[\overline{\psi}_{f}\gamma_{\mu}\overleftrightarrow{D_{\nu}}\psi_{f}\right]$. With the lattice operators corresponding to this quark bilinear denoted by $\overline{\mathcal{Q}}^{\text{latt}}_{l,m,n}$ with subscripts defined as above (where operators transforming irreducibly under H(4) are constructed in the same way as those for the corresponding gluonic operator), this mixing under renormalisation can be expressed as

[TABLE]

In Ref. Alexandrou et al. (2016) it is shown numerically that this mixing, i.e., the magnitude of $Z^{qg}_{l,m}$, is at the few-percent level for a similar action to the one used here, and that the renormalisation $Z^{gg}_{l,m}$ is approximately unity, with several levels of stout smearing used on the operator. These small mixing effects are neglected in the present calculation.

III.2 Determination of matrix elements

Matrix elements of the operators discussed in the previous section in the $\phi$ meson can be extracted from ratios of two and three-point correlation functions. With $\eta_{j}(\vec{p},t)$ denoting the vector of $\phi$ interpolating operators and where $\epsilon_{j}$ are Euclidean polarisation vectors related to the Minkowski expression in Eq. (4) by $\epsilon_{j}(\vec{p},\lambda)=E^{j}(\vec{p},\lambda)$, such that333Note that cubic symmetry guarantees the polarisation-independence of $Z_{\phi}(\vec{p})$.

[TABLE]

where $\vec{p}$ is the momentum of a state and $\lambda$ labels its polarisation, the two-point function can be expressed as

[TABLE]

Contributions from excited states (which are exponentially suppressed) are omitted from this expression. In analysis, care is taken to restrict to time ranges where such contamination is negligible.

Three-point correlation functions are constructed by taking the correlated product, configuration-by-configuration and source-location–by–source-location, of these two-point functions444For three-point functions with off-diagonal polarisations in the helicity basis, the two-point functions used are zero when an ensemble average is taken, but signals emerge through their correlation with the gluonic operators. The gluon transversity operators considered here are themselves zero on ensemble average. with the gluonic operators calculated as described in the previous section. While the only case for which the vacuum expectation value $\left\langle\eta_{k}(\vec{p}\,^{\prime},t)\,\eta^{\dagger}_{j}(\vec{p},0)\right\rangle\left\langle\vphantom{\eta^{\dagger}_{j}}\mathcal{O}(\vec{p}\,^{\prime}-\vec{p},\tau)\right\rangle$ is non-zero is where $\vec{p}=\vec{p}\,^{\prime}$ and $\mathcal{O}$ is a spin-independent gluon operator, a vacuum subtraction is performed for every operator and all momenta in this calculation. This correlated subtraction of zero improves the signal-to-noise ratio significantly. Inserting complete sets of states, the subtracted three-point correlators can thus be expressed as

[TABLE]

for $0\ll\tau\ll t\ll T$ (where $T$ denotes the time extent of the lattice). For the case $0\ll t\ll\tau\ll T$, $t$ is replaced by $(T-t)$ in the final line of the above expression and there is an additional multiplicative factor of $(-1)^{n_{4}}$, where $n_{4}$ is the number of temporal indices in the operator.

The two and three-point correlation functions were constructed from propagators computed using a bare quark mass $m=-0.2450$ and 5 iterations of gauge-invariant Gaussian smearing in the spatial directions at both source and sink, with interpolating operators of the form $\eta_{j}(x)=\overline{\psi}(x)\gamma_{i}\psi(x)$ in terms of smeared quark fields. On each of 1042 configurations, spaced by 10 trajectories, 96 source locations were used, and measurements were averaged over these source locations before a bootstrap analysis was performed to assess statistical uncertainties.

The leading exponential time-dependence in Eq. (13), as well as factors of $Z_{\phi}$, can be eliminated by forming the ratio:

[TABLE]

which is proportional, through factors of $m_{\phi}$ and momentum components, to the matrix elements of interest, $\langle\vec{p}\,^{\prime},\lambda^{\prime}|\mathcal{O}|\vec{p},\lambda\rangle$. For each lattice operator $\mathcal{O}$, this ratio was constructed for all diagonal and off-diagonal polarisation combinations $jk$, and for momenta up to $\vec{p}\,^{2}=4$ and $\vec{p}\,^{\prime 2}=4$, taken in all combinations that give a resultant momentum transfer up to $|\vec{\Delta}|^{2}=(\vec{p}\,^{\prime}-\vec{p}\,)^{2}=6$. Ratios for the cases with $\tau\ll t$ and $t\ll\tau$ are averaged, with the appropriate signs included as discussed following Eq. (13).

III.3 Extraction of GFFs

The Euclidean operators used here are given explicitly in Appendix A and discussed in Sec. III.1. Matrix elements of these operators, encoded in the ratios $R_{jk}(\vec{p},\vec{p}\,^{\prime},t,\tau,\mathcal{O})$ described in the previous section, are matched to GFFs by applying Eqs. (7) and (8) to the corresponding Minkowski-space operator expressions (also made explicit in Appendix A). For each basis of operators, at each value of the momentum transfer, this generates systems of equations for the GFFs, $B_{i,m}^{(2)}(\Delta^{2})$ or $A_{i,m}^{(2)}(\Delta^{2})$ depending on the operator, where each equation corresponds to one choice of operator in the basis and one set of momentum and polarisation vectors. In general, these systems are large, the relation between the GFFs and the matrix elements is not simply invertible, and the systems can not be solved for all GFFs simultaneously. This is discussed in detail below.

The extraction of the GFFs proceeds in four of steps:

Construct averages of the ratios $R_{jk}(\vec{p},\vec{p}\,^{\prime},t,\tau,\mathcal{O})$ for equivalent choices of polarisation, momenta, and operators in a given basis; 2. 2.

Fit constants to the averaged ratios in their plateau regions; 3. 3.

Determine the GFFs which are dominant in the analysis for each operator; 4. 4.

Solve the (possibly over or under-determined) system for the dominant GFFs at each $\Delta^{2}$;

each of which will be described in detail.

For each operator under consideration, the numerical values of the ratios $R_{jk}(\vec{p},\vec{p}\,^{\prime},t,\tau,\mathcal{O})$ are averaged, at the bootstrap level, over all choices of momentum and polarisation that give the same linear combination of GFFs up to a sign (by Eq. (7) or (8) as appropriate). This procedure defines a reduced set of unique, but not linearly independent, linear equations for each operator and momentum transfer.

For each averaged ratio $\overline{R}(t,\tau)$, the maximal connected plateau region in $t$–$\tau$ space is determined, where this is defined as the region where the bootstrap-level differences between all pairs of adjacent points are consistent with zero. If this maximal plateau region for an averaged ratio consists of less than 10 $(t,\tau)$ pairs, that ratio is discarded from the analysis (typical fits include many more points in the plateau region). Given this maximal plateau region, the variation in central values of fits to all sub-regions is taken as a measure of the fitting uncertainty, while the bootstrap-level fit of a constant over the maximal region gives the central value and statistical uncertainty of the fit. These uncertainties are propagated into the subsequent analysis as described below. In this analysis, all $(t,\tau)$ combinations are available, so a comprehensive elimination of excited states can be achieved; this aspect of the calculation is better controlled than for studies of quark operators where each $\tau$ (or $t$) value requires additional computation and so typically only a few values can be used. Figure 1 shows an example of such a plateau fit to an averaged ratio in the $t$–$\tau$ plane.

Because of the large number of GFFs that contribute to the off-forward matrix elements, not all can be determined from the LQCD calculations presented here. A complete extraction would require precise data from many different sets of initial and final momenta giving the same momentum transfer. This could be achieved either with new techniques allowing high-precision data to be obtained at large momenta Chambers et al. (2017); Bali et al. (2016); Syritsyn , or with very large lattices having allowed values of momentum transfer that are sufficiently closely spaced in physical units to allow binning. Given the sets of momenta available with good precision in these calculations, the linear systems generated by the matching of the LQCD results to the corresponding matrix elements in terms of GFFs do not contain enough independent equations to constrain all GFFs for some operators at some momentum transfers. In other cases, the contributions from a number of the GFFs are suppressed by several orders of magnitude relative to others, again making the extraction of these quantities impossible with the current statistical precision. Moreover, for some bases of operators, symmetries relate the coefficients of two or more GFFs for every choice of momentum and polarisation, meaning that those GFFs can not be separated by any fit, regardless of the precision of the results or the number of momenta available.

For each basis of operators, a subset of the GFFs can, however, be extracted. This set of dominant GFFs is found by inspection of the relative weights of each GFF in the system of linear equations to be solved. For example, for all gluon transversity basis operators considered, at all momentum transfers, the majority of equations in the linear systems have the coefficient of $A_{1,0}^{(2)}(\Delta^{2})$ significantly larger than the coefficients of the other GFFs (see Appendix C). For this reason, $A_{1,0}^{(2)}(\Delta^{2})$ is considered to be the dominant GFF, and its extraction is the focus of this work. Similar arguments lead to three GFFs, namely $B^{(2)}_{1,0}(\Delta^{2})$, $B^{(2)}_{2,0}(\Delta^{2})$ and $B^{(2)}_{4,0}(\Delta^{2})$, being targeted in the spin-independent case. More GFFs are resolvable at some particular momenta, but the aim of this work is to obtain a subset of the GFFs that can be determined consistently at all momenta, as described below.

To achieve this, least-$\chi^{2}$ fits at the bootstrap level are performed to the systems of equations for each basis of operators, with fitting uncertainties assigned to each averaged ratio $\overline{R}$ as described above. Such fits are performed multiple times, fitting to every subset of the GFFs that includes those classed as dominant, with the GFFs which are not fit set to 0 and assigned an uncertainty of 10 in the $\chi^{2}$ fit555While there are no current bounds on the magnitudes of the GFFs to support this choice, an order of magnitude variation would be surprising. Moreover, as fits to all subsets of GFFs are included in the analysis, significant contributions outside of this bound from a ‘sub-dominant’ GFF would become apparent from inconsistencies between fits to different subsets. This is not observed.. There are thus 16 sets of fits for the unpolarised operator and 128 sets of fits for the transversity operator. The variation in central values over the fits to different subsets with acceptable values of $\chi^{2}$ is included in quadrature as a second fitting systematic uncertainty on the final results, while the central values and statistical and plateau fitting uncertainties are taken from the bootstrap fits over the minimal set of dominant GFFs.

Precisely, for a fit to the subset of GFFs $f_{i\in S,j}$, where $S$ denotes the set of GFFs which are fit over and the subscript $j$ labels the discrete values of the momentum transfer $\Delta^{2}$, the first contribution to the $\chi^{2}$ function can be expressed as:

[TABLE]

Here $\overline{R}(b)$ denotes the set of averaged plateau values extracted from the ratios discussed previously, for a given bootstrap, and $\Delta\overline{R}$ represents the statistical and fitting uncertainties on these quantities, determined as described earlier. The label $b$ indicates that this $\chi^{2}$ is formed for each bootstrap. $M_{\overline{R}}$ represents the expectations, from Eq. (7) or (8), for the averaged ratio $\overline{R}$, in terms of the GFFs, where those GFFs not in the subset $S$ are set to zero. The quantity $\hat{\Delta}$ assigns an uncertainty of 10 to those GFFs which are not fit to:

[TABLE]

Choosing an uncertainty of 100 gives entirely consistent results, albeit with larger uncertainties.

In the fits that are performed, an additional contribution is added to the $\chi^{2}$ function, representing a dipole function fit to each of the dominant GFFs as a whole. This addition has the effect of correlating the solutions of the systems of overdetermined equations at different values of $\Delta^{2}$, which results in GFFs that are somewhat more smoothly behaved and reduces the uncertainty on the $\left|\Delta^{2}\right|=\{2,3\}$ (in lattice units) points in particular. The particular $\Delta^{2}$ points that see improvement, fit alone, are less well constrained than others as there are less combinations of polarisations/momenta available. Precisely, this additional contribution to the total $\chi^{2}$ can be expressed as

[TABLE]

where $\mu_{i}$ and $a_{i}$ are the parameters of the dipole fits to the GFF labelled by $i$. Note that $\chi^{2}_{A}$ is a quadratic function of the GFFs $f_{i,j}$, so the denominator weights each term by how well the relevant GFF is determined. A correlation matrix is not used, since $\chi^{2}_{A}$ does not include cross-terms between different momenta.

Results are obtained by minimization of the total:

[TABLE]

to determine the dominant GFFs at each momentum transfer. The mean and standard deviation over the GFFs determined on each bootstrap are used as the central value and first uncertainty on the quoted results. In addition, the uncertainties on the GFFs determined by minimizing $\sum_{b}\chi^{2}_{\text{tot}}(b,f_{i,j},\mu_{i},a_{i})$ are included in quadrature. In most cases, where the bootstrap uncertainties accurately reflect the fitting uncertainties, this addition has little effect. In the few cases where sub-dominant form factors are significant at certain momenta (and setting them to be centered around 0 distorts the bootstrap fits), this addition inflates the fitting uncertainties considerably. Figure 2 illustrates one of the critical aspects of the fitting procedure, namely solving the linear system of constraints. The results in Figs. 3 and 4 are shown both with and without the additional smoothing constraint discussed above.

III.4 Results: Gluon GFFs

The procedure described in the previous section allows the determination of the $\Delta^{2}$-dependence of one of the eight gluon transversity GFFs, and three of the seven unpolarised gluon GFFs. At some momentum transfers, additional GFFs (or linear combinations thereof) can be determined. This study is, however, focussed on determining the $\Delta^{2}$-dependence of GFFs as a whole.

The single gluon transversity GFF that can be determined, which is also the only transversity GFF that contributes in the forward limit, is shown in Fig. 3. The forward limit of this quantity defines a transverse momentum asymmetry, which was previously determined on the same gauge ensemble as used here Detmold and Shanahan (2016). While the results labelled basis 1 are more precise than those for basis 2, as there are fewer momentum and polarisation combinations that give non-zero matrix elements of basis-2 operators, they agree within uncertainties.

Results for the three spin-independent gluon GFFs that can be determined are shown in Fig. 4. Again, the results from the two distinct bases are broadly consistent at the $1\sigma$ level, with comparable precision achieved in both data sets. While the spin-independent gluon GFFs are clearly quantitatively different from the gluon transversity GFF, this difference is as yet hard to interpret meaningfully, given the limited subset of GFFs which are determined. Within the uncertainties of this calculation, the effects of renormalisation, which have been neglected, are not significant Alexandrou et al. (2016).

III.5 Quark GFFs

To interpret the gluonic observables obtained in this study, it is interesting to compare them with the analogous quark GFFs and the electromagnetic form factors. As this first lattice determination of the gluon GFFs is performed at a single unphysical value of the quark masses, and no extrapolation to the physical point is performed, it is natural to compare gluon and quark GFFs at the same unphysical parameters. While the transversity gluon GFFs have no direct quark analogues, the spin-independent gluon GFFs have a one-to-one correspondence with the spin-independent quark GFFs. These quantities, as well as the quark electromagnetic form factors, are calculated using the same lattice setup and analysis procedures described previously, with minor differences as detailed below.

For the spin-independent strange quark GFFs and the electromagnetic form factors of the $\phi$ meson, the relevant three-point functions, $C_{jk}^{\text{3pt}}$ (Eq. (13)), are determined using the quark bilinear operators

[TABLE]

respectively. The calculations omit the disconnected couplings of the sea quarks to the currents, and $\cal{O}(\alpha)$ mixing with the gluonic operators is ignored. In addition, a single contraction is considered, so the system under consideration should be thought of as an $\overline{s}s^{\prime}$ meson.

Matrix elements of the spin-independent quark operators have precisely the same form as those of the spin-independent gluon operators, given explicitly in Eq. (7). The GFFs $B_{i,m}^{q(n)}$ are given the additional superscript ‘$q$’ to identify them as the quark analogues. The decomposition of the electromagnetic current matrix elements, which have no gluonic analogues, into form factors for spin-1 particles is Hedditch et al. (2007); Owen et al. (2015)

[TABLE]

where, as before, $P=(p+p^{\prime})/2$ is the average momentum, and the momentum transfer is defined as $\Delta=p^{\prime}-p$. The initial and final-state energies are $\mathcal{E}$ and $\mathcal{E}^{\prime}$ (and the polarisations, as before, are $E$ and $E^{\prime}$).

The three-point functions needed for the LQCD determination of these matrix elements are constructed using sequential propagators with fixed sink momentum $\vec{p}\,^{\prime}=0$, since each additional momentum used carries additional computational cost. As a result, there is less information (in the form of fewer independent equations in the system that determines the GFFs or FFs) than in the gluonic study discussed above. Moreover, unlike in the gluon case, the construction of three-point functions with different sink times requires separate sequential propagator computation. Three sink times, $t_{\text{sink}}\in\{12,14,16\}$ in lattice units, are used here. The statistical behaviour of these matrix elements is, however, less noisy than that of the gluon operator correlators.

As was discussed previously for the gluon case, calculations are performed using lattice operators which transform irreducibly under the hypercubic group H(4). The vector current is implemented using the naive local current, related to the continuum current by a renormalisation factor $Z_{V}$ that is determined by demanding unit charge. This transforms in the $\tau_{1}^{(4)}$ representation. For the spin-independent quark operator in Eq. (19), the irreducible representations have exactly the same form as those for the spin-independent gluon operator, which are given explicitly in Appendix A. A single lattice representation of the spin-independent quark operator is used, namely $\overline{Q}^{E}_{2,14}$, analogous to the gluon quantity labelled as $\overline{O}^{E}_{2,14}$ in Appendix A.

Given the three point functions constructed as described above, the analysis to extract the spin-independent quark GFFs proceeds in the same way as was detailed for the spin-independent gluon GFFs above. The only modification is that, with only three sink times available, identified plateaus in the $(t$–$\tau)$ plane are considered acceptable if they include a minimum of 5 timeslices. For the electromagnetic form factors, where a complete extraction of all FFs is possible, both this method (with all three FFs considered ‘dominant’) and a separate analysis, in which ratios of two and three-point functions are combined to give direct extractions of the individual FFs, are used. This second method is discussed in detail in Appendix D and gives results consistent with the more general method that must be used in the cases where larger numbers of form factors contribute. This comparison also gives confidence that this latter method is reliable.

This procedure allows the determination of the $\Delta^{2}$-dependence of three of the seven unpolarised quark GFFs. These three quantities are the direct analogues of the three spin-independent gluon GFFs that were determined as described in the previous sections, and can therefore be compared with these one-to-one. The results are shown in Fig. 5. The forward-limit quantities $B_{1,0}^{q(2)}(0)$ and $B_{2,0}^{q(2)}(0)$ satisfy the Soffer-type bounds derived on their relation in Ref. Bacchetta and Mulders (2001). These quantities were previous calculated for a heavy $\rho$ meson in a quenched calculation Best et al. (1997). Translating the results presented here to the notation in that work, the linear combination of $B_{1,0}^{q(2)}(0)$ and $B_{2,0}^{q(2)}(0)$ named $a_{1}$ has a comparable value, but the combination $d_{1}$ has a different sign and magnitude, possibly due to the effects of unquenching. The results for the form factors of the vector current are shown in Fig. 6 (and with a different choice of decomposition in Fig. 8). As discussed above, in this case all three form factors can be determined.

IV Discussion

While many observables related to the quark structure of hadrons and nuclei have now been both measured experimentally and calculated at some level from QCD, including electromagnetic form factors and quark momentum distributions, the gluon structure of these particles remains far more mysterious. One question of fundamental interest is the spatial distributions of gluons relative to that of quarks; is the ‘gluonic radius’ of a hadron larger, smaller, or of a similar size to the corresponding quark radius. This, while interesting, is a somewhat nebulous question. The observables considered in this work, namely the spin-independent and transversity gluon GFFs, each define a gluonic radius, each of which could, in principle, be different, and there is no unique basis for the GFF decomposition. The spin-independent gluon GFFs have direct quark analogues and can be compared with these on a one-to-one basis. While there are clear quantitative differences between each of the the quark and gluon GFFs in the three matched pairs that were studied here, these differences vary (i.e., the spin-independent quark and gluon GFFs are not universally related by an approximate scaling or sign change, as shown in Fig. 7). The gluon transversity distribution is, again, quantitatively different, but this difference is hard to interpret in a physical sense. Resolving a full three-dimensional picture of the gluon structure of the $\phi$ meson will require more precise calculations that are able to resolve the entire basis of GFFs for the first moments of the gluon distributions, and also extend to higher moments.

An EIC will, for the first time, allow experimental measurements of gluon GFFs in nucleons and nuclei. The present work represents a demonstration that QCD predictions of gluonic structure quantities can be obtained using LQCD. Future LQCD studies of nucleons and nuclei with fully-controlled uncertainties will inform the design and targets of an EIC experimental program, guide the interpretation of first measurements, and, for some quantities, act as theory benchmarks for the EIC.

Acknowledgements

This work is supported by the U.S. Department of Energy under Early Career Research Award DE-SC0010495 and under Grant No. DE-SC0011090. The work of WD is supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, within the framework of the TMD Topical Collaboration. The Chroma software library Edwards and Joo (2005) was used in the data analysis. We thank Ian Cloët, Bob Jaffe, James Maxwell, Richard Milner, and Ross Young for helpful discussions, and Kostas Orginos for the production of the gauge configurations used in this work. We thank W. Cosyn, A. Freese and B. Pire for bringing errors in a previous version of this work to our attention.

Appendix A Explicit expressions for lattice operators

Bases of Euclidean operators that transform irreducibly under the hypercubic group H(4) are defined, for different symmetry classes of operators, in Ref. Göckeler et al. (1996). The same notation as in that work is used here, with irreducible representations labelled by $\tau_{i}^{(n)}$ where $n$ denotes the dimension of the representation and $i$ enumerates representations.

The Euclidean analogue of the transversity operator defined in Eq. (3), for $n=2$, is built from the quantity

[TABLE]

where the clover definition of the gluon field strength tensor is used for $G_{\mu\nu}^{(E)}$ in the numerical calculations. While there are three irreducible representations that are safe from mixing with lower or same-dimensional quark bilinear operators, namely $\tau_{1}^{(2)}$, $\tau_{2}^{(6)}$ and $\tau_{2}^{(2)}$ Detmold and Shanahan (2016), only a subset of the basis operators of these representations, which provide the cleanest signals, are considered here. In particular, three of the six basis operators from $\tau_{2}^{(6)}$, and both basis operators from $\tau_{2}^{(2)}$, have vanishing matrix elements in states that have small momenta in lattice units (for example, with zero momentum transfer, these operators only have nonvanishing matrix elements for boosts with $\vec{p}\,^{2}\geq 3$ Detmold and Shanahan (2016)). The five remaining operators from representations $\tau_{1}^{(2)}$ and $\tau_{2}^{(6)}$ are non-vanishing in a larger number of external states, and so provide the most information. These are considered in this work.

For $\tau_{1}^{(2)}$, the basis vectors are Göckeler et al. (1996); Detmold and Shanahan (2016):

[TABLE]

The $\tau_{2}^{(6)}$ vectors which are used here are:

[TABLE]

The Euclidean analogue of the spin-independent gluonic operator in Eq. (1), for $n=0$, is constructed from

[TABLE]

Basis operators from two irreducible representations are considered here. A basis of operators for the $\tau_{1}^{(3)}$ representation is:

[TABLE]

For $\tau_{3}^{(6)}$ the vectors are:

[TABLE]

where only the operators with $\nu=4$ are used here, as they provide the cleanest signals at most momenta.

The Minkowski-space analogue of each basis operator is determined by applying the relations

[TABLE]

to the Euclidean-space form. These Minkowski operators are used, as described in Sec. III.2, to match the numerical LQCD results for the matrix elements of operators to the expressions for these quantities in terms of GFFs.

Explicit expressions for the Euclidean and Minkowski-space quark operators are given in Ref. Gockeler et al. (1996). The structure of H(4) irreducible representations constructed from the Euclidean operators is identical to that given in Eqs. (29) through (32) for the gluon case.

Appendix B GFFs of the spin-dependent gluon operator in spin-1 states

Although the spin-dependent gluon distributions are not studied numerically in this work, the GFF decomposition of matrix elements of the spin-dependent gluon operator in Eq. (2) is derived here. It can be expressed as (equivalent to the decomposition in Ref. Cosyn et al. (2018) after application of the Schouten identities)

[TABLE]

Here $S$ denotes symmetrisation and trace-subtraction in all free indices. The average momentum is defined as $P=(p+p^{\prime})/2$, and the momentum transfer is $\Delta=p^{\prime}-p$. Note that, because of the symmetries of the operator, $\tilde{B}_{i,0}^{(2)}(\Delta^{2})=0$ for all $i\in\{1\ldots 7\}$. Furthermore, none of the GFFs $\tilde{B}_{i,m}^{(n)}(\Delta^{2})$ contribute to forward-limit matrix elements, which must vanish.

Appendix C Example of a system of linear equations for the GFFs

This section gives an explicit example of the systems of equations that are solved to extract the GFFs as described in Section III.3. The case shown is for the gluon transversity operator in the first basis ($\tau_{1}^{(2)}$), at the first non-zero momentum transfer. This is the simplest case at non-zero momentum, as this basis includes only two operators.

Matching plateau-fits to the ratios $R_{jk}(\vec{p},\vec{p}\,^{\prime},t,\tau,\mathcal{O})$ (Eq. (14)) to linear combinations of the GFFs $A_{i,0}^{(2)}(\Delta^{2}=1)$ as described in Section III.3, for all choices of polarisation and momenta (up to $\vec{p}\,^{2}=4$ and $\vec{p}\,^{\prime 2}=4$) that give $\Delta^{2}=1$ in lattice units and for all transversity operators in basis 1 ($\tau_{1}^{(2)}$), defines 154 linear equations. This system is reduced to 39 equations by averaging the ratios (before fitting in $t$ and $\tau$) over all choices of momentum and polarisation that give the same linear combination of GFFs up to a sign. The polarisations and momenta defining one member of each reduced set are given in Table 2.

The system of equations determining $A_{i,0}^{(2)}(\Delta^{2}=1)$ for this basis, at this momentum, can be expressed as:

$\begin{pmatrix}0.604&0.0424&0&0&0&0&0.0588&0\\ 0.592&-2.45\times 10^{-3}&0.0785&-0.0785&6.58\times 10^{-3}&-0.0992&-0.103&-4.15\times 10^{-3}\\ 0.485&0.0429&0&0&0&0&0.0379&0\\ 0.481&0.0431&-3.02\times 10^{-5}&3.02\times 10^{-5}&-2.53\times 10^{-6}&-4.03\times 10^{-7}&0.0374&-1.69\times 10^{-8}\\ 0.475&-3.29\times 10^{-3}&0.0791&-0.0791&6.59\times 10^{-3}&-0.0791&-0.0824&-3.29\times 10^{-3}\\ 0.353&-7.97\times 10^{-4}&0.0385&-0.0385&3.28\times 10^{-3}&-0.0598&-0.0631&-2.54\times 10^{-3}\\ 0.347&-0.0382&0&0&0&0&0.0962&0\\ 0.258&0.0806&0&0&0&0&-0.0374&0\\ 0.258&0.0808&0&0&0&0&-0.0379&0\\ 0.253&0.101&-8.60\times 10^{-4}&8.60\times 10^{-4}&-7.20\times 10^{-5}&6.32\times 10^{-7}&-0.0588&2.65\times 10^{-8}\\ 0.239&-1.66\times 10^{-3}&0.0401&-0.0401&3.29\times 10^{-3}&-0.0393&-0.0402&-1.61\times 10^{-3}\\ 0.238&-1.65\times 10^{-3}&0.0396&-0.0396&3.29\times 10^{-3}&-0.0396&-0.0412&-1.65\times 10^{-3}\\ 0.228&-0.0581&8.30\times 10^{-4}&-8.30\times 10^{-4}&6.94\times 10^{-5}&-1.04\times 10^{-6}&0.0962&-4.33\times 10^{-8}\\ 0.228&-0.0379&0&0&0&0&0.0758&0\\ 0.0590&-0.0109&0.139&-0.139&0.0112&-4.97\times 10^{-3}&-3.94\times 10^{-4}&-8.24\times 10^{-6}\\ 0.0578&-2.56\times 10^{-4}&9.42\times 10^{-3}&-9.42\times 10^{-3}&3.89\times 10^{-4}&-4.65\times 10^{-3}&2.51\times 10^{-4}&5.25\times 10^{-6}\\ 0.0338&1.59\times 10^{-3}&-0.128&0.128&-0.0107&3.18\times 10^{-4}&0.0154&1.33\times 10^{-5}\\ 0.0183&6.36\times 10^{-3}&-1.29\times 10^{-4}&1.29\times 10^{-4}&3.84\times 10^{-4}&4.84\times 10^{-3}&5.99\times 10^{-3}&5.18\times 10^{-6}\\ 0.0155&-4.78\times 10^{-3}&-0.128&0.128&-0.0111&-4.52\times 10^{-3}&9.41\times 10^{-3}&8.14\times 10^{-6}\\ 1.19\times 10^{-3}&-0.0106&0.129&-0.129&0.0108&-3.22\times 10^{-4}&-6.45\times 10^{-4}&-1.35\times 10^{-5}\\ 0.549&2.44\times 10^{-3}&0&0&0&0&0.0895&0\\ 0.546&-1.88\times 10^{-3}&0.0676&-0.0676&5.69\times 10^{-3}&-0.0918&-0.0960&-3.86\times 10^{-3}\\ 0.498&0.0710&0&0&0&0&0.0123&0\\ 0.480&-2.37\times 10^{-3}&0.0685&-0.0685&5.70\times 10^{-3}&-0.0799&-0.0828&-3.33\times 10^{-3}\\ 0.429&0.0714&0&0&0&0&0&0\\ 0.424&0.0834&-5.14\times 10^{-4}&5.14\times 10^{-4}&-4.30\times 10^{-5}&1.33\times 10^{-7}&-0.0123&5.55\times 10^{-9}\\ 0.412&2.85\times 10^{-3}&0&0&0&0&0.0657&0\\ 0.412&-2.85\times 10^{-3}&0.0685&-0.0685&5.70\times 10^{-3}&-0.0685&-0.0714&-2.85\times 10^{-3}\\ 0.409&-8.65\times 10^{-3}&4.61\times 10^{-4}&-4.61\times 10^{-4}&3.86\times 10^{-5}&-8.30\times 10^{-7}&0.0771&-3.47\times 10^{-8}\\ 0.0674&-6.43\times 10^{-3}&0.0856&-0.0856&6.70\times 10^{-3}&-5.55\times 10^{-3}&-8.26\times 10^{-5}&-1.73\times 10^{-6}\\ 0.0656&4.96\times 10^{-4}&-9.21\times 10^{-4}&9.21\times 10^{-4}&-6.37\times 10^{-6}&-0.0119&-0.0132&-5.32\times 10^{-4}\\ 0.0514&-0.0685&0&0&0&0&0.0771&0\\ 0.0347&-0.0124&0.155&-0.155&0.0127&-3.05\times 10^{-3}&-6.00\times 10^{-4}&-1.26\times 10^{-5}\\ 0.0327&5.99\times 10^{-3}&-0.0692&0.0692&-6.03\times 10^{-3}&-2.50\times 10^{-3}&5.17\times 10^{-4}&1.08\times 10^{-5}\\ 0.0301&4.59\times 10^{-3}&-0.0738&0.0738&-5.95\times 10^{-3}&2.98\times 10^{-3}&0.0123&1.07\times 10^{-5}\\ 0.0285&-1.84\times 10^{-3}&-0.147&0.147&-0.0126&-2.43\times 10^{-3}&0.0143&1.24\times 10^{-5}\\ 0.0171&0.0685&0&0&0&0&-0.0657&0\\ 0.0146&0.0920&-9.75\times 10^{-4}&9.75\times 10^{-4}&-8.17\times 10^{-5}&9.63\times 10^{-7}&-0.0895&4.03\times 10^{-8}\\ 1.59\times 10^{-3}&6.43\times 10^{-3}&0.0736&-0.0736&6.61\times 10^{-3}&5.40\times 10^{-3}&-1.97\times 10^{-3}&-1.71\times 10^{-6}\end{pmatrix}$ $\begin{pmatrix}A^{(2)}_{1,0}(1)\\[2.0pt] A^{(2)}_{2,0}(1)\\[2.0pt] A^{(2)}_{3,0}(1)\\[2.0pt] A^{(2)}_{4,0}(1)\\[2.0pt] A^{(2)}_{5,0}(1)\\[2.0pt] A^{(2)}_{6,0}(1)\\[2.0pt] A^{(2)}_{7,0}(1)\\[2.0pt] A^{(2)}_{8,0}(1)\end{pmatrix}$ $=\begin{pmatrix}0.179(36)\\ 0.150(38)\\ 0.152(30)\\ 0.154(37)\\ 0.129(32)\\ 0.056(31)\\ 0.067(41)\\ 0.056(35)\\ 0.069(21)\\ 0.093(36)\\ 0.028(32)\\ 0.041(27)\\ 0.012(33)\\ 0.029(30)\\ 0.024(11)\\ -0.005(21)\\ -0.0056(96)\\ -0.002(11)\\ 0.009(16)\\ 0.0162(91)\\ 0.086(26)\\ 0.131(31)\\ 0.155(33)\\ 0.086(33)\\ 0.098(16)\\ 0.094(17)\\ 0.088(27)\\ 0.114(25)\\ 0.075(27)\\ 0.034(25)\\ -0.006(22)\\ -0.001(31)\\ 0.022(11)\\ 0.014(16)\\ 0.0010(16)\\ 0.0008(85)\\ 0.018(23)\\ 0.001(29)\\ 0.005(18)\\ \end{pmatrix}$

[TABLE]

where the numbers and uncertainties on the right hand side of the equation come from the plateau fits to averaged ratios obtained as described in the main text. The ordering of the rows is as in Table 2.

Appendix D Direct solution of form factor decomposition for electromagnetic current

Since only three form factors contribute to matrix elements of the electromagnetic current, a direct solution of the constraint equations relating ratios of three-point and two-point functions to the form factors is straightforward Hedditch et al. (2007); Owen et al. (2015). This extraction is performed as a check on the more general method discussed in Section III.

For each momentum transfer, ${\Delta^{2}}$, three ratios of two-point and three-point functions are required to extract the form factors at that momentum. At zero momentum transfer, only the $G_{1}$ form factor can be determined. In terms of the ratios

[TABLE]

for the currents $J^{i}=\bar{\psi}\gamma^{i}\psi$, where $R_{jk}(\vec{p},\vec{p}\,^{\prime},t,\tau,\mathcal{O})$ is defined in Eq. (14) and dependence on the current and sink times is suppressed, the generic form of the solution for the FFs can be expressed as

[TABLE]

Here $X=C,M,Q$ labels the Sachs form factors, which are related to the basis used in Eq. (21) by

[TABLE]

One choice of the combinations $N_{X,f}$ for each momentum used, given that only zero sink momentum sequential propagators were computed, is given in Table 3.