aainstitutetext: Physics Department, William & Mary, Williamsburg, VA 23187,
USAbbinstitutetext: Thomas Jefferson National Accelerator Facility, Newport
News, VA 23606, USAccinstitutetext: Department of Physics, Old Dominion
University, Norfolk, Virginia, USAddinstitutetext: Aix Marseille Univ,
Université de Toulon, CNRS, CPT, Marseille, France.

#  Towards Unpolarized GPDs from Pseudo-Distributions

Hervé Dutrieux b Robert G. Edwards b Colin Egerer b Joseph Karpie a
Christopher Monahan a,b Kostas Orginos c,b Anatoly Radyushkin b David Richards
b Eloy Romero d Savvas Zafeiropoulos
(on behalf of the HadStruc Collaboration<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>mrs.fr

###### Abstract

We present an exploration of the unpolarized isovector proton generalized
parton distributions (GPDs) $H^{u-d}(x,\xi,t)$ and $E^{u-d}(x,\xi,t)$ in the
pseudo-distribution formalism using distillation. Taking advantage of the
large kinematic coverage made possible by this approach, we present results on
the moments of GPDs up to the order $x^{3}$ – including their skewness
dependence – at a pion mass $m_{\pi}=358$ MeV and a lattice spacing $a=0.094$
fm.

††preprint: JLAB-THY-24-4059

##  1 Introduction

Generalized parton distributions (GPDs) Müller  _et al._ (1994); Ji (1997a);
Radyushkin (1996a, 1997) characterize the three-dimensional internal structure
of hadrons in terms of the parton’s longitudinal momentum fraction $x$, the
fraction of longitudinal momentum transfer between the incoming and outgoing
hadrons, or skewness $\xi$, and the invariant momentum transfer $t$. The
three-dimensional nature of the distributions, in contrast to the one-
dimensional parton distribution functions (PDFs), enables insights into
crucial features such as the orbital angular momentum carried by the quarks
and gluons Ji (1997b, a). Generalizing the elastic form factors, GPDs can also
be used to define radial profiles of energy or pressure distribution of the
partonic matter Polyakov (2003); Polyakov and Schweitzer (2018).

The principal experimental means of gaining insight into GPDs is through
exclusive processes, notably deeply virtual Compton scattering (DVCS) Ji
(1997a); Radyushkin (1996a) and deeply virtual meson production (DVMP)
Radyushkin (1996b); Collins  _et al._ (1997). There is a worldwide
experimental program dedicated to these processes, including experiments at
Jefferson Lab at both 6 and 12 GeV Dudek  _et al._ (2012) and the upcoming
Electron-Ion Collider (EIC) Accardi  _et al._ (2016); Abdul Khalek  _et al._
(2022). However, these experiments only give indirect access to GPDs. In the
framework of collinear factorization, the amplitude of these processes depends
on GPDs through a convolution with a perturbative kernel, which gives a
practical access to some specific features of the GPDs (notably the diagonal
region $x\approx\xi$ or some specific integrals of the GPDs). Therefore the
insights they provide are primarily two-dimensional, or “moment-like” Qiu and
Yu (2023a, b), the three-dimensional behavior being obscured. A manifestation
of this inverse problem is encapsulated in “shadow GPDs” Bertone  _et al._
(2021), wherein a range of GPDs can give rise to extremely similar DVCS cross-
sections. In the limit of small Bjorken-$x$, corresponding to small $x$ and
$\xi$, the enhanced effects of perturbative evolution may give a hope of a
better control of the uncertainty propagation in this inverse problem Freund
(2000); Dutrieux  _et al._ (2023); Moffat  _et al._ (2023). However, in the
regime of moderate to large $x$, which precisely corresponds to the
kinematical region accessible to lattice QCD, a model independent extraction
of GPDs from those experimental processes seems difficult to achieve. It is
especially the case in the so-called ERBL region ($|x|<|\xi|$) where the lack
of theoretical constraints like positivity Radyushkin (1999); Pire  _et al._
(1999); Diehl  _et al._ (2001); Pobylitsa (2002) makes the deconvolution
problem particularly poorly behaved Dutrieux  _et al._ (2022). For example,
the very different results obtained when extracting gravitational form factors
from DVCS data using different parametrizations Burkert  _et al._ (2018);
Kumerički (2019); Dutrieux  _et al._ (2021) demonstrate a relative lack of
sensitivity of the current DVCS dataset to this quantity. In contrast, these
form factors have been probed with great precision on the lattice using local
composite operators Hagler  _et al._ (2003, 2008); Bali  _et al._ (2019);
Alexandrou  _et al._ (2020a); Hackett  _et al._ (2024).

GPDs provide therefore an ideal setup where, at least for the time-being, only
lattice QCD can provide systematically controlled information. Besides DVCS
and DVMP, an increasing number of exclusive processes with enhanced
sensitivity to GPDs are being studied, but their experimental exploration is
only at its beginning: double DVCS Guidal and Vanderhaeghen (2003); Belitsky
and Mueller (2003), where a recent study of measurability looks promising Deja
_et al._ (2023), di-photon production Pedrak  _et al._ (2020); Grocholski  _et
al._ (2022), photon meson pair production Boussarie  _et al._ (2017);
Duplančić  _et al._ (2018), or the general category of single diffractive hard
exclusive processes Qiu and Yu (2022, 2023a). The maturity of this new
phenomenology of GPDs in the perspective of the EIC therefore complements the
progress of the extraction of GPDs on the lattice. Preliminary joint lattice-
experimental phenomenology of GPDs are already being explored e.g. Guo  _et
al._ (2022, 2023); Riberdy  _et al._ (2024).

GPDs are defined as the matrix elements of operators separated along the light
cone. As in the case of the PDFs, their direct calculation is precluded on a
Euclidean lattice, except for low-order moments which can be computed from
local operators. New developments in the last decade have led to the
possibility of computing the $x$-dependence of parton distributions using non-
local operators with a space-like separation Liu and Dong (1994); Detmold and
Lin (2006); Braun and Müller (2008); Ji (2013, 2014); Radyushkin (2017);
Chambers  _et al._ (2017); Ma and Qiu (2018); Bali  _et al._ (2018). The
application of some of those approaches to GPDs was performed in Ji  _et al._
(2015); Liu  _et al._ (2019); Radyushkin (2019); Ma  _et al._ (2022) and
pioneering calculations have been conducted in Chen  _et al._ (2020); Lin
(2021); Alexandrou  _et al._ (2020b); Lin (2022); Alexandrou  _et al._ (2022);
Bhattacharya  _et al._ (2022, 2023a, 2024a, 2023b); Lin (2023); Bhattacharya
_et al._ (2024b); Hannaford-Gunn  _et al._ (2024).

### Main phenomenological results

We present a formalism for the computation of pseudo-GPDs using distillation.
This technique allows us to probe the GPDs on a large kinematic range in terms
of skewness and momentum transfer as presented in Fig. 1. For numerical
application, we use 186 combinations of initial and final momenta giving 116
kinematic values of $(\xi,t)$. This grants us access to Mellin moments of the
GPDs sensitive to the skewness on top of those only sensitive to the forward
limit $\xi=0$. We use a lattice ensemble with pion mass $m_{\pi}=358$ MeV,
lattice spacing $a=0.094$ fm, and lattice volume 32${}^{3}\times$64\. For the
extraction of the GPD moments, we have used hadron momenta up to 1.4 GeV,
yielding maximal Ioffe time values of $\nu\leq 0.6\,z/a$ where $z$ is the
space-like separation of the non-local operator. This corresponds to
$\nu_{\max}=3.5$ using $z=6a$. As demonstrated in the context of PDFs in
Egerer  _et al._ (2021), we can reach momenta at least twice as large when
using the distillation framework with momentum smearing, allowing for a study
of the $x$-dependence of the parton distributions. However, we reserve this
effort for a subsequent paper, as preliminary exploration at large momentum
hinted at the necessity of a better control of excited state contamination and
other lattice systematic uncertainties.

Figure 1: Kinematic coverage presented in this study, consisting of 186 pairs
$(\vec{p}_{f},\vec{p}_{i})$ and 116 points in the $(\xi,t)$ plane. The values
of the Ioffe time $\nu=z\cdot P$ are quoted for $z=a=0.094$ fm. By using up to
$z=6a$, $\nu$ reaches up to 3.5 in some kinematics.

For the reader interested in the phenomenology of GPDs beyond the framework of
lattice calculations, we summarize in Fig. 2 the main physics results of this
study, in the form of moments of the isovector GPDs $H^{u-d}$ and $E^{u-d}$:

$\displaystyle\int_{-1}^{1}\mathrm{d}x\,x^{n-1}\begin{pmatrix}H^{u-d}\\\
E^{u-d}\end{pmatrix}(x,\xi,t)=\sum_{k=0\textrm{
even}}^{n-1}\begin{pmatrix}A_{n,k}(t)\\\
B_{n,k}(t)\end{pmatrix}\xi^{k}\pm\textrm{mod}(n+1,2)\xi^{n}C_{n}(t)\,,$ (1)

where $C_{n}(t)$ are moments of the isovector $D$-term with $+$ for the $H$
GPD and $-$ for the $E$ GPD. For brevity, the results are quoted at the scale
$\mu=2$ GeV in the form of a dipole fit with the value at $t=0$ and a dipole
mass $\Lambda_{n,k}$ following the example:

$\displaystyle
A_{n,k}(t)=A_{n,k}(t=0)\left(1-\frac{t}{\Lambda_{n,k}^{2}}\right)^{-2}\,.$ (2)

Figure 2: Short summary of the extraction of GPD Mellin moments using a dipole
representation. The quoted uncertainty contains an evaluation of statistical
and excited state uncertainty. More details starting from section 5.

More sophisticated treatment of the data is presented starting from section 5
of this document, in particular using a more flexible parametrization of the
$t$-dependence. The quoted uncertainty is constructed by merging the 68%
confidence intervals of two different extractions to include an evaluation of
excited state contamination. This systematic uncertainty is mostly responsible
for the asymmetric errors. A slight asymmetry also arises for statistical
reasons for the most poorly determined moments (notably $A_{4,2}$ and
$B_{4,2}$) due to the non-linearity of the dipole fit.

These results are obtained at an unphysical pion mass of 358 MeV without an
attempt at extrapolating to the physical pion mass. Likewise, as we will
discuss in detail further in the document, we observe clear signs of
systematic effects, either attributable to the lack of continuum limit
(discretization effects) or the uncertainty in the light-cone limit of our
space-like matrix elements (higher-twist effects). We use data with
separations up to $z=0.6$ fm, which presents some challenges with respect to
the procedure of perturbative matching as discussed later. With these
limitations in mind, we hope that the results, and in particular the novel
characterization of the skewness-dependent moments $A_{3,2},B_{3,2}$ and
$A_{4,2},B_{4,2}$ can provide insights for theoretical models of GPDs and
extractions from the data.

The paper is organized as follows. In section 2, we review the formalism of
GPDs in the short-distance factorization, notably deriving the matching
formulas for the moments sensitive to the skewness. In section 3, we present
our calculation strategy on the lattice, stressing the extensive coverage of
the kinematic domain enabled by distillation. In section 4, we describe our
numerical analysis of the correlation functions. Then we present our results,
first for the elastic form factors (section 5) and second for the generalized
form factors (section 6).

##  2 Theoretical formalism

###  2.1 Building light-cone GPDs from space-like matrix elements

GPDs parameterize off-forward matrix elements of quark and gluon operators
with a light-like separation. In the convention of Ji (1998), the leading
twist-2 vector quark GPDs of the nucleon are defined according to:

$\displaystyle F^{q}\left(x,p_{f},p_{i}\right)$
$\displaystyle=\frac{1}{2}\int\frac{{\rm d}z^{-}}{2\pi}e^{ixP^{+}z^{-}}$ (3)
$\displaystyle\quad\times\bra{N\left(p_{f},\lambda_{f}\right)}\bar{\psi}^{q}\left(-\frac{z}{2}\right)\gamma^{+}\hat{W}\left(-\frac{z}{2},\frac{z}{2};A\right)\psi^{q}\left(\frac{z}{2}\right)\ket{N\left(p_{i},\lambda_{i}\right)}\mid_{z^{+}=0,\mathbf{z}_{\perp}=\mathbf{0}_{\perp}}\,,$
$\displaystyle=\frac{1}{2P^{+}}\overline{u}\left(p_{f},\lambda_{f}\right)\left[\gamma^{+}H^{q}\left(x,\xi,t\right)+\frac{i\sigma^{+\nu}q_{\nu}}{2m}E^{q}\left(x,\xi,t\right)\right]u\left(p_{i},\lambda_{i}\right)\,,$
(4)

where $\hat{W}$ is the Wilson line in the fundamental representation, $m$ the
nucleon mass, $\lambda_{i},\lambda_{f}$ are the initial and final state
helicities, and $p_{i},p_{f}$ are the initial and final state momenta that
define the variables

$P\equiv\frac{1}{2}(p+p^{\prime})\,,\ \ q\equiv p^{\prime}-p\,,\ \ t\equiv
q^{2}\,,\ \ \xi\equiv-\frac{q^{+}}{2P^{+}}\,.$ (5)

We use the spinor normalization
$\bar{u}(p,\lambda)u(p,\lambda^{\prime})=2m\delta_{\lambda\lambda^{\prime}}$.

For a Euclidean lattice QCD calculation, the light-cone definition is
unsuitable. Since real-time is inaccessible, we use instead equal-time non-
local operators with a space-like separation. Let us consider more general
fundamental matrix elements defining GPDs of spin-1/2 nucleon without
restriction on $z^{\mu}$:

$M^{\mu}\left(p_{f},p_{i},z\right)=\bra{N\left(p_{f},\lambda_{f}\right)}\mathcal{O}^{\mu}\left(z;A\right)\ket{N\left(p_{i},\lambda_{i}\right)}.$
(6)

This matrix element has multiplicative ultraviolet divergences Ishikawa  _et
al._ (2017) that only depend on the length of the Wilson line, so it is useful
to use instead the ratio:

$\mathcal{M}^{\mu}\left(p_{f},p_{i},z\right)=\frac{M^{\mu}\left(p_{f},p_{i},z\right)}{M^{0}\left(0,0,z\right)}\,,$
(7)

which has a finite continuum limit. In fact, this ratio is renormalization
group invariant (RGI), which means that it is independent of the scheme that
is used to define the renormalized quantities in the numerator and
denominator. This denominator is the same matrix element used in the analogous
ratio of the forward matrix elements in Egerer  _et al._ (2021). In the case
of isovector unpolarized quarks, the denominator also directly normalizes the
Dirac elastic form factor at 1 when $t=0$.

Introducing the Ioffe time $\nu=z\cdot P$, the light-cone GPD definition (3)
can be expressed as:

$\displaystyle F^{q}\left(x,p_{f},p_{i}\right)$
$\displaystyle=\lim_{z\rightarrow(0,z^{-},0_{\perp})}\frac{1}{2P^{+}}\int\frac{{\rm
d}\nu}{2\pi}e^{ix\nu}\mathcal{M}^{+}(p_{f},p_{i},z)\,,$ (8)
$\displaystyle=\lim_{z^{2}\rightarrow 0}\frac{1}{2}\int\frac{{\rm
d}\nu}{2\pi\nu}e^{ix\nu}z_{\mu}\mathcal{M}^{\mu}(p_{f},p_{i},z)\,.$ (9)

This expression is particularly convenient as it is directly Lorentz
invariant. It should be noted that the limit in Eq. (9) suffers from collinear
divergences, but perturbation theory allows us to predict the behavior in
$z^{2}$ of the leading twist contribution when $z$ is small enough compared to
hadronic scales. However, since all our calculations are performed at
$z^{2}<0$, there exists an uncertainty associated with the light-cone limit
$z^{2}\rightarrow 0$.

The matrix element admits the Lorentz decomposition Bhattacharya  _et al._
(2022):

$\displaystyle\mathcal{M}^{\mu}\left(p_{f},p_{i},z\right)$
$\displaystyle=\langle\langle\gamma^{\mu}\rangle\rangle\mathcal{A}_{1}\left(\nu,\xi,t,z^{2}\right)+z^{\mu}\langle\langle\mathds{1}\rangle\rangle\mathcal{A}_{2}\left(\nu,\xi,t,z^{2}\right)+i\langle\langle\sigma^{\mu
z}\rangle\rangle\mathcal{A}_{3}\left(\nu,\xi,t,z^{2}\right)$
$\displaystyle\qquad\qquad+\frac{i}{2m}\langle\langle\sigma^{\mu
q}\rangle\rangle\mathcal{A}_{4}\left(\nu,\xi,t,z^{2}\right)+\frac{q^{\mu}}{2m}\langle\langle\mathds{1}\rangle\rangle\mathcal{A}_{5}\left(\nu,\xi,t,z^{2}\right)$
$\displaystyle\quad+\frac{i}{2m}\langle\langle\sigma^{zq}\rangle\rangle\left[P^{\mu}\mathcal{A}_{6}\left(\nu,\xi,t,z^{2}\right)+q^{\mu}\mathcal{A}_{7}\left(\nu,\xi,t,z^{2}\right)+z^{\mu}\mathcal{A}_{8}\left(\nu,\xi,t,z^{2}\right)\right]\,.$
(10)

The Lorentz structures preceding the amplitudes employ the following
abbreviations: $\sigma^{\mu z}\equiv\sigma^{\mu\rho}z_{\rho}$, $\sigma^{\mu
q}\equiv\sigma^{\mu\rho}q_{\rho}$,
$\sigma^{zq}\equiv\sigma^{\rho\lambda}z_{\rho}q_{\lambda}$, and
$\langle\langle\Gamma\rangle\rangle\equiv\overline{u}\left(p_{f},\lambda_{f}\right)\Gamma
u\left(p_{i},\lambda_{i}\right)$. $\xi$ may be generalized to space-like
separations via:

$\xi=-\frac{q\cdot z}{2P\cdot z}=-\frac{q\cdot z}{2\nu}\,.$ (11)

This expression coincides exactly with the definition of $\xi$ in the light-
cone limit (5), but interestingly the kinematic bound satisfied by the light-
cone $\xi$

$|\xi_{\textrm{light-cone}}|\leq\frac{\sqrt{-t}}{\sqrt{-t+4m^{2}}}\leq 1\,,$
(12)

no longer holds for space-like $z$. In fact, with this definition of $\xi$,
one can choose kinematics such that $\nu=0$ and $\xi$ is infinite. This is,
however, not an issue as all amplitudes and matching relations appear
naturally as functions of $\bar{\nu}=\nu\xi=-q\cdot z/2$, which is always
finite. The decomposition of the GPD matrix elements in terms of $\nu$ and
$\bar{\nu}$ emerges clearly in the double distribution (DD) framework
Radyushkin (2024), as shown in the Lorentz invariant relationship Eq. (55).

The matrix element contracted with $z_{\mu}$ that enters the light-cone limit,
Eq. (9), reads:

$\displaystyle z_{\mu}\mathcal{M}^{\mu}\left(p_{f},p_{i},z\right)$
$\displaystyle=\langle\langle\not{z}\rangle\rangle\mathcal{A}_{1}\left(\nu,\xi,t,z^{2}\right)+z^{2}\langle\langle\mathds{1}\rangle\rangle\mathcal{A}_{2}\left(\nu,\xi,t,z^{2}\right)+0\times\mathcal{A}_{3}\left(\nu,\xi,t,z^{2}\right)$
$\displaystyle\qquad\qquad+\frac{i}{2m}\langle\langle\sigma^{zq}\rangle\rangle\mathcal{A}_{4}\left(\nu,\xi,t,z^{2}\right)-\frac{\nu\xi}{m}\langle\langle\mathds{1}\rangle\rangle\mathcal{A}_{5}\left(\nu,\xi,t,z^{2}\right)$
$\displaystyle\quad+\frac{i}{2m}\langle\langle\sigma^{zq}\rangle\rangle\left[\nu\mathcal{A}_{6}\left(\nu,\xi,t,z^{2}\right)-2\xi\nu\mathcal{A}_{7}\left(\nu,\xi,t,z^{2}\right)+z^{2}\mathcal{A}_{8}\left(\nu,\xi,t,z^{2}\right)\right]\,,$
(13)

where we have used that $\sigma^{\mu\nu}z_{\mu}z_{\nu}=0$. One can immediately
identify that $\mathcal{A}_{3}$ disappears in the contraction and
$\mathcal{A}_{2,8}$ are $\mathcal{O}(z^{2})$ contaminations since in QCD
$\mathcal{A}_{2,8}$ can only have logarithmic colinear divergences in the
limit of $z^{2}\rightarrow 0$. If we introduce Ioffe-time GPDs as:

$\displaystyle\begin{pmatrix}H^{q}\\\
E^{q}\end{pmatrix}(x,\xi,t)=\int\frac{\mathrm{d}\nu}{2\pi}e^{ix\nu}\begin{pmatrix}H^{q}\\\
E^{q}\end{pmatrix}(\nu,\xi,t)\,,$ (14)

then we can identify which amplitudes contribute to the lightcone limit by
comparison with Eq. (4):

$\displaystyle H^{q}(\nu,\xi,t)$ $\displaystyle=\lim_{z^{2}\rightarrow
0}\bigg{[}\mathcal{A}_{1}-\xi A_{5}\bigg{]}\,,$ (15) $\displaystyle
E^{q}(\nu,\xi,t)$ $\displaystyle=\lim_{z^{2}\rightarrow
0}\bigg{[}\mathcal{A}_{4}+\nu\mathcal{A}_{6}-2\xi\nu\mathcal{A}_{7}+\xi
A_{5}\bigg{]}\,,$ (16)

where scale dependence from regulating the collinear divergences at the
$z^{2}\to 0$ limit has been suppressed. We have used the Gordon identity to
distribute the contribution of $\mathcal{A}_{5}$ between $H$ and $E$. This
relation is equivalent to the one quoted in Bhattacharya  _et al._ (2022) up
to choices of conventions. The discussion in Radyushkin (2024) using the DD
formalism Radyushkin (1996a, b, 1997) demonstrates that the $\mathcal{A}_{5}$
amplitude is sensitive to the $D$-term in the limit $z^{2}\rightarrow 0$,
while the other amplitudes $\mathcal{A}_{1,4,6,7}$ build up the double
distributions in a similar fashion as here. The right-hand sides of Eqs. (15)
and (16) are very advantageous since they can be evaluated immediately at any
value of $z^{2}$ with the assurance that they converge towards the correct
light-cone limit for arbitrarily small values of $z^{2}$. However, there is
clearly no unique way to construct a Lorentz invariant object at $z^{2}<0$
which converges analytically towards the correct light-cone limit.

In the case of unpolarized quark PDFs, a single amplitude contributes leading
twist information at lowest powers of $z^{2}$ and $\alpha_{s}$. For helicity
distributions a linear combination of two is required Edwards  _et al._
(2023); Balitsky  _et al._ (2022); Egerer  _et al._ (2022a), not dissimilar to
Eqs. (15) and (16). In the case of GPDs, different choices of the quasi-GPD
were compared in terms of the Lorentz invariant amplitudes Bhattacharya  _et
al._ (2022). In the light-cone limit $z^{2}\to 0$ at fixed $\nu,\xi$, these
definitions converge to the same final GPD. The practical differences between
the definitions for realistic scales would be dominated by twist-3
contributions Braun (2023); Radyushkin (2024) in the amplitudes that come from
terms proportional to the transverse momentum transfer
$\Delta_{\mu}^{T}=q^{\mu}+2\xi P^{\mu}$. When considering the light-cone
matrix element in Eq. (9) with a generic $\gamma^{\mu}$, when $\mu\neq+$ the
dominant deviations from the leading twist matrix elements again occur in the
OPE from twist-3 operators that are total derivatives of leading twist
operators or from genuine quark-gluon-antiquark correlations Anikin  _et al._
(2000); Penttinen  _et al._ (2000); Belitsky and Mueller (2000); Radyushkin
and Weiss (2000, 2001); Kivel and Polyakov (2001). In Braun (2023), the effect
of such twist-3 contributions is to partially restore the translational
symmetry of the off-light-cone matrix element, which is broken in a leading
twist calculation. As shown in Radyushkin (2024) adding the
$\mathcal{A}_{6,7}$ terms to Eq. (16) is explicitly canceling the contribution
to two such twist-3 contributions. If one could determine higher twist
combinations, of any order and from the same or a different calculation, then
they could be added similarly to Eq. (15) and (16) to improve the light-cone
limit $z^{2}\to 0$, even if not all twist-3 contribution are separable using
Lorentz invariant functions. An evaluation of the plausible size and
analytical properties of the genuine higher-twist power corrections is carried
out in Braun  _et al._ (2024), where the corrections are found to be fairly
mild when using the ratio normalization of Eq. (7) in the range of Ioffe time
relevant for this study.

###  2.2 Perturbative matching

Extending the pseudo-PDF’s one-loop matching Radyushkin (2018); Zhang  _et
al._ (2018); Izubuchi  _et al._ (2018), the matching of unpolarized non-
singlet Ioffe-time GPDs in the $\overline{\rm MS}$ scheme to $z^{2}$-dependent
pseudo-distributions in coordinate space was first presented at one-loop
accuracy in Radyushkin (2019) for the pion’s analog of $\mathcal{A}_{1}$. For
clarity, we denote with a bar the $\overline{\rm MS}$ quantity. The matching
is given by:

$F\left(\nu,\xi,t,z^{2}\right)=\overline{F}\left(\nu,\xi,t,\mu^{2}\right)-\frac{\alpha_{s}C_{F}}{2\pi}\left\\{\ln\left(-\frac{e^{2\gamma_{E}+1}}{4}z^{2}\mu^{2}\right)B+L\right\\}\otimes\overline{F}+\mathcal{O}\left(z^{2}\Lambda_{\rm
QCD}^{2}\right)+\mathcal{O}\left(z^{2}t\right),$ (17)

where:

$\displaystyle B\otimes\overline{F}$ $\displaystyle=\int_{0}^{1}{\rm d}\alpha\
\left[\left(\frac{2\alpha}{1-\alpha}\right)_{+}\cos\left(\bar{\alpha}\xi\nu\right)+\frac{\sin\left(\bar{\alpha}\xi\nu\right)}{\xi\nu}-\frac{\delta\left(1-\alpha\right)}{2}\right]\overline{F}\left(\alpha\nu,\xi,t,\mu^{2}\right)$
(18) $\displaystyle L\otimes\overline{F}$ $\displaystyle=\int_{0}^{1}{\rm
d}\alpha\
\left[4\left(\frac{\ln\left(1-\alpha\right)}{1-\alpha}\right)_{+}\cos\left(\bar{\alpha}\xi\nu\right)-2\frac{\sin\left(\bar{\alpha}\xi\nu\right)}{\xi\nu}+\delta\left(1-\alpha\right)\right]\overline{F}\left(\alpha\nu,\xi,t,\mu^{2}\right),$
(19)

with $\bar{\alpha}=1-\alpha$. The standard plus-prescription is defined as

$\int_{0}^{1}{\rm d}\alpha\ G\left(\alpha\right)_{+}f\left(\alpha
x\right)=\int_{0}^{1}{\rm d}\alpha\ G\left(\alpha\right)\left[f\left(\alpha
x\right)-f\left(x\right)\right]\,.$ (20)

Here $B$ is simply the ordinary leading-order evolution kernel of non-singlet
GPDs, except expressed in Ioffe time instead of $x$. Compared to its
expression in $x$-space, the evolution operator is much simpler. For instance,
it is analytical in $\xi\nu$, whereas the $x$-space evolution operator has
non-analytic structures in $\xi/x$. Furthermore, the evolution/matching of
Ioffe-time GPDs shares with Ioffe-time PDFs the characteristic of only
involving smaller Ioffe times (the $\alpha$ variable in the integral runs from
0 to 1). Thus, we do not need to know the entire Ioffe-time dependence of the
GPD for evolution. In contrast, the $x$-space evolution of GPDs will typically
require knowledge of the entire $x$-dependence of the GPD. It is only when
$|x|>|\xi|$ that a partial knowledge of the $x$-dependence of GPDs is enough
to evolve them111Depending on which kinematic region in $(x,\xi)$ of the GPDs
one exactly knows, it is possible to imagine a somewhat cumbersome strategy
where one would use double distributions as an intermediate step before
performing the evolution Dall’Olio  _et al._ (2024). It is however potentially
numerically challenging, and certainly less straightforward than the evolution
of a Ioffe-time GPD.. As a consequence, it seems that Ioffe-time GPDs are
generally speaking a much friendlier object to manipulate than their
$x$-dependent counterpart, even though they contain the same information if
the entire range in $\nu$ and $x$ is known.

The kernels governing the evolution, Eq. (18), and scheme matching, Eq. (19),
connect the scales $\mu^{2}$ in $\overline{\rm{MS}}$ and $z^{2}$ at fixed
values of Ioffe time, skewness, and momentum transfer. In the limit
$\xi\rightarrow 0$, the kernels reduce to the PDF case. It is interesting to
observe the effects of non-zero skewness. To isolate the skewness dependence
attached purely to the perturbative matching, we introduce a $\xi$-independent
GPD resembling a normalized isovector PDF:

$\overline{F}\left(x,\xi,t,\mu^{2}\right)=\frac{15}{16}\,x^{-0.5}\left(1-x\right)^{2}\,.$
(21)

The associated Ioffe-time GPD

$\overline{F}\left(\nu,\xi,t,\mu^{2}\right)=\int_{-1}^{1}\mathrm{d}x\,e^{-ix\nu}\overline{F}\left(x,\xi,t,\mu^{2}\right)$
(22)

is presented as the red curve in the right panels of Fig. 3 for the real part
and Fig. 4 for the imaginary part. The blue curves on the same plots represent
the matched pseudo-GPD $F(\nu,\xi,t,z^{2})$ where we have assumed that $\mu=2$
GeV, $\alpha_{s}=0.28$ and $z=a=0.094$ fm. This yields
$\ln(-z^{2}\mu^{2}e^{2\gamma_{E}+1}/4)\approx 0.67$. It can be observed that
the matching effects are quite small, especially for larger values of $\xi$.

Figure 3: (left) Convolution of
$\mathrm{Re}\,\overline{F}\left(\nu,\xi,t,\mu^{2}\right)$ obtained from Eq.
(21) with the evolution kernel $B$ and the matching kernel $L$ for various
values of $\xi$. (right) The original Ioffe-time GPD (red) and the result of
the matching for $\mu=2$ GeV, $\alpha_{s}=0.28$, $z=0.094$ fm and various
values of $\xi$. Figure 4: Same as Fig. 3, but for the imaginary part.

To understand how evolution and matching separately contribute to the effects
of non-zero skewness, we show the separate convolutions $B\otimes\overline{F}$
and $L\otimes\overline{F}$ on the left panels of Figs. 3 and 4. A striking
feature is that the curves for various skewnesses seem to become similar as
$\nu$ gets larger, but the smaller the skewness, the longer they take to fall
towards this “universal” limit. Although it is difficult to assess precisely,
a general rule can be observed from the plots that the universal limit is
reached when $\nu\approx 5/\xi$, and it is unclear whether the curve for
$\xi=0$ participates in this universal limit at all. We have checked that this
remains consistently true for much larger values of $\nu$.

In fact, this is not surprising. The $x$-space evolution kernel has a cusp
when $x=\xi$, which leaves us to expect some kind of trouble in Ioffe time
when $\nu=\mathcal{O}(1/\xi)$. As we have noticed before, the relevant
variable that appears consistently in our calculations is not $\xi$, but
rather $\nu\xi$. A similar study is conducted in $x$-space in Bertone  _et
al._ (2022) for the evolution kernel. It demonstrates the formation of non-
analytic behavior at $x=\xi$ through the evolution operator and shows how
evolution effects are considerably increased in this region (see Fig. 3 in
that paper, for instance). We observe here the same effect in Ioffe-time
space, although in the range of $\nu$ available for realistic calculations,
the $\xi$ dependence of the evolution kernels is quite benign.

Additionally, one can notice in Fig. 4 that the small $\nu$ behavior of the
imaginary part of the convolutions is independent of the skewness. As we will
see in the next section, this stems from the fact that the first derivative
with respect to $\nu$ of the imaginary part of the non-singlet Ioffe-time GPD
at $\nu=0$ is independent of $\xi$ due to the polynomiality property of GPDs.
This also explains why the small $\nu$ behavior of the real part of the
convolutions in Fig. 3 is of order $\mathcal{O}(\nu^{2})$, but this time
dependent on $\xi$.

Finally, a feature of the matching already noticed before in the case of PDFs
or Good Lattice Cross Sections (e.g. Joó  _et al._ (2019a); Sufian  _et al._
(2020); Dutrieux  _et al._ (2024)) is the fact that the evolution and matching
convolutions have generally opposite signs in the regime of $\nu$ where data
exists. The result is a reduction or even near cancellation of the dominant
perturbative effect. Of course, we did not include any initial skewness
dependence in the model, so such results are qualitative at best. It will be
interesting to study if these canceling effects continue to hold with higher
orders in perturbation theory and more sophisticated resummation techniques.

In Dutrieux  _et al._ (2024), the perturbatively calculated kernel for the
resummed $z^{2}$ PDF step-scaling from one scale to another was studied and
compared to DGLAP evolution. The kernel for the forward pseudo-PDF was found
to be smaller than the $\overline{\rm MS}$ DGLAP evolution kernel. Further
comparisons should be made in future studies to other scheme-dependent
methods, such as RI-Mom renormalization, where the scale-independent structure
is different.

###  2.3 Moments of GPDs

By differentiating the defining relationship between the Ioffe-time
distribution and the $x$-space distribution, Eq. (22), one finds that:

$\frac{d^{n}}{d\nu^{n}}F(\nu,\xi,t,z^{2})\bigg{|}_{\nu=0}=(-i)^{n}\int_{-1}^{1}dx\,x^{n}F(x,\xi,t,z^{2})\,.$
(23)

We introduce the short-hand notation:

$F_{n}(\xi,t,z^{2})\equiv\int_{-1}^{1}dx\,x^{n-1}F(x,\xi,t,z^{2})\,,$ (24)

where $n{-}1$ is a widely encountered convention. Then:

$F(\nu,\xi,t,z^{2})=\sum_{n=0}^{\infty}\frac{(-i\nu)^{n}}{n!}F_{n+1}(\xi,t,z^{2})\,.$
(25)

Introducing the matching relation:

$F(\nu,\xi,t,z^{2})=\int_{0}^{1}d\alpha\,K(\alpha,\xi\nu,\alpha_{s},z^{2}\mu^{2})\overline{F}(\alpha\nu,\xi,t,\mu^{2})\,,$
(26)

and differentiating it $n$ times yields:

$\displaystyle(-i)^{n}F_{n+1}(\xi,t,z^{2})$
$\displaystyle=\sum_{k=0}^{n}\begin{pmatrix}n\\\
k\end{pmatrix}\int_{0}^{1}d\alpha\,\frac{d^{n-k}}{d\nu^{n-k}}K(\alpha,\xi\nu,\alpha_{s},z^{2}\mu^{2})\bigg{|}_{\nu=0}\frac{d^{k}}{d\nu^{k}}\overline{F}(\alpha\nu,\xi,t,\mu^{2})\bigg{|}_{\nu=0}\,,$
(27)
$\displaystyle=\sum_{k=0}^{n}\overline{F}_{k+1}(\xi,t,\mu^{2})\begin{pmatrix}n\\\
k\end{pmatrix}\int_{0}^{1}d\alpha\,(-i\alpha)^{k}\frac{d^{n-k}}{d\nu^{n-k}}K(\alpha,\xi\nu,\alpha_{s},z^{2}\mu^{2})\bigg{|}_{\nu=0}$
(28)
$\displaystyle\equiv(-i)^{n}\sum_{k=0}^{n}\overline{F}_{k+1}(\xi,t,\mu^{2})\xi^{n-k}c_{n+1,n-k}(\alpha_{s},z^{2}\mu^{2})\,,$
(29)

where

$c_{n+1,k}(\alpha_{s},z^{2}\mu^{2})\equiv i^{k}\begin{pmatrix}n\\\
k\end{pmatrix}\int_{0}^{1}d\alpha\,\alpha^{n-k}\frac{d^{k}}{d(\xi\nu)^{k}}K(\alpha,\xi\nu,\alpha_{s},z^{2}\mu^{2})\bigg{|}_{\nu=0}\,.$
(30)

In this paper we are concerned with the spin-1/2 unpolarized GPDs, so
$K(\alpha,\xi\nu)$ is even in $\xi\nu$, and thus only the $c_{n+1,k}$ where
$k$ is even are non-zero. In particular, $c_{n+1,k}$ is always real.

We now list the first few values of $c_{n+1,k}$ at order
$\mathcal{O}(\alpha_{s})$, using the convention that

$L=\ln\left(-\frac{e^{2\gamma_{E}+1}}{4}z^{2}\mu^{2}\right)\,:$ (31)

skewness independent | quadratic skewness
---|---
$c_{1,0}=1$ |
$c_{2,0}=1-\frac{\alpha_{s}C_{F}}{2\pi}\left(-\frac{4}{3}L+\frac{14}{3}\right)$ |
$c_{3,0}=1-\frac{\alpha_{s}C_{F}}{2\pi}\left(-\frac{25}{12}L+\frac{47}{6}\right)$ | $c_{3,2}=-\frac{\alpha_{s}C_{F}}{2\pi}\left(\frac{5}{12}L-\frac{7}{6}\right)$
$c_{4,0}=1-\frac{\alpha_{s}C_{F}}{2\pi}\left(-\frac{157}{60}L+\frac{931}{90}\right)$ | $c_{4,2}=-\frac{\alpha_{s}C_{F}}{2\pi}\left(\frac{11}{20}L-\frac{53}{30}\right)$

The coefficient of $-\alpha_{s}C_{F}L/2\pi$ in $c_{n+1,0}$ is simply the
leading-order DGLAP anomalous dimension. For $\xi\neq 0$, the GPD Mellin
moments mix among themselves through the coefficients $c_{n+1,k}$ with $k>0$.
It is convenient to write the matching in matrix form:

$\begin{pmatrix}F_{1}(\xi,t,z^{2})\\\ F_{2}(\xi,t,z^{2})\\\
F_{3}(\xi,t,z^{2})\\\ F_{4}(\xi,t,z^{2})\\\
\cdots\end{pmatrix}=\begin{pmatrix}c_{1,0}&0&0&0&\\\ 0&c_{2,0}&0&0&\\\
\xi^{2}c_{3,2}&0&c_{3,0}&0&\\\ 0&\xi^{2}c_{4,2}&0&c_{4,0}&\\\
&&&&\ddots\end{pmatrix}\times\begin{pmatrix}\overline{F}_{1}(\xi,t,\mu^{2})\\\
\overline{F}_{2}(\xi,t,\mu^{2})\\\ \overline{F}_{3}(\xi,t,\mu^{2})\\\
\overline{F}_{4}(\xi,t,\mu^{2})\\\ \cdots\end{pmatrix}\,.$ (32)

Then, by isolating the logarithmic part of the kernel, one can introduce an
expression resumed at leading-logarithmic accuracy which generalizes the well-
known PDF case (for the LL resummation applied to the pseudo-PDF matching, see
Gao  _et al._ (2021)) by using the appropriate set of moments diagonalizing
the GPD evolution Ohrndorf (1982); Chase (1980):

$\begin{pmatrix}\overline{F}_{1}(\xi,t,\mu^{2})\\\
\overline{F}_{2}(\xi,t,\mu^{2})\\\ \overline{F}_{3}(\xi,t,\mu^{2})\\\
\overline{F}_{4}(\xi,t,\mu^{2})\\\
\cdots\end{pmatrix}=G^{-1}(\xi)A(\mu^{2},z^{2})G(\xi)\times\bigg{(}I-B(\xi,z^{2})\bigg{)}\times\begin{pmatrix}F_{1}(\xi,t,z^{2})\\\
F_{2}(\xi,t,z^{2})\\\ F_{3}(\xi,t,z^{2})\\\ F_{4}(\xi,t,z^{2})\\\
\cdots\end{pmatrix}\,,$ (33)

where $A$ is the matrix of the leading-logarithmic resummed DGLAP evolution

$A(\mu^{2},z^{2})=\begin{pmatrix}1&0&0&0&\\\ 0&r^{-4/3}&0&0&\\\
0&0&r^{-25/12}&0&\\\ 0&0&0&r^{-157/60}&\\\ &&&&\ddots\end{pmatrix}\textrm{
with
}r=\left(\frac{\alpha_{s}(-1/z^{2})}{\alpha_{s}(\mu^{2})}\right)^{C_{F}/(2\pi\beta_{0})}\,.$
(34)

Here $\beta_{0}=(33-2n_{f})/(12\pi)\approx 0.716$, $G$ is the matrix
diagonalizing the evolution of moments of GPDs (a specific type of Gegenbauer
moments for leading-order evolution),

$G(\xi)=\begin{pmatrix}1&0&0&0&\\\ 0&1&0&0&\\\ -\xi^{2}/5&0&1&0&\\\
0&-3\xi^{2}/7&0&1&\\\ &&&&\ddots\end{pmatrix}\,,$ (35)

and $I$ is the identity matrix of the relevant dimension. The non-logarithmic
part of the matching kernel is

$B(\xi,z^{2})=\begin{pmatrix}0&0&0&0&\\\
0&-a_{s}\left(-\frac{4}{3}l+\frac{14}{3}\right)&0&0&\\\
-\xi^{2}a_{s}\left(\frac{5}{12}l-\frac{7}{6}\right)&0&-a_{s}\left(-\frac{25}{12}l+\frac{47}{6}\right)&0&\\\
0&-\xi^{2}a_{s}\left(\frac{11}{20}l-\frac{53}{30}\right)&0&-a_{s}\left(-\frac{157}{60}l+\frac{931}{30}\right)&\\\
&&&&\ddots\\\ \end{pmatrix}\,,$ (36)

where $a_{s}=\alpha_{s}(-1/z^{2})C_{F}/(2\pi)$ and
$l=\ln(e^{2\gamma_{E}+1}/4)\approx 0.768$.

We have taken the “natural” scale of non-local lattice calculations to be
$-1/z^{2}$, but one can easily choose to vary this scale. It is important to
notice that the expression of Eq. (33) is far from unique, and several
sensible perturbative choices would result in different formulations,
differing by contributions of order $\mathcal{O}(\alpha_{s}^{k}L^{m})$ where
$k\geq 2$ and $m\leq k-1$. In the simpler context of PDFs, more complicated
choices of perturbative approach and some of their properties were explored in
Dutrieux  _et al._ (2024), where the existence of this ambiguity motivates the
desire to compute at least the logarithmic part of the kernel directly on the
lattice.

Now we introduce the polynomiality property of GPDs. It is useful here to
separate in the analysis the $D$-term from the rest of the GPD. As we have
already mentioned, the $D$-term is contained in a separate Lorentz amplitude
$\mathcal{A}_{5}$ Radyushkin (2024), while the rest of the $H$ and $E$ GPDs is
constructed out of a combination of the $\mathcal{A}_{1,4,6,7}$ amplitudes.
Furthermore, the matching does not mix the $D$-term with the rest of the GPD.
To see it, we notice that the polynomiality property writes:

$\overline{F}_{n+1}(\xi,t,\mu^{2})=\sum_{\begin{subarray}{c}k=0\\\ k\textrm{
even}\end{subarray}}^{n}\xi^{k}A_{n+1,k}(t,\mu^{2})+\textrm{mod}(n,2)\xi^{n+1}C_{n+1}(t,\mu^{2})\,.$
(37)

where

$C_{n+1}(t,\mu^{2})=\int_{-1}^{1}d\alpha\,\alpha^{n}D(\alpha,t,\mu^{2})\,$
(38)

are the moments of the $D$-term, while the $A_{n,k}$ do not depend on the
$D$-term. Then, the matching formula of Eq. (29) yields:

$\displaystyle\sum_{k=0}^{n}\bigg{[}\sum_{\begin{subarray}{c}i=0\\\ i\textrm{
even}\end{subarray}}^{k}\xi^{i}A_{k+1,i}(t,\mu^{2})+\textrm{mod}(k,2)\xi^{k+1}C_{k+1}(t,\mu^{2})\bigg{]}\xi^{n-k}c_{n+1,n-k}$
$\displaystyle\hskip
50.0pt=\sum_{k=0}^{n}\bigg{[}\sum_{\begin{subarray}{c}i=0\\\ i\textrm{
even}\end{subarray}}^{k}\xi^{n+i-k}A_{k+1,i}(t,\mu^{2})c_{n+1,n-k}\bigg{]}+\xi^{n+1}\textrm{mod}(n,2)\sum_{k=0}^{n}C_{k+1}(t,\mu^{2})c_{n+1,n-k}\,.$
(39)

The first term only contains powers of $\xi$ up to $\xi^{n}$, while the second
term has the same form as a $D$-term. This means that the $z^{2}$ dependent
pseudo-GPD obeys a similar polynomiality property as the $\overline{MS}$ GPD,
and that the $D$-term of the pseudo-GPD can be matched to the $D$-term of the
$\overline{MS}$ GPD independently from the rest of the GPD. It is convenient
to write the polynomiality property of the rest of the GPD, noted here as
$\overline{F}_{n+1,\textrm{ no D-term}}$ as a simple matrix multiplication:

$\begin{pmatrix}\overline{F}_{1,\textrm{ no D-term}}(\xi,t,\mu^{2})\\\
\overline{F}_{2,\textrm{ no D-term}}(\xi,t,\mu^{2})\\\
\overline{F}_{3,\textrm{ no D-term}}(\xi,t,\mu^{2})\\\
\overline{F}_{4,\textrm{ no D-term}}(\xi,t,\mu^{2})\\\
\cdots\end{pmatrix}=\begin{pmatrix}1&0&0&0&0&0&\\\ 0&1&0&0&0&0&\\\
0&0&1&\xi^{2}&0&0&\\\ 0&0&0&0&1&\xi^{2}&\\\
&&&&&&\ddots\end{pmatrix}\times\begin{pmatrix}A_{1,0}(t,\mu^{2})\\\
A_{2,0}(t,\mu^{2})\\\ A_{3,0}(t,\mu^{2})\\\ A_{3,2}(t,\mu^{2})\\\
A_{4,0}(t,\mu^{2})\\\ A_{4,2}(t,\mu^{2})\\\ \cdots\end{pmatrix}\,,$ (40)

and the independent matching of this piece is obtained by inserting this
expression into Eq. (33), and evaluating the result for a sufficient range of
different values of $\xi$ to obtain invertibility (one value of $\xi$ for
$F_{1,2}$ and two values for $F_{3,4}$ are enough).

As for the $D$-term, from:

$\displaystyle\overline{H}(x,\xi,t,\mu^{2})=\overline{H}_{\textrm{ no
D-term}}(x,\xi,t,\mu^{2})+\textrm{sgn}(\xi)\overline{D}\left(\frac{x}{\xi},t,\mu^{2}\right)\Theta(|\xi|-|x|)\,,$
(41)

one gets:

$\displaystyle\overline{H}(\nu,\xi,t,\mu^{2})=\overline{H}_{\textrm{ no
D-term}}(\nu,\xi,t,\mu^{2})+\xi\int_{-1}^{1}d\alpha\,e^{-i\alpha\xi\nu}\overline{D}(\alpha,t,\mu^{2})$
(42)

and comparing with Eq. (15) gives that:

$\displaystyle\lim_{z^{2}\rightarrow 0}\mathcal{A}_{5}(\nu,\xi,t,z^{2})$
$\displaystyle=-\int_{-1}^{1}d\alpha\,e^{-i\alpha\xi\nu}D(\alpha,t,z^{2})\,,$
(43) $\displaystyle=i\xi\nu\int_{-1}^{1}d\alpha\,\alpha
D(\alpha,t,z^{2})-i(\xi\nu)^{3}\int_{-1}^{1}d\alpha\,\alpha^{3}D(\alpha,t,z^{2})+...\,,$
(44)

where we have used the oddness in $\alpha$ of the $D$-term which implies that
$\mathcal{A}_{5}$ should be purely imaginary in the limit $z^{2}\rightarrow
0$. Deviations from this might arise both because of lattice systematics and
higher twist contaminations due to the non-zero value of $z^{2}$. Finally, the
introduction of the matching yields:

$\displaystyle\lim_{z^{2}\rightarrow
0}\mathcal{A}_{5}(\nu,\xi,t,z^{2})=i\xi\nu
C_{2}(t,\mu^{2})c_{2,0}(\alpha_{s},z^{2}\mu^{2})$ $\displaystyle\hskip
70.0pt-i(\xi\nu)^{3}\bigg{[}C_{4}(t,\mu^{2})c_{4,0}(\alpha_{s},z^{2}\mu^{2})+C_{2}(t,\mu^{2})c_{4,2}(\alpha_{s},z^{2}\mu^{2})\bigg{]}+...$
(45)

###  2.4 Discrete and Lorentz symmetries of Ioffe-time GPDs

Let us first understand the symmetries of the amplitudes $\mathcal{A}_{k}$
entering the Lorentz decomposition of Eq. (10) with respect to $\nu\to-\nu$
and $\xi\to-\xi$. In this section, for brevity we choose a gauge where the
link connecting the quark and anti-quark fields is unity. Since it is a gauge
invariant quantity the discrete symmetries will not be different in other
gauges. The parity (P), charge conjugation (C), and time reversal (T)
properties of the operator
$O^{\mu}(z)=\bar{\psi}(-\frac{z}{2})\gamma^{\mu}\psi(\frac{z}{2})$ in Eq. (3)
and its (anti)-hermitian combination
$O^{\mu}_{\pm}(z)=\frac{1}{2}\big{(}O^{\mu}(z)\pm O^{\mu}(-z)\big{)}$ are
given by:

$\displaystyle PO^{\mu}(z)P=$
$\displaystyle(-1)^{\mu}\bar{\psi}\left(\frac{z}{2}\right)\gamma^{\mu}\psi\left(-\frac{z}{2}\right)=$
$\displaystyle(-1)^{\mu}O^{\mu}(-z)\quad$ $\displaystyle\to\quad
PO_{\pm}^{\mu}(z)P=$ $\displaystyle\pm(-1)^{\mu}O_{\pm}^{\mu}(z)$
$\displaystyle CO^{\mu}(z)C=$
$\displaystyle-\bar{\psi}\left(\frac{z}{2}\right)\gamma^{\mu}\psi(-\frac{z}{2})=$
$\displaystyle-O^{\mu}(-z)\quad$ $\displaystyle\to\quad CO_{\pm}^{\mu}(z)C=$
$\displaystyle\mp O_{\pm}^{\mu}(z)$ $\displaystyle TO^{\mu}(z)T=$
$\displaystyle(-1)^{\mu}\bar{\psi}\left(-\frac{z}{2}\right)\gamma^{\mu}\psi(\frac{z}{2})=$
$\displaystyle(-1)^{\mu}O^{\mu}(z)\quad$ $\displaystyle\to\quad
TO^{\mu}_{\pm}(z)T=$ $\displaystyle(-1)^{\mu}O^{\mu}_{\pm}(z)\,,$ (46)

where $(-1)^{\mu}=(1,-1,-1,-1)$ and the operators have the following
properties with respect to hermitian conjugation:

$O^{\mu}(z)^{\dagger}=\bar{\psi}\left(\frac{z}{2}\right)\gamma^{\mu}\psi\left(-\frac{z}{2}\right)=O^{\mu}(-z)\quad\to\quad
O^{\mu}_{\pm}(z)^{\dagger}=\pm O^{\mu}_{\pm}(z)\,.$ (47)

By using the combined $PT$ symmetries, the matrix element in Eq. (6), written
now with explicit helicity arguments, obeys

$M^{\mu}(p_{f},p_{i},z,\lambda_{f},\lambda_{i})=M^{\mu}(p_{i},p_{f},z,-\lambda_{i},-\lambda_{f})\,,$
(48)

which swaps the sign of $\xi$. The relation between amplitudes with $\pm\xi$
can be determined term by term from this relationship once the changes in the
kinematic factors in the Lorentz decomposition of Eq. (10) are known. With the
spinors employed in this study, the bilinears
$\bar{u}(p_{f},\lambda_{f})\Gamma u(p_{i},\lambda_{i})$ equal
$\lambda_{f}\lambda_{i}[\bar{u}(p_{f},-\lambda_{f})\Gamma
u(p_{i},-\lambda_{i})]^{\ast}$ for $\Gamma=1,\gamma^{\mu},\sigma^{\mu\nu}$.
Note that all structures with the tensor bilinear are accompanied by an
additional factor of $i$ that generates an additional sign when complex
conjugated. Under these transformations, we can identify the following
relation obeyed by the kinematic factors labeled using $k\in\\{1,..,8\\}$:

$K_{k}^{\mu}(p_{f},p_{i},z,\lambda_{f},\lambda_{i})=\lambda_{f}\lambda_{i}Z^{PT}_{k}K_{k}^{\mu}(p_{i},p_{f},z,-\lambda_{i},-\lambda_{f})\,,$
(49)

where $Z^{PT}_{k}$ is the parity of the amplitudes $\mathcal{A}_{k}$ under
$\xi\to-\xi$, which takes the values $Z^{PT}_{3,5,7}=-1$ and $Z^{PT}_{k\neq
3,5,7}=+1$.

Next, consider the impact of hermiticity:

$M^{\mu}(p_{f},p_{i},z,\lambda_{f},\lambda_{i})^{\ast}=M^{\mu}(p_{i},p_{f},-z,\lambda_{i},\lambda_{f})\,,$
(50)

which swaps the sign of $\xi$ and $\nu$ simultaneously. The parity of the
kinematic factors under hermiticity is given by $Z^{H}_{2,5,6,8}=-1$ and
$Z^{H}_{k\neq 2,5,6,8}=+1$. Therefore, the amplitudes are related to their
complex conjugates as:

$\mathcal{A}_{k}(\nu,\xi)=Z^{H}_{k}\mathcal{A}_{k}(-\nu,-\xi)^{\ast}\,.$ (51)

Combining this result with the previous rule, we find that the real part of
the amplitudes changes with the sign of $\nu$ following
$Z^{\nu}_{k}=Z^{H}_{k}Z^{PT}_{k}$ which is $Z^{\nu}_{k}=-1$ for $k=2,3,6,7,8$
and $Z^{\nu}_{k}=+1$ otherwise. The opposite is true for the imaginary parts.

To understand the implications of this behavior, consider the (anti-)hermitian
operators $O_{\pm}^{\mu}(z)$, whose matrix elements are purely real or
imaginary in the forward case, $(p_{f},\lambda_{f})=(p_{i},\lambda_{i})$. In
the off-forward case the matrix elements
$M_{\pm}^{\mu}(p_{f},p_{i},z)=\frac{1}{2}\big{(}M^{\mu}(p_{f},p_{i},z)\pm
M^{\mu}(p_{f},p_{i},-z)\big{)}$ are given by

$\displaystyle M_{\pm}^{\mu}(p_{f},p_{i},z)$
$\displaystyle=\frac{1}{2}\sum_{k}\left[K_{k}^{\mu}(p_{f},p_{i},z)\mathcal{A}_{k}(\nu,\xi,t,z^{2})\pm
K_{k}^{\mu}(p_{f},p_{i},-z)\mathcal{A}_{k}(-\nu,\xi,t,z^{2})\right]$ (52)
$\displaystyle=\frac{1}{2}\sum_{k}K_{k}^{\mu}(p_{f},p_{i},z)\left[\mathcal{A}_{k}(\nu,\xi,t,z^{2})\pm
Z^{z}_{k}\mathcal{A}_{k}(\nu,\xi,t,z^{2})^{\ast}\right]\,,$ (53)

where $Z^{z}_{k}$ is the parity of the kinematic factors in the sign of $z$.
As can be seen, the hermitian and anti-hermitian operators expose the real and
imaginary components of the amplitudes, but the relation depends on the
amplitude. The phases in the matrix elements are entirely controlled by the
kinematic factors and choices in defining spinors within the calculations.

Finally, let us observe how these discrete symmetries interact with the
Lorentz symmetry, which materializes through the polynomiality property, Eq.
(37). The feature of polynomiality is an intricate relationship in the $\xi$
dependence of the moments of GPDs. It arises naturally in the double
distribution representation, which describes GPDs in $x$-space as (we use only
the $D$-term independent part here):

$H(x,\xi)=\int_{-1}^{1}\mathrm{d}\beta\int_{-1+|\beta|}^{1-|\beta|}\mathrm{d}\alpha\,\delta(x-\beta-\alpha\xi)h(\beta,\alpha)\,.$
(54)

We have omitted the $t$-dependence and the scale of the double distribution
and GPD. It is easy to observe that the parity in $\alpha$ and $\beta$ of the
double distribution controls respectively the $\xi$ and $x$ parity of the GPD.
The relation between the $(\nu,\xi)$ space of Ioffe-time distributions and
$(\beta,\alpha)$ of double distributions is given by integrals of the form:

$\mathcal{A}(\nu,\bar{\nu}=\xi\nu)=\int_{-1}^{1}\mathrm{d}\beta\,e^{i\nu\beta}\int_{-1+|\beta|}^{1-|\beta|}\mathrm{d}\alpha\,e^{i\bar{\nu}\alpha}h(\beta,\alpha)\,.$
(55)

Therefore, definite parity in $\alpha$ of the double distribution is
equivalent to the same definite parity in $\xi$ of the amplitude
$\mathcal{A}$. Provided the double distribution $h$ is real, evenness in
$\beta$ is equivalent to $\mathcal{A}$ even in $\nu$ and real, whereas oddness
in $\beta$ is equivalent to $\mathcal{A}$ odd in $\nu$ and imaginary.

##  3 Numerical Implementation

Previous calculations from the HadStruc collaboration isolated the isovector
twist-2 quark PDFs Orginos  _et al._ (2017); Karpie  _et al._ (2018); Joó  _et
al._ (2019a, b, 2020); Egerer  _et al._ (2021); Edwards  _et al._ (2023);
Egerer  _et al._ (2022b) and unpolarized and helicity gluon PDFs Khan  _et
al._ (2021); Egerer  _et al._ (2022a) of the nucleon using a combination of
distillation Peardon  _et al._ (2009) and the pseudo-distribution formalism.
We use the same gauge ensemble, featuring an $m_{\pi}=358\ {\rm MeV}$ pion
mass within a $32^{3}\times 64$ lattice volume with an $a=0.094\ {\rm fm}$
lattice spacing–we denote this ensemble as a094m358. The a094m358 ensemble is
an isotropic $2+1$ flavor Wilson clover fermion ensemble generated by the
JLab/W&M/LANL/MIT collaboration Edwards  _et al._ (2016), with the strange
quark fixed to its physical value and Wilson clover fermion sea quarks. A
total of $348$ configurations with four distillation “sources” evenly
separated in the time extent per configuration (detailed below) are utilized
to determine all matrix elements. Aspects of the a094m358 ensemble are
summarized in Table 1–further details can be found in Refs. Yoon  _et al._
(2016, 2017).

ID |  $a$ (fm) |  $m_{\pi}$ (MeV) | $\beta$ | $m_{\pi}L$ | $L^{3}\times N_{T}$ | $N_{\rm cfg}$ | $N_{\rm srcs}$ | ${\rm rk}\left(\mathcal{D}\right)$
---|---|---|---|---|---|---|---|---
a094m358 | $0.094(1)$ | $358(3)$ | $6.3$ | $5.4$ | $32^{3}\times 64$ | $348$ | $4$ | $64$
Table 1: Lattice spacing, pion mass, inverse coupling, control parameter
$m_{\pi}L$ for finite-volume effects, lattice volume. The total number of
configurations $N_{\rm cfg}$, temporal sources $N_{\rm srcs}$ per
configuration, and rank of the distillation space ${\rm
rk}\left(\mathcal{D}\right)$ are also given.

To isolate the bare isovector matrix elements of the space-like quark bilinear
in Eq. (3) we require two-point functions and connected three-point functions
with the insertion of the operator
$\overline{\psi}\left(-\frac{z}{2}\right)\gamma^{\mu}W^{(f)}\left(-\frac{z}{2},\frac{z}{2};A\right)\frac{\tau^{3}}{2}\psi\left(\frac{z}{2}\right)$,
where the flavor isovector projector $\tau^{3}/2$ is included. This projector
eliminates statistically noisy disconnected diagrams and gives the difference
of the $u$ and $d$ quark distributions. Placing the operator symmetric about
the origin, with quark and anti-quark fields at $(z/2)$ and $(-z/2)$, ensures
the operator has transparent $C,P,T$ and hermiticity properties, reviewed in
Section 2.4.

Inserting complete sets of energy eigenstates, the spectral content of the
two- and three-point correlation functions reads

$\displaystyle C_{2{\rm pt}}\left(\vec{p},T\right)$
$\displaystyle=\langle\mathcal{N}\left(-\vec{p},T\right)\overline{\mathcal{N}}\left(\vec{p},0\right)\rangle=\sum_{n}\frac{\left|\mathcal{Z}_{n}\left(\vec{p}\right)\right|^{2}}{2E_{n}\left(\vec{p}\right)}e^{-E_{n}\left(\vec{p}\right)T}$
(56) $\displaystyle C_{3{\rm
pt}}^{\left[\gamma^{\mu}\right]}\left(\vec{p}_{f},\vec{p}_{i},T;z,\tau\right)$
$\displaystyle=\langle\mathcal{N}\left(-\vec{p}_{f},T\right)\overline{\psi}\left(-\frac{z}{2},\tau\right)\gamma^{\mu}W^{(f)}\left(-\frac{z}{2},\frac{z}{2};A\right)\psi\left(\frac{z}{2},\tau\right)\overline{\mathcal{N}}\left(\vec{p}_{i},0\right)\rangle$
$\displaystyle=\sum_{n^{\prime},n}\frac{\mathcal{Z}_{n^{\prime}}\left(\vec{p}_{f}\right)\mathcal{Z}_{n}^{\dagger}\left(\vec{p}_{i}\right)}{4E_{n^{\prime}}\left(\vec{p}_{f}\right)E_{n}\left(\vec{p}_{i}\right)}\bra{n^{\prime}}\mathring{\mathcal{O}}^{\mu}\left(-\frac{z}{2},\frac{z}{2};A\right)\ket{n}e^{-E_{n^{\prime}}\left(\vec{p}_{f}\right)\left(T-\tau\right)}e^{-E_{n}\left(\vec{p}_{i}\right)\tau}\,,$
(57)

with the source and sink interpolating fields $\mathcal{N}$ separated by a
Euclidean time $T$, momentum-dependent interpolator-state overlaps given by
$\mathcal{Z}_{n}\left(\vec{p}_{i,f}\right)$. The bare quark bilinear,
abbreviated by $\mathring{\mathcal{O}}^{\mu}\left(z,\tau;A\right)$, is
introduced between the source and sink interpolators for Euclidean times
$1\leq\tau\leq T-1$. It is the ground state, or $n,n^{\prime}=0$ in Eqs. (57),
that corresponds to the physical nucleon state of interest.

Our computation of the correlation functions in Eqs. (56) and (57) using
distillation Peardon  _et al._ (2009) provides potential cost-saving benefits
when mapping $x$-dependent GPDs because distillation enables us to efficiently
calculate many combinations of source and sink momenta while reusing the most
expensive components of the calculation. The Wick contracted quark fields
smeared with distillation at both source and sink reduce Eqs. (56) and (57)
into traces over products of distinct, reusable computational units in the
distillation space. The independently generated hadron elementals and
perambulators encode the operator construction at source/sink and quark
propagation between time slices, respectively. The elementals and
perambulators are shared between both two-point and three-point correlation
functions (see Egerer  _et al._ (2021) Fig. 1). The color and spatial degrees
of freedom (a $N_{c}\times L^{3}$ dimensional space for quarks) are projected
into the lowest modes of the Laplacian (a $N_{\textrm{vec}}$ dimensional
space). Since hadronic physics is dominated by the long-distance modes of QCD,
which the lowest modes of the Laplacian are a proxy for, the number of modes
required to reproduce a ground state hadronic correlation function is much
smaller than triple the volume of modern lattices. However, correlators
computed in distillation effectively sample the whole source and sink
timeslices and therefore are equivalent to a large number of smeared point-
source propagators. As a result, our computation with four time sources on 348
configurations results in much more precise matrix elements than traditional
approaches.

A substantial amount of the computational resources for this calculation are
expended in computing generalized perambulators, or genprops for short:

$\Xi_{\alpha\beta;ab}^{{\mathbf{\Gamma}}(i,j)}\left(T_{f},T_{i};\tau,z\right)=\sum_{\vec{y}}\xi_{a}^{(i)\dagger}\left(T_{f}\right)D^{-1}_{\alpha\sigma;ac}\left(T_{f};\tau,\vec{y}\right)\mathbf{\Gamma}\left(\tau\right)D^{-1}_{\rho\beta;db}\left(\tau,\vec{x};T_{i}\right)\xi_{b}^{(j)}\left(T_{i}\right)\,,$
(58)

where color (Lorentz) indices are denoted by Latin (Greek) letters, $D^{-1}$
are inverses of the Wilson-clover fermion operator, and the Dirac structure of
the insertion is abbreviated
$\mathbf{\Gamma}=\gamma_{\sigma\rho}^{\mu}e^{-i\vec{q}\cdot\vec{y}}\delta_{\vec{y},\vec{x}+z\hat{z}}W_{cd}\left(\vec{x}+z\hat{z},\vec{x}\right)$.
The $\xi^{(k)}$ are eigenvectors of the three-dimensional gauge-covariant
Laplacian at a time slice $T$,

$-\nabla^{2}\left(T\right)\xi^{(k)}\left(T\right)=\lambda^{(k)}\left(T\right)\xi^{(k)}\left(T\right)\,,$
(59)

that defines the distillation operator. In practice, we compute inversions
against eigenvectors at both source and sink and use $\gamma_{5}$-hermiticity
to assemble each genprop
$\Xi_{\alpha\beta;ab}^{{\mathbf{\Gamma}}(i,j)}\left(T_{f},T_{i};\tau,z\right)$.
Within the distillation approach, we only need to produce genprops with all
desired Dirac matrices, Wilson line lengths, and $3$-momentum transfers
$\\{\vec{q}\\}$. We are then free to choose any off-forward matrix element,
defined by final (initial) momenta $\vec{p}_{f}$ ($\vec{p}_{i}$), that is
consistent with the pre-determined set of momentum transfers $\\{\vec{q}\\}$,
and compute that matrix element via comparatively cheap tensor contractions.

The distillation framework affords two particular advantages over other
methods generally employed. The first is that it enables a momentum projection
to be performed at each point in our two- and three-point functions. Thus it
provides a far better sampling of the gauge configurations. Second, the
factorization of correlation functions into the convolution of the
perambulators describing the quark propagation and the elementals representing
the interpolating operators admits the use of a variational method; we defer
the use of this for later work.

In this manner, distillation provides straight-forward access to what has been
called in the literature “symmetric” and “asymmetric” matrix elements
Bhattacharya  _et al._ (2022) without repeating nearly identical inversions
when varying the sources/sinks, which are otherwise required for standard
Gaussian-smeared correlation functions computed with the help of sequential
inversion methods Bernard  _et al._ (1986); Kilcup  _et al._ (1985).

###  3.1 Interpolator Construction

Mass eigenstates in the continuum are labeled according to angular momentum
and parity quantum numbers $J^{P}$, each of which corresponds to an
irreducible representation (irrep) of the continuous rotation group
$O\left(3\right)$. Fermionic and bosonic states are subsequently described by
double and single-valued representations, respectively. The projection of $J$
along a standard axis, typically taken to be $J_{z}$, provides a further
quantum number and is used to label the rows of the representation.

On a finite isotropic lattice the continuum rotational symmetry reduces to the
finite-dimensional cubic group or its dual, the octahedral group $O_{h}$.
Continuum mass eigenstates in this environment are characterized by their
patterns of subduction across the finite number of irreps $\Lambda$ of
$O_{h}$. As the nucleon is our numerical arena in this study, we will be
concerned with irreps of the double-cover octahedral group $O_{h}^{D}$. Due to
this implied many-to-one mapping, the spectrum of the ground-state nucleon
will be contaminated not only by excited states of the same continuum $J^{P}$
quantum numbers but by higher spin states as well.

For continuum states in motion, $J_{z}$ ceases to be a good quantum number and
is replaced by helicity $\lambda$; $J^{\lambda}$ then labels the irreps of the
little group, or subgroup, $O\left(3\right)$ that leaves the momentum vector
invariant. Unlike the continuum little group, which is independent of momentum
direction, the $O_{h}^{D}$ symmetry is broken further into little groups that
depend on ${}^{*}\left(\vec{p}\right)$, the group of rotations that leave
$\vec{p}$ invariant Moore and Fleming (2006), further compounding the mixing
of (continuum) mass eigenstates.

The reduced symmetry of an isotropic cubic lattice and the mixing of distinct
continuum $J^{P}$ states motivates our construction of nucleon interpolators
that transform irreducibly under the irreps of $O_{h}^{D}$ and its little
groups. Our strategy for isolating the $1/2^{+}$ ground-state nucleon with
distillation hinges on first building continuum interpolators
$\mathcal{O}^{J,P,M}(\vec{p}=\vec{0})$ that possess definite $J^{P}$ and
flavor quantum numbers Edwards  _et al._ (2011); Dudek and Edwards (2012).
This starting point reads symbolically

$\bra{\vec{p}=\vec{0};J^{\prime},P^{\prime},M^{\prime}}\left[\mathcal{O}^{J,P,M}(\vec{p}=\vec{0})\right]^{\dagger}\ket{\vphantom{\vec{p}=\vec{0};J^{\prime},P^{\prime},M^{\prime}}0}=Z_{J}\delta_{JJ^{\prime}}\delta_{PP^{\prime}}\delta_{MM^{\prime}}\,,$
(60)

where the overlap of a continuum spin-$J^{\prime}$ state with the continuum
interpolator $\mathcal{O}^{J,P,M}$ is denoted by $Z_{J}$. To reliably isolate
nucleons with a momentum $\vec{p}$ within our isotropic lattice, which breaks
$O_{h}^{D}$ into its little groups, we implement an algorithm Thomas  _et al._
(2012) first envisioned for mesons in flight. The starting point is continuum
helicity operators, which are obtained by boosting our continuum operators at
rest by an amount $\left|\vec{p}\right|$ along their quantization axis
$\hat{z}$, followed by a series of active rotations that enforce a change of
basis from $J_{z}$ to $\lambda$. This process reads symbolically as:

$\left[\mathbb{O}^{J,P,\lambda}\left(\vec{p}\right)\right]^{\dagger}=\sum_{M}\mathcal{D}_{M\lambda}^{(J)}\left(R\right)\left[\mathcal{O}^{J,P,M}\left(\left|\vec{p}\right|\hat{z}\right)\right]^{\dagger}\,,$
(61)

where a continuum helicity operator
$\mathbb{O}^{J,P,\lambda}\left(\vec{p}\right)$ is obtained via application of
a Wigner-$\mathcal{D}$ matrix $\mathcal{D}_{M\lambda}^{(J)}\left(R\right)$
onto the boosted $J_{z}$-quantized operator
$\mathcal{O}^{J,P,M}\left(\left|\vec{p}\right|\hat{z}\right)$. Ref. Thomas
_et al._ (2012) highlights the ambiguity within a finite cubic volume in
rotating from $\left|\vec{p}\right|\hat{z}$ to $\vec{p}$ in this manner. A
consistent convention is established by splitting the rotation $R$ into two
distinct rotations,

$R=R_{\rm lat}R_{\rm ref}\,,$ (62)

where $R_{\rm ref}$ rotates $\left|\vec{p}\right|\hat{z}$ to a reference
vector $\vec{p}_{\rm ref}$222Each $O_{h}^{D}$ little group defined by
$*\left(\vec{p}\right)$ is assigned a reference direction such that momenta
within the same little group (e.g. $\left(0,1,0\right)$ and
$\left(-2,0,0\right)$) can be treated consistently. For example, we assign
momenta within the ${\rm Dic}_{2}$ little group, such as $\left(1,0,1\right)$,
the reference direction $\vec{p}_{\rm ref}=\left(0,1,1\right)$. Euler angles
within the $zyz$-convention then rotate $\left|\vec{p}\right|\hat{z}$ to
$\vec{p}_{\rm ref}$., and $R_{\rm lat}$ is a lattice rotation that rotates
$\vec{p}_{\rm ref}$ to $\vec{p}$. In general, $R_{\rm ref}$ may not be allowed
within a cubic volume, however, its action is permissible given that it merely
produces the aforementioned basis change. We have validated this procedure for
all possible momenta within a finite cubic volume, with the exception of
momenta of the form $\vec{p}=\left(n,m,l\right)$, with distinct
$n,m,l\in\mathbb{Z}\backslash\\{0\\}$, which are left invariant under the
$C_{2}$ little group of $O_{h}^{D}$. Hence, we will exclude $C_{2}$-type
momenta from this study. The resulting continuum helicity operators respect
the following overlap criterion:

$\bra{\vec{p};J^{\prime},P^{\prime},\lambda^{\prime}\vphantom{\left[\mathbb{O}^{J,P,\lambda}\left(\vec{p}\right)\right]^{\dagger}}}\left[\mathbb{O}^{J,P,\lambda}\left(\vec{p}\right)\right]^{\dagger}\ket{\vphantom{\vec{p};J^{\prime
P^{\prime}},\lambda^{\prime}}\\!\
0}=Z_{J^{\prime}J;P^{\prime}P;\lambda}\delta_{\lambda\lambda^{\prime}}\,.$
(63)

In motion, parity $P$ and spin $J$ are no longer conserved quantities, because
their symmetries are broken by the explicit direction $\vec{p}$, while
helicity remains conserved in flight. The construction of our interpolators,
or equivalently subduced helicity operators, is finalized by specifying how a
continuum helicity $\lambda$ subduces into a little group irrep $\Lambda$ with
rows $\mu\in\\{1,\cdots,{\rm dim}\left(\Lambda\right)\\}$:

$\left[\mathbb{O}^{[J,P,\left|\lambda\right|]}_{\Lambda,\mu}\left(\vec{p}\right)\right]^{\dagger}=\sum_{\hat{\lambda}=\pm\lambda}\left[S_{\Lambda,\mu}^{J,\hat{\lambda}}\right]^{\ast}\left[\mathbb{O}^{J,P,\hat{\lambda}}\left(\vec{p}\right)\right]^{\dagger}\,,$
(64)

where the subduced helicity operators used in our study are linear
combinations of each value of continuum helicity, the precise combinations of
which are set by so-called subduction coefficients
$S_{\Lambda,\mu}^{J,\hat{\lambda}}$. As a consequence, the continuum energy
eigenstates as part of the spectral content of the three-point functions in
Eq. (57) are instead subduced eigenstates of helicity.

Interpolators constructed in this manner are assigned a compact spectroscopic
notation $X^{2S+1}L_{\mathcal{P}}J^{P}$, with $X$ denoting the hadronic state,
$S$ the overall Dirac spin, $L$ the angular momentum induced by covariant
derivatives with permutation symmetry $\mathcal{P}$, and $J^{P}$ the continuum
spin/parity quantum numbers. We select the simplest non-relativistic and
spatially-local operator, denoted $N^{2}S_{s}\frac{1}{2}^{+}$, to interpolate
the ground-state nucleon from the vacuum. A follow-on study could make use of
an expanded operator basis (such as in Refs. Khan  _et al._ (2021); Egerer
_et al._ (2019, 2022a)) to improve excited state control within each
$\vec{p}_{i}/\vec{p}_{f}$ momentum channel.

###  3.2 Computing the kinematic matrix

The connection between the physically relevant matrix elements and the
subduced matrix elements we extract from the three-point correlation functions
was established in Edwards  _et al._ (2023) for the isovector quark helicity
PDFs of the nucleon. In the interest of self-containment, we highlight the
points that are key to the extraction of the GPDs.

To isolate the Lorentz invariant amplitudes, the subduced matrix elements can
be related to the helicity matrix elements, which requires inverting the
effect of the subduction coefficients $S^{J,\hat{\lambda}}_{\Lambda,\mu}$, Eq.
(64), that describe the way the different helicity $\lambda$ states are
related (subduced) into the various rows $\mu\in\\{1,{\rm
dim}\left(\Lambda\right)\\}$ of the irreps $\Lambda$. Following the
conventions of Thomas  _et al._ (2012), the creation operators require the
complex conjugate of the coefficient. After, the Lorentz decomposition can be
applied to the continuum matrix element. Finally, the subduced matrix elements
are given by a linear combination of invariant amplitudes:

$\bra{\vec{p},\Lambda,\mu_{f}}\mathcal{J}^{\Lambda_{\Gamma},\mu_{\Gamma}}\ket{\vec{p},\Lambda,\mu_{i}}=\sum_{k}\sum_{\lambda_{f},\lambda_{\Gamma},\lambda_{i}}S^{\Lambda,J}_{\mu_{f},\lambda_{f}}\left[S_{\mu_{i},\lambda_{i}}^{\Lambda,J}\right]^{*}\mathcal{K}_{k}\left(\lambda_{f}\left[J,\vec{p}\right];\lambda_{i}\left[J,\vec{p}\right]\right)\mathcal{A}_{k}\left(\nu,z^{2}\right)\,.$
(65)

Here each term is weighted by a helicity-dependent kinematic prefactor
$\mathcal{K}_{k}\left(\lambda_{f}\left[J,\vec{p}\right];\lambda_{i}\left[J,\vec{p}\right]\right)$,
given as the Lorentz covariant coefficients of $\mathcal{A}_{k}$ in Eq. (10),
whose irreducibility under the relevant cubic (sub)group is encoded via
subduction coefficients.

Following the procedure in Edwards  _et al._ (2023), one need not relate the
subduced matrix element to the helicity matrix element as an intermediate
step. Instead one can calculate the kinematic matrices,

$K_{k}(\mu_{f},[\Lambda,\vec{p}])=\sum_{\lambda_{f},\lambda_{\Gamma},\lambda_{i}}S^{\Lambda,J}_{\mu_{f},\lambda_{f}}\left[S_{\mu_{i},\lambda_{i}}^{\Lambda,J}\right]^{*}\mathcal{K}_{k}\left(\lambda_{f}\left[J,\vec{p}\right];\lambda_{i}\left[J,\vec{p}\right]\right)\,,$
(66)

and invert the relationship to obtain the amplitudes directly. As discussed in
Edwards  _et al._ (2023), this viewpoint exposes the ever-present potential
issues with mixing between states of different parity or continuum spin, all
of which belong to the same finite volume irrep $\Lambda$. In QCD, the
spectrum is ordered in spin, and ultimately the ground state which will be the
desired is the $\frac{1}{2}^{+}$ nucleon. These other states contribute to
excited state contamination that can be modeled and removed. As such, the
kinematic matrices will use only spin-$1/2$ spinors.

The construction of the subduced spinors was outlined in detail in Edwards
_et al._ (2023). The algorithm begins with constructing a continuum spinor
with spin and momentum in the $z$ direction. This spinor is then rotated, as
the interpolating field was in Eq. (61), to align with the correct lattice-
allowed momenta. Finally, the spinor is subduced into the finite volume
$\Lambda$ irrep as the interpolator was in Eq. (64). The convenience of
parameterizing a rotation with Euler angles is that these rotations are
independent of the spin representation of the object being considered. Armed
with subduced spinors, we can construct the Lorentz covariant coefficients of
each of the desired invariant amplitudes in their relation to the subduced
matrix elements.

###  3.3 From the matrix element to the amplitudes

To isolate each amplitude individually we need to consider many combinations
of initial/final rows as well as $\Gamma$ structures. We expand upon the
strategy of Edwards  _et al._ (2023) to address this challenge.

In the case of the unpolarized quark forward matrix element, the nucleon
matrix element $M^{0}\left(p,p,z\right)$ was sufficient to isolate
$\mathcal{A}_{1}(\nu,\xi=0,t,z^{2})$ and obtain a PDF. With so many amplitudes
in the off-forward case, Eq. (10), this is impossible, just as it was in the
case of helicity quark and gluon matrix elements Edwards  _et al._ (2023);
Balitsky  _et al._ (2022); Egerer  _et al._ (2022a). While fixing the momenta
and separation, the matrix element is considered in the four possible initial
and final spin combinations333Regardless of source or sink momentum, each
cubic irrep relevant for the (continuum) $J^{P}=1/2^{+}$ nucleon is two-
dimensional. and for each of the four $\gamma^{\mu}$ matrices. These combined
sixteen matrix elements can be used to isolate the eight unknown
$\mathcal{A}_{k}$ given the matrix relation

$\underset{16\times 1}{\mathbf{M}}=\underset{16\times
8}{\mathbf{K}}\times\underset{8\times 1}{\mathbf{A}}\,,$ (67)

where we define a column vector of fitted subduced matrix elements by

$\underset{16\times
1}{\mathbf{M}}=\begin{bmatrix}\begin{pmatrix}M^{0}\left(p_{f},p_{i},z\right)_{11}\\\
M^{0}\left(p_{f},p_{i},z\right)_{12}\\\
M^{0}\left(p_{f},p_{i},z\right)_{21}\\\
M^{0}\left(p_{f},p_{i},z\right)_{22}\end{pmatrix}^{T}\begin{pmatrix}M^{1}\left(p_{f},p_{i},z\right)_{11}\\\
M^{1}\left(p_{f},p_{i},z\right)_{12}\\\
M^{1}\left(p_{f},p_{i},z\right)_{21}\\\
M^{1}\left(p_{f},p_{i},z\right)_{22}\end{pmatrix}^{T}\begin{pmatrix}M^{2}\left(p_{f},p_{i},z\right)_{11}\\\
M^{2}\left(p_{f},p_{i},z\right)_{12}\\\
M^{2}\left(p_{f},p_{i},z\right)_{21}\\\
M^{2}\left(p_{f},p_{i},z\right)_{22}\end{pmatrix}^{T}\begin{pmatrix}M^{3}\left(p_{f},p_{i},z\right)_{11}\\\
M^{3}\left(p_{f},p_{i},z\right)_{12}\\\
M^{3}\left(p_{f},p_{i},z\right)_{21}\\\
M^{3}\left(p_{f},p_{i},z\right)_{22}\end{pmatrix}^{T}\end{bmatrix}^{T}\,.$
(68)

The kinematic matrix is built from each Lorentz covariant structure in Eq.
(10) sandwiched between initial/final-state subduced spinors:

$\underset{16\times
8}{\mathbf{K}}=\begin{bmatrix}\langle\langle\gamma^{0}\rangle\rangle_{11}&z^{0}\langle\langle\mathds{1}\rangle\rangle_{11}&i\langle\langle\sigma^{0z}\rangle\rangle_{11}&\cdots&\frac{i}{2m}z^{0}\langle\langle\sigma^{z\Delta}\rangle\rangle_{11}\\\
{\LARGE{\color[rgb]{1,0,0}\vdots}}&{\LARGE{\color[rgb]{1,0,0}\vdots}}&{\LARGE{\color[rgb]{1,0,0}\vdots}}&&{\LARGE{\color[rgb]{1,0,0}\vdots}}\\\
\langle\langle\gamma^{1}\rangle\rangle_{11}&z^{1}\langle\langle\mathds{1}\rangle\rangle_{11}&i\langle\langle\sigma^{1z}\rangle\rangle_{11}&\cdots&\frac{i}{2m}z^{1}\langle\langle\sigma^{z\Delta}\rangle\rangle_{11}\\\
{\LARGE{\color[rgb]{1,0,0}\vdots}}&{\LARGE{\color[rgb]{1,0,0}\vdots}}&{\LARGE{\color[rgb]{1,0,0}\vdots}}&&{\LARGE{\color[rgb]{1,0,0}\vdots}}\\\
\langle\langle\gamma^{2}\rangle\rangle_{11}&z^{2}\langle\langle\mathds{1}\rangle\rangle_{11}&i\langle\langle\sigma^{2z}\rangle\rangle_{11}&\cdots&\frac{i}{2m}z^{2}\langle\langle\sigma^{z\Delta}\rangle\rangle_{11}\\\
{\LARGE{\color[rgb]{1,0,0}\vdots}}&{\LARGE{\color[rgb]{1,0,0}\vdots}}&{\LARGE{\color[rgb]{1,0,0}\vdots}}&&{\LARGE{\color[rgb]{1,0,0}\vdots}}\\\
\langle\langle\gamma^{3}\rangle\rangle_{11}&z^{3}\langle\langle\mathds{1}\rangle\rangle_{11}&i\langle\langle\sigma^{3z}\rangle\rangle_{11}&\cdots&\frac{i}{2m}z^{3}\langle\langle\sigma^{z\Delta}\rangle\rangle_{11}\\\
{\LARGE{\color[rgb]{1,0,0}\vdots}}&{\LARGE{\color[rgb]{1,0,0}\vdots}}&{\LARGE{\color[rgb]{1,0,0}\vdots}}&&{\LARGE{\color[rgb]{1,0,0}\vdots}}\\\
\end{bmatrix}\,,$ (69)

and the vector

$\displaystyle\underset{8\times
1}{\mathbf{A}}=\begin{bmatrix}\mathcal{A}_{1}&\mathcal{A}_{2}&\mathcal{A}_{3}&\mathcal{A}_{4}&\mathcal{A}_{5}&\mathcal{A}_{6}&\mathcal{A}_{7}&\mathcal{A}_{8}\end{bmatrix}^{T}$
(70)

contains the required amplitudes with omitted $(\nu,\xi,t,z^{2})$ arguments.
The subscripts in each entry of Eqs. (68) and (69) are integer pairs
designating the sink/source irrep row combination
$\left(\mu_{f}\mu_{i}\right)$; red dots are employed for brevity to indicate
repeated Lorentz covariant structures sandwiched between distinct
$\left(\mu_{f}\mu_{i}\right)$ combinations.

Given the linear system of sixteen equations for eight amplitudes (67), our
goal is to obtain the most probable $\mathcal{A}_{k}$. The traditional
solution to an over-constrained linear system of equations is to find the
minimum of $\chi^{2}=||\mathbf{M}-\mathbf{K}\mathbf{A}||_{2}^{2}$ where
$||.||_{2}$ denotes the $l^{2}$ norm. The solution is given by
$\bar{\mathbf{A}}=\mathbf{K}^{+}\mathbf{M}$, where $\mathbf{K}^{+}$ is the
pseudoinverse of $\mathbf{K}$ Penrose (1956). The pseudoinverse can be
calculated via singular value decomposition (SVD). Numerically, we use a
cutoff on the relative size of the smallest to the largest singular value of
$\epsilon=10^{-15}$ in the `linalg.pinv` function of the publicly available
`NumPy` package. This cutoff formally shifts the matrix from the original
definition of the Moore-Penrose pseudoinverse but increases the numerical
stability.

Half of the equations in Eq. (67) are related to the other half by various
symmetries. Ultimately the four helicity combinations only give two
independent equations. For most kinematics, this is enough to reliably
reconstruct the eight amplitudes $\mathcal{A}_{k}$. A few specific situations
arise however:

  1. 1.

If $z=0$, the kinematic factors of all amplitudes except $\mathcal{A}_{1,4,5}$
cancel. Those are therefore the only amplitudes that we can hope to extract.
We will discuss the results at $z=0$ specifically in section 5.

  2. 2.

If the spatial momentum transfer $\vec{q}=\vec{p}_{f}-\vec{p}_{i}$ is exactly
along the $z$ axis, several issues may appear. For some kinematic pairs and
$z\neq 0$ (respectively $z=0$), there are fewer than eight (respectively
three) singular values in the kinematic matrix, meaning that there exist exact
(up to $\epsilon=10^{-15}$) linear relations between some amplitudes that
cannot be separated from one another. In other cases, in spite of eight
different singular values, the unitary matrices entering the SVD present large
similarities in their columns and result in a very bad case of uncertainty
propagation. As we observe in a consistent fashion that these kinematics
behave poorly, we remove them altogether. They are not included in the summary
plot of Fig. 1.

If $z\neq 0$ and $\vec{q}$ is not along the $z$ axis ($q_{1}$ or $q_{2}\neq
0$), then we observe that, due to our choice of $z^{\mu}=(0,0,0,z_{3})$, only
the $M^{3}(p_{f},p_{i},z)$ matrix elements are sensitive to $\mathcal{A}_{2}$
and $\mathcal{A}_{8}$. In practice, this means that the $M^{3}$ matrix
elements are entirely and exclusively used to characterize those two
amplitudes. Although it has been argued that using $M^{3}$ introduces a finite
mixing Constantinou and Panagopoulos (2017), using it does not change the
extraction of any other amplitude, and therefore does not impact our light-
cone GPDs, which do not involve either $\mathcal{A}_{2}$ or $\mathcal{A}_{8}$.

When $z=0$, the $\mathcal{A}_{2}$ and $\mathcal{A}_{8}$ amplitudes are absent
from the decomposition and the $M^{3}$ matrix elements have a constraining
power on the amplitudes $\mathcal{A}_{1,4,5}$. We observe, however, that the
difference is typically of the order of a few percent at most on
$\mathcal{A}_{1,4}$, which contain the elastic form factors (EFFs). As we will
see later, we perform the extraction of the EFFs using both the $z=0$ and
$z\neq 0$ data, and obtain very similar results up to a $\mathcal{O}(z^{2})$
contamination which cannot be attributed to the $M^{3}$ data (since it plays
no role when $z\neq 0$). We give more details on the question of the use of
$M^{3}$ in appendix B.

As an extension to this approach, one could consider the effect of additional
momentum combinations that select the same values of $\nu$, $\xi$, and $t$ to
produce an even more constrained system of equations. This could especially
help in reducing lattice systematic effects, such as discretization errors of
$\mathcal{O}\left(aq\right)$. Using the relationships in section 2.4 between
amplitudes of $\pm\nu$ and $\pm\xi$, even more combinations can be included.
In this paper, we limit ourselves to averaging the $\pm z_{3}$ separations and
perform an external amplitude extraction for each pair
$(\vec{p}_{f},\vec{p}_{i})$. Nonetheless, we have a certain redundancy,
because we have data for $\pm\xi$ and various $\nu$ at the same or very nearby
$t$, which will enables us to reduce some of the lattice artifacts.

An alternative possibility to perform the extraction is to remove or average
directly the trivial duplicate constraints and form a true matrix inverse. The
solution of this linear inverse would exactly reproduce the analog of the
analytic formulas presented in Bhattacharya  _et al._ (2022) to extract the
amplitudes from the matrix elements. In our calculation, the spinors required
to calculate $\langle\langle\Gamma\rangle\rangle_{\mu_{f}\mu_{i}}$ are
complicated functions of the momenta, depending upon its magnitude, the
specific finite volume irreducible representation, and choices of convention.
With many momenta combinations, solving the analytic formulas would require a
numerical algebra implementation for practical use. The approach outlined in
this section replaces the step of re-deriving these analytic solutions for
varying considerations such as off-axis separations or new momentum frames
with a generic linear algebra problem.

Finally, we have cross-checked the stability of our SVD extraction by using an
$l^{1}$-norm minimization of $||\mathbf{M}-\mathbf{K}\mathbf{A}||_{1}$, also
known as least absolute deviation (LAD) estimator LAD (2008). This procedure
has disadvantages compared to SVD because it requires an iterative
minimization for each jack-knife sample of the data, which is hundreds of
times computationally more expensive than a simple matrix operation. The
results are consistent, as one can check in appendix A, where we present the
extraction of a full set of amplitudes in some kinematics.

##  4 Analysis of correlation functions

We isolate the desired bare matrix element
$M^{\mu}\left(\vec{p}_{f},\vec{p}_{i},z\right)$ using two strategies to
evaluate possible systematics linked to excited state contamination. Our first
procedure is the summation method Maiani  _et al._ (1987); Capitani  _et al._
(2012). We form an optimized ratio of three-point and two-point correlation
functions,

$R^{\mu}\left(\vec{p}_{f},\vec{p}_{i},z;T,\tau\right)=\frac{C^{\mu}_{3{\rm
pt}}\left(\vec{p}_{f},\vec{p}_{i},z;T,\tau\right)}{C_{2{\rm
pt}}\left(\vec{p}_{f};T\right)}\sqrt{\frac{C_{2{\rm
pt}}\left(\vec{p}_{i};T-\tau\right)C_{2{\rm
pt}}\left(\vec{p}_{f};\tau\right)C_{2{\rm
pt}}\left(\vec{p}_{f};T\right)}{C_{2{\rm
pt}}\left(\vec{p}_{f};T-\tau\right)C_{2{\rm
pt}}\left(\vec{p}_{i};\tau\right)C_{2{\rm pt}}\left(\vec{p}_{i};T\right)}}\,,$
(71)

where $T$ is the temporal separation between the source/sink interpolating
fields and $0\leq\tau\leq T$ is the current insertion time slice. The three-
point correlation functions are computed for the values of
$T\in\\{4,6,8,10,12,14\\}$ and every $1\leq\tau\leq T-1$. For asymptotically
large temporal separations ($0\ll\tau\ll T$), the ratio
$R^{\mu}\left(\vec{p}_{f},\vec{p}_{i},z;T,\tau\right)$ will become
proportional to the desired matrix element,

$R^{\mu}\left(\vec{p}_{f},\vec{p}_{i},z;T,\tau\right)\overset{\tau,T\rightarrow\infty}{\longrightarrow}\frac{1}{\sqrt{2}}\times\frac{1}{\sqrt{4E_{0}\left(\vec{p}_{f}\right)E_{0}\left(\vec{p}_{i}\right)}}M^{\mu}\left(p_{f},p_{i},z\right)\,.$
(72)

The additional factor of $\sqrt{2}$ is due to the normalization of the
isovector vector current.444 The electromagnetic current in the light-quark
sector reads
$\mathcal{J}^{\mu}=\frac{2}{3}\overline{u}\gamma^{\mu}u-\frac{1}{3}\overline{d}\gamma^{\mu}d$,
which can be expressed in terms of isovector $\rho^{\mu}$ and isoscalar
$\omega^{\mu}$ components according to
$\mathcal{J}^{\mu}=\frac{\rho^{\mu}}{\sqrt{2}}+\frac{\omega^{\mu}}{3\sqrt{2}}$,
where
$\rho^{\mu}=\left(\overline{u}\gamma^{\mu}u-\overline{d}\gamma^{\mu}d\right)/\sqrt{2}$
and
$\omega^{\mu}=\left(\overline{u}\gamma^{\mu}u+\overline{d}\gamma^{\mu}d\right)/\sqrt{2}$.
We use the current $\rho^{\mu}$ to probe the isovector flavor structure of the
off-forward nucleon matrix elements. We form
$R^{\mu}\left(\vec{p}_{f},\vec{p}_{i},z;T,\tau\right)$ by treating both three-
point and all two-point correlators as complex-valued functions, as the
construction of the subduced helicity operators we use as interpolators,
detailed in section 3.1, involves a delicate interplay of complex phases. A
depiction of the ratio for a specific choice of $\vec{p}_{f}$ and
$\vec{p}_{i}$ is given in Fig. 5.

Then the matrix elements are extracted by summing the ratio
$R^{\mu}\left(\vec{p}_{f},\vec{p}_{i},z;T,\tau\right)$ over the insertion time
slice $\tau$, excluding all time separations below a threshold $t_{s}$,

$\mathcal{R}_{t_{s}}^{\mu}\left(\vec{p}_{f},\vec{p}_{i},z;T\right)\equiv\sum_{\tau/a=t_{s}}^{T-t_{s}}R^{\mu}\left(\vec{p}_{f},\vec{p}_{i},z;T,\tau\right)\,.$
(73)

The summed ratio
$\mathcal{R}_{t_{s}}^{\mu}\left(\vec{p}_{f},\vec{p}_{i},z;T\right)$ is a
geometric series that depends linearly on the bare matrix element
$M^{\mu}\left(p_{f},p_{i},z\right)$. We apply the functional form

$\mathcal{R}^{\mu}_{\rm
fit}\left(\vec{p}_{f},\vec{p}_{i},z;T\right)=C+M^{\mu}\left(p_{f},p_{i},z\right)T$
(74)

to extract the bare matrix element. This procedure should yield errors of the
order $\mathcal{O}\left(e^{-\Delta ET}\right)$, where $\Delta E$ is the energy
gap between the lowest-lying effective excited state and the ground state. We
show in the right panel of Fig. 6 the impact of the exclusion of short time
separations, both in the ratio method (green points) and in the exponential
fits that we will describe now.

Figure 5: Ratio $R^{0}(\vec{p}_{f},\vec{p}_{i},z=a;T,\tau)$ for
$\vec{p}_{f}=(2,1,0)$ and $\vec{p}_{i}=(0,1,2)$ (lattice units) with proton
rows $1,1$ for various values of $T$ and $\tau$. We show the result of the
matrix element fit using either the summation method or an exponential fit
with one excited state, with the exclusion of all $\tau$ and $T-\tau<3$. The
data at $T=14$ is dotted for legibility purposes.

Figure 6: (left) Dispersion relation $E_{0}\left(|\vec{p}|\right)$ obtained
from non-linear fits of the ground state and first excited state on the two-
point correlation functions. The dotted line is the continuum dispersion
relation. – (right) Impact of the exclusion of short-time separations below
$t_{s}$ on the extraction of the same bare matrix element as in Fig. 5. We
show the result of a ground state exponential fit (blue), exponential fit with
one excited state (red), and ratio method (green). The standard result in this
analysis will be using both the red point at $t_{s}=3$ and at $t_{s}=4$ to
evaluate excited state uncertainty.

Our second procedure extracts the bare matrix elements using exponential fits,
following the spectral decompositions in Eqs. (56) and (57).

We first perform a non-linear fit of the ground state and first excited state
on the two-point correlation functions, excluding small source-sink
separations. For each value of $|\vec{p}|$, we use several rotationally
equivalent two-point functions to increase the robustness of our fit. As a
result, we obtain for each jack-knife sample the fitted values of
$\mathcal{Z}_{n}\left(\vec{p}\right)$ and $E_{n}\left(\vec{p}\right)$ for
$n=0$ and $1$. The dispersion relation of the ground state is shown on the
left panel of Fig. 6.

Then we use the fitted values of $\mathcal{Z}_{n}\left(\vec{p}\right)$ and
$E_{n}\left(\vec{p}\right)$ in the decomposition of the three-point
correlation function. We first fit the three-point function assuming pure
ground state dominance, that is, using only the $n=0$ fits of the two-point
functions. This gives the blue points in the right panel of Fig. 6, where we
sequentially exclude the small separations $\tau$ and $T-\tau<t_{s}$. It is
clear that the points follow a non-constant trend, a sign of sizeable excited
state contamination. Based on these results, we use both the results for $n=0$
and $n=1$ (red points in the right panel of Fig. 6), meaning that we fit the
three-point function with four free parameters corresponding, in principle, to
$\langle 0|\mathcal{O}|0\rangle$, $\langle 0|\mathcal{O}|1\rangle$, $\langle
1|\mathcal{O}|0\rangle$ and $\langle 1|\mathcal{O}|1\rangle$. The results are
now generally in very good agreement with the ratio fits. In general, we
observe that the exponential fit with one excited state seems less prone to
fluctuations compared to the ratio method when increasing the cut in $t_{s}$.
Therefore, we will use the excited state exponential fit with a cut $t_{s}=3$
and $t_{s}=4$ to include a measure of excited state contamination.

Finally, we form the ratio with the matrix element for
$\vec{p}_{f}=\vec{p}_{i}=0$, $\Lambda_{f}=\Lambda_{i}=1$ and $\mu=0$ to form
the RGI reduced matrix element, Eq. (7).

The process of subduction detailed in section 3.1 leads to ${\rm
dim}\left(\Lambda_{f}\right)\times{\rm dim}\left(\Lambda_{i}\right)$
determinations of each matrix element
$M^{\mu}\left(\vec{p}_{f},\vec{p}_{i},z\right)$, where $\Lambda_{f}$ and
$\Lambda_{i}$ are the irreps at sink and source, respectively. The fitted
matrix elements will be determined for each combination of initial and final
interpolator rows $\mu_{i}/\mu_{f}$.

##  5 Elastic form factors

###  5.1 Local matrix elements

Let us first consider the local limit $z=0$. Then, Eq. (10) reduces to a
decomposition resembling the standard form factors parametrizing an
electromagnetic interaction with a spin-1/2 nucleon:

$\displaystyle M^{\mu}\left(p_{f},p_{i},0\right)$
$\displaystyle=\langle\langle\gamma^{\mu}\rangle\rangle\mathcal{A}_{1}\left(0,\xi,t,0\right)+\frac{i}{2m}\langle\langle\sigma^{\mu
q}\rangle\rangle\mathcal{A}_{4}\left(0,\xi,t,0\right)+\frac{q^{\mu}}{2m}\langle\langle\mathds{1}\rangle\rangle\mathcal{A}_{5}\left(0,\xi,t,0\right),$
(75) $\displaystyle=\langle\langle\gamma^{\mu}\rangle\rangle
F_{1}\left(t\right)+\frac{i}{2m}\langle\langle\sigma^{\mu q}\rangle\rangle
F_{2}\left(t\right)+\frac{q^{\mu}}{2m}\langle\langle\mathds{1}\rangle\rangle\mathcal{A}_{5}\left(0,\xi,t,0\right)\,,$
(76)

where $F_{1}\left(t\right)$ and $F_{2}\left(t\right)$ are the familiar Dirac
and Pauli form factors of the nucleon. The Ward identity $q_{\mu}M^{\mu}=0$
requires that $\mathcal{A}_{5}$ vanishes in the local limit.555The terms
$q_{\mu}\langle\langle\gamma^{\mu}\rangle\rangle$ and
$q_{\mu}\langle\langle\sigma^{\mu q}\rangle\rangle$ evaluate to zero via the
Dirac equation, while the term
$q_{\mu}q^{\mu}\langle\langle\mathds{1}\rangle\rangle$ does not vanish. Any
departure of $\mathcal{A}_{5}\left(\nu,\xi,t,z^{2}\right)$ from zero in the
local limit can thus be interpreted as the degree to which lattice artifacts
systematically affect a calculation of Eq. (10) since we use the local vector
current which is not the conserved vector current on the lattice.

Figure 7: Elastic form factors $F_{1}(t)$, $F_{2}(t)$ and $\mathcal{A}_{5}(t)$
extracted from the local matrix elements using all skewness data. The
$\mathcal{A}_{5}(t)$ data are plotted with a color scheme that reflects the
skewness at each point. The 15 grey areas represent the 15 bins in $t$ which
will come into play later on.

We present in Fig. 7 the results of the elastic form factors on the 186 pairs
of momenta $(\vec{p}_{f},\vec{p}_{i})$ whose kinematic coverage is displayed
on Fig. 1. We have merged the errorbars provided by the cut of lattice time
separation $t_{s}=3$ and 4 (see section 4). As expected, $F_{1}(t)$ and
$F_{2}(t)$ are globally independent of $\xi$. The amplitude
$\mathcal{A}_{5}(0,\xi,t,0)$ is generally at most a few sigmas away from zero
but there seems to be a non-vanishing signal, loosely odd in $\xi$ as expected
using the discrete symmetries of section 2.4. The equivalent term in
Bhattacharya  _et al._ (2022) was found to be small, but nonetheless non-zero
for $z=0$. App. B describes the effect on this systematic error from the
inclusion or exclusion of $M^{3}$ data.

When using the extraction of bare matrix elements from the summation of ratios
instead of the exponential fit, an even starker signal of non-vanishing
$\mathcal{A}_{5}(t)$ appears. The fact that $\mathcal{A}_{5}(t)$ at $z=0$
presents an enhanced sensitivity to the treatment of excited states suggests
that some underestimation of excited state contamination may play a role in
addition to the discretization errors in the non-vanishing value of
$\mathcal{A}_{5}(t)$. On the other hand, $F_{1}(t)$ and $F_{2}(t)$ seem
generally unaltered by the behavior of $\mathcal{A}_{5}(t)$ and its skewness
dependence. In fact, performing the SVD while ignoring altogether the
$\mathcal{A}_{5}(t)$ contribution does not produce a significant change in the
extraction of $F_{1}(t)$ and $F_{2}(t)$.

A recent expansive calculation of nucleon form factors by the Nucleon Matrix
Elements Collaboration Park  _et al._ (2022) included determinations of the
electric and magnetic Sachs form factors

$G_{E}\left(t\right)=F_{1}\left(t\right)+\frac{t}{4m^{2}}F_{2}\left(t\right)\,,\qquad
G_{M}\left(t\right)=F_{1}\left(t\right)+F_{2}\left(t\right)\,,$ (77)

on a Wilson clover fermion ensemble characterized by a pion mass ($m_{\pi}\sim
270$ MeV) that is slightly lighter than the a094m358 ensemble we consider in
this work. The result is shown in Fig. 8. As expected, the extraction with a
lighter pion mass has a more pronounced decrease in $t$.

Figure 8: Comparison between our extraction of the $G_{E}$ and $G_{M}$ form
factors from local matrix elements and the determination in Park  _et al._
(2022) on a slightly lighter Wilson clover fermion ensemble (pion mass of 270
Mev vs. 358 MeV in this work).

The next step of our analysis consists in forming bins of data in $t$, as
indicated by the grey regions in Fig. 7. This is beneficial for several
reasons. First, by averaging over data with various
$(\vec{p}_{f},\vec{p}_{i})$ in a single bin, we hope to mitigate some lattice
systematic uncertainty. For instance, we have many cases of symmetry
$\xi\longleftrightarrow-\xi$ in our data, but also different magnitudes of
$(\vec{p}_{f},\vec{p}_{i})$ in the same bin producing different excited state
behavior. Furthermore, the binning results in a decrease in the number of
degrees of freedom for the final fit of the $t$-dependence, meaning that the
empirical covariance matrix is a better estimate of the true covariance of the
data. Indeed, to estimate the covariance matrix, one should typically have
many more measurements than random variables. Here, with 186 kinematic values
for 348 gauge configurations, the stability of the empirical covariance matrix
is reduced. Trimming down the 186 kinematic values to fifteen bins improves
the confidence in this estimation. Finally, the binning gives us access to a
significant number of values of $(\nu,\xi\nu)$ in each bin when $z\neq 0$,
which will allow for the extraction of generalized form factors in the next
section.

We perform the binning in two steps. First, we fit each amplitude by a
constant within each bin. Then we take the fifteen data points obtained in
this fashion, and fit an overall $t$-dependence using a dipole fit,

$F(t)=A\left(1-\frac{t}{\Lambda^{2}}\right)^{-2}\,.$ (78)

8 Next, we displace the points in each bin towards the center of their bin
using the $t$-slope indicated by the dipole fit. With this method, we expect
to correct the intrabin $t$-dependence of order
$\mathcal{O}(t-t_{\textrm{center}})$, where $t_{\textrm{center}}$ is the value
at the center of the bin to which each point belongs. In effect, the
modification is generally very small as $|t-t_{\textrm{center}}|$ is typically
below 0.02 GeV2. A few reach 0.04 GeV2, but only for larger values of $t$,
where the $t$ dependence is reduced. The only noticeable effect of intrabin
$t$-correction arises in the first bin, where the slope in $t$ is pronounced
and the data very numerous, and we consider this therefore as a negligible
effect in the overall uncertainty budget. For this reason, we will not quote
an uncertainty in $t$ in our results.

We fit the elastic form factor using either the dipole form of Eq. (78) or a
$z$-expansion Lee  _et al._ (2015) with three free-parameters,

$F(t)=\sum_{k=0}^{2}a_{k}z(t)^{k}\,,\ \
z(t)=\frac{\sqrt{t_{\textrm{cut}}-t}-\sqrt{t_{\textrm{cut}}-t_{0}}}{\sqrt{t_{\textrm{cut}}-t}+\sqrt{t_{\textrm{cut}}-t_{0}}}\,,$
(79)

where $t_{\textrm{cut}}=4m_{\pi}^{2}$ and
$t_{0}=t_{\textrm{cut}}(1-\sqrt{1-t_{\max}/t_{\textrm{cut}}})$. We use
$t_{\max}=-1.4$ GeV2. We observe that the specific value of $t_{\textrm{cut}}$
and $t_{\max}$ does not produce a significant effect on the extraction. The
final fit results are produced in Table 2 and displayed in Fig. 9. We quote in
Table 2 both a statistical uncertainty and an excited state uncertainty
obtained by comparing the results with a cut of $t_{s}=3$ and 4 on the
correlation functions. The uncertainty band in Fig. 9 is constructed as the
union of the bands obtained using the two cuts. As one can observe, both the
statistical and excited state uncertainties are extremely narrow for the
elastic form factor, to the point that the model dependence in the
extrapolation to $t=0$ becomes a significant source of uncertainty as well.
The $\chi^{2}$ per degree of freedom of those fits is bad, between 6 and 30.
We attribute that to the fact that the precision of the measurements is such
that lattice spacing errors become very noticeable and make it apparent that
the Lorentz decomposition is violated on the lattice. As a result, the points
do not fall exactly on a single curve, and increasing the statistics would
probably only result in even worse fits.

|  $A$ (dipole at $t=0$) |  $\Lambda$ (MeV) |  $z$-exp order 2 at $t=0$
---|---|---|---
$F_{1}(t)$ local | $0.97\pm 0.004\pm 0.006$ | $1246\pm 3\pm 6$ | $0.989\pm 0.004\pm 0.008$
$F_{1}(t)$ non-local | $0.969\pm 0.005\pm 0.010$ | $1248\pm 4\pm 13$ | $0.988\pm 0.005\pm 0.012$
$F_{2}(t)$ local | $3.53\pm 0.01\pm 0.06$ | $966\pm 2\pm 7$ | $3.38\pm 0.01\pm 0.03$
$F_{2}(t)$ non-local | $3.44\pm 0.01\pm 0.03$ | $982\pm 2\pm 4$ | $3.30\pm 0.01\pm 0.02$
Table 2: Values of the fit of the $F_{1}$ and $F_{2}$ form factors both from
local and non-local matrix elements with a dipole and a $z$-expansion up to
order 2 (three free parameters). The statistical and an estimate of the
systematic uncertainty linked to excited state contamination are provided.
Figure 9: Fit of the elastic form factors from local matrix elements after
binning into 15 values. The uncertainty bands are not only obtained from the
statistical uncertainty of the data but also by merging the fit results using
two cuts in lattice time separation to evaluate excited state contaminations.
The results of the fit are produced in the main body.

###  5.2 Elastic form factors from non-local data

Although elastic form factors are typically computed using the local operator
$z=0$, they can be accessed using non-local operators as well. In fact, some
kinematics $(\vec{p}_{f},\vec{p}_{i})$ are particularly convenient for that
purpose: if $p_{f,3}=p_{i,3}=0$ and $z=(0,0,0,z_{3})$, then both $\nu$ and
$\nu\xi$ are identically zero whatever the value of $z^{2}$, and one has:

$\displaystyle F_{1}(t)=\lim_{z^{2}\rightarrow
0}\mathcal{A}_{1}(\nu=0,\nu\xi=0,t,z^{2})\,,$ (80) $\displaystyle
F_{2}(t)=\lim_{z^{2}\rightarrow 0}\mathcal{A}_{4}(\nu=0,\nu\xi=0,t,z^{2})\,.$
(81)

Using non-local operators exposes us to the risk of introducing additional
uncertainty due to the evaluation of the right-hand sides at $z^{2}<0$. We
plot on the left panel of Fig. 10 the value of
$\mathcal{A}_{1}(\nu=0,\nu\xi=0,t=-1.35$ GeV${}^{2},z^{2})$ using
$\vec{p}_{f}=(1,-1,0)$ and $\vec{p}_{i}=(-1,1,0)$ for every value of $z$ from
0 to $6a$ (we average the $\pm z$ data) using the cut $t_{s}=3$. Evidently, a
noticeable $z$ dependence can be observed.

Figure 10: (left) Value of
$\mathcal{A}_{1}(\nu=0,\nu\xi=0,-1.35\,\textrm{GeV}^{2},z^{2})$ for the pair
of momenta $\vec{p}_{f}=(1,-1,0)$ and $\vec{p}_{i}=(-1,1,0)$ in lattice units,
using the cut $t_{s}=3$. We show both fits by a constant (blue and green) and
a quadratic model $\alpha+\beta(z/a)^{2}$ (orange and red), using either all
data points (solid area) or excluding the local data $z=0$ (hatches). –
(right) To understand why the constant fit ends up below every single point,
we show the full set of measurements in the diagonalized space of the
covariance matrix. See explanations in the text.

It is important to notice that the various points on the plot are highly
correlated with one another, as they all stem from the same external proton
states computed on the same gauge configurations. The correlation matrix of
the data is given by:

$\textrm{Corr}=\begin{pmatrix}1&0.96&0.93&0.86&0.76&0.64&0.51\\\
0.96&1&0.98&0.93&0.84&0.73&0.60\\\ 0.93&0.98&1&0.98&0.91&0.82&0.70\\\
0.86&0.93&0.98&1&0.98&0.91&0.80\\\ 0.76&0.84&0.91&0.98&1&0.97&0.90\\\
0.64&0.73&0.82&0.91&0.97&1&0.97\\\
0.51&0.60&0.70&0.80&0.90&0.97&1\end{pmatrix}.$ (82)

Neighboring points on the plot are correlated by more than 96%, and the
furthest are still over 50%. Those circumstances are common in lattice QCD
plots where the axis variable is the non-local separation $z$. This can
produce counter-intuitive results.

For pedagogical purpose on the counter-intuitive effect of very high
correlation, despite the clear trend in $z$ that appears visually, one could
wish to fit our dataset by a constant. Intuitively, one would expect the
constant to fall somewhere in the middle of the plot, around $F_{1}\approx
0.31$. Such would indeed be the case if our points were uncorrelated. But the
best fit by a constant is in fact $F_{1}\approx 0.28$, strangely close to the
result at $z=0$, and in good agreement with the projection to $z=0$ of the
quadratic fit by $\alpha+\beta(z/a)^{2}$. Excluding the local data at $z=0$
produces exactly the same fit results, either for the constant or quadratic
fit. However, now the constant fit is incompatible with every single point of
the fitted dataset. The remarkable agreement of the constant fit with the
quadratic one at $z=0$ is quite surprising and happens in a general fashion in
other kinematics as well. Although we will use the quadratic fit for the non-
local elastic form factor, it is interesting to try to understand why a value
that is significantly below every single data point can be the best constant
fit of the dataset. This illustrates how such plots can be quite misleading if
the very high degree of correlation is not taken into account.

One needs to remove correlations from the data. First, we diagonalize the
covariance matrix $\textrm{Cov}=P^{-1}DP$ where $D$ is the diagonal matrix of
eigenvalues. Then the vector

$P\times\begin{pmatrix}\mathcal{A}_{1}(z=0)\\\ \cdots\\\
\mathcal{A}_{1}(z=6a)\end{pmatrix}$ (83)

is an uncorrelated random vector. However, at this stage, its central values
are difficult to interpret. In our case, it is interesting to compare these
uncorrelated data with what would happen if $\mathcal{A}_{1}$ was indeed
independent of $z$. Therefore, we form:

$R=P\times\begin{pmatrix}\mathcal{A}_{1}(z=0)\\\ \cdots\\\
\mathcal{A}_{1}(z=6a)\end{pmatrix}\bigg{/}P\times\begin{pmatrix}1\\\ \cdots\\\
1\end{pmatrix}\,,$ (84)

where the division is understood as a term-by-term ratio of the elements in
both vectors. If $\mathcal{A}_{1}$ was indeed constant in $z$, its
uncorrelated version $R$ would remain constant with the same value. The values
of the seven uncorrelated random variables contained in R are shown on the
right panel of Fig. 10.

Now the origin of our problem appears clearly. There exist uncorrelated linear
combinations of our measurements which favor very small constant fits. In
fact, it is easy to see from this example that, by carefully crafting the
correlations between the data, the best fit from a given model can be
arbitrarily far from any of the correlated measurements.

With this in mind, we can extract the elastic form factor using only non-local
matrix elements, and performing quadratic fits over $z\in\\{a,...,6a\\}$. When
using kinematics that are not the special $p_{f,3}=p_{i,3}=0$, since $\nu\neq
0$ and $\xi\neq 0$, the matrix element contains more than just the elastic
form factor. But when the higher-order moments are fitted jointly, as we will
describe in the next section, the entire kinematic dataset can be used to
constrain the elastic form factor using the non-local data. The result is
shown in Fig. 11. In light grey, the drifting value of $\mathcal{A}_{1,4}$
with increasing $z$ is depicted, while its value at $z=0$ is represented by
the red dots. The fit of the grey dots by a quadratic is depicted as the blue
dots, barely distinguishable from the $z=0$ data, or the result of a constant
fit. The final fit by a $z$ expansion of either the local (red) or non-local
(blue) elastic form factor extraction is completely compatible, although the
$\chi^{2}$ per d.o.f. are typically significantly smaller when using the non-
local extraction. We report the results of the non-local extraction of elastic
form factors in Table 2 as well, and the dipole result serves as the standard
result reported at the beginning of this document in Fig. 2. One will notice
that the pattern of $z^{2}$ dependence is clearly different between the
$F_{1}$ and $F_{2}$ elastic form factors.

Figure 11: In grey, the series of 6 points $\mathcal{A}_{1,4}(\nu=0,t,z^{2})$
for $z\in\\{a,...,6a\\}$, in blue a quadratic (correlated) fit of the grey
points, in red the data at $z=0$. The grey points are shifted gradually in $t$
to allow us to observe their trend. For instance, at large $t$, there is no
$z^{2}$ drift in $F_{2}$ on the right panel. The fit by a $z$-expansion (79)
is very similar whether one uses the local or non-local extraction. We merge
the two cuts in lattice time separation.

##  6 Generalized Form Factors

From the discussion in section 2.3, we expect that

$\displaystyle\textrm{Re}\,\mathcal{A}_{1}(\nu,\xi,t,z^{2})$
$\displaystyle=F_{1}(t)-\frac{\nu^{2}}{2}\bigg{(}A_{3,0}(t,z^{2})+\xi^{2}A_{3,2}(t,z^{2})\bigg{)}+\mathcal{O}(\nu^{4})+\mathcal{O}(\Lambda_{\rm
QCD}^{2}z^{2},tz^{2})\,,$ (85)
$\displaystyle\textrm{Im}\,\mathcal{A}_{1}(\nu,\xi,t,z^{2})$
$\displaystyle=-\nu
A_{2,0}(t,z^{2})+\frac{\nu^{3}}{6}\bigg{(}A_{4,0}(t,z^{2})+\xi^{2}A_{4,2}(t,z^{2})\bigg{)}+\mathcal{O}(\nu^{5})+\mathcal{O}(\Lambda_{\rm
QCD}^{2}z^{2},tz^{2})\,,$ (86)

up to power corrections and lattice spacing systematic uncertainties. A
similar relation is obtained by substituting the appropriate combination to
form the $E$ GPD $\mathcal{A}_{4}+\nu\mathcal{A}_{6}-2\xi\nu\mathcal{A}_{7}$
and its $B_{n,k}$ generalized form factors. The question of the $D$-term
extracted from $\mathcal{A}_{5}$ is discussed in section 6.4. We fit these
relation at a fixed value of $z$ in each bin, with an intrabin $t$ correction
proceeding in the same fashion as what we have described for the elastic form
factors. The result using all points in the bin $t\approx-0.33$ GeV2, $z=6a$,
and the cut $t_{s}=3$ is displayed in Fig. 12.

Figure 12: The values of $\mathcal{A}_{1}(\nu,|\xi|,t\approx-0.33$
GeV${}^{2},z=6a)$ with a cut $t_{s}=3$ inside a single bin for the real (left)
and imaginary (right) parts, along with the results of fits by the functional
forms of Eqs. (85) and (86). We take the absolute value of $\xi$ to simplify
the plot as we expect $\mathcal{A}_{1}$ to be even in $\xi$.

The $\chi^{2}$ per d.o.f. of this fit of 28 points with 3 parameters is close
to 2 both for the real and imaginary parts. Note that although it may seem
that the special kinematic point at $\nu=0$ corresponding to
$\vec{p}_{f}=(0,0,0)$ and $\vec{p}_{i}=(1,1,0)$ is dragging the value of
$F_{1}(t)$ upwards, this is not the case. Excluding this point from the fit
would give $F_{1}(t,z=6a)=0.687\pm 0.008$, whereas with it included, we find
the very similar $F_{1}(t,z=6a)=0.689\pm 0.007$. As we have already discussed,
the data favor a value of $F_{1}(t)$ at $z=6a$ that is noticeably larger than
the data at $z=0$ prefers, namely $0.652\pm 0.007$, that is $5\sigma$ below.
One will also notice that data at larger $\xi$ exists at smaller values of
$\nu$ since part of the average $P_{3}$ momentum is sacrificed to create the
momentum transfer. The coverage in $(\xi,\nu)$ available for our study is
depicted in Fig. 13, with the exception of a few points where
$p_{f,3}=-p_{i,3}\neq 0$, corresponding to $\nu=0$ and $\xi$ infinite (but
$\nu\xi\neq 0$). As one can see from Eqs. (85) and (86), those kinematics
probe exclusively the elastic form factor and the skewness-dependent
generalized form factors, so they are particularly useful.

Figure 13: The distribution of available values of $(\xi,\nu)$ where $\nu$ is
presented at $z=a$.

###  6.1 Matching

At higher order and in contrast to the elastic form factors, the generalized
form factors are scale-dependent, and we have derived in section 2.3 a
perturbative prediction for their $z^{2}$ dependence. The matching of the
gravitational form factors $A_{2,0},B_{2,0}$ is particularly interesting to
observe because our data are very precise. Let us use again the same bin as in
the previous discussion. We try three different leading-order/leading-
logarithmic matching procedures:

  1. 1.

The procedure that arises immediately from the fixed order matching in Eq.
(17):

$\displaystyle A_{2,0}(t,\mu^{2})$
$\displaystyle=A_{2,0}(t,z^{2})\bigg{[}1+\frac{\alpha_{s}(\mu^{2})C_{F}}{2\pi}\left(\gamma_{1}L+d_{1}\right)\bigg{]}\,,$
(87)

where $\gamma_{1}=-4/3$ is the DGLAP anomalous dimension, $d_{1}=14/3$ the
moment of the non-logarithmic part of the matching kernel, and
$L=\ln(e^{2\gamma_{E}+1}/4\times(z^{2}/\lambda)\times\mu^{2})$. Here $\lambda$
should be close to $-1$ and will be varied in the interval $[-1/2,-2]$ as is
usual to evaluate uncertainty in scale fixing.

  2. 2.

A variant in which we evaluate $\alpha_{s}$ using its leading-logarithmic
3-flavor evolution at the “natural” scale of the data $\lambda/z^{2}$:

$\displaystyle A_{2,0}(t,\mu^{2})$
$\displaystyle=A_{2,0}(t,z^{2})\bigg{[}1+\frac{\alpha_{s}(\lambda/z^{2})C_{F}}{2\pi}\left(\gamma_{1}L+d_{1}\right)\bigg{]}\,.$
(88)

  3. 3.

A full leading-logarithmic matching as derived in section 2.3:

$\displaystyle A_{2,0}(t,\mu^{2})$
$\displaystyle=A_{2,0}(t,z^{2})\left(\frac{\alpha_{s}(\lambda/z^{2})}{\alpha_{s}(\mu^{2})}\right)^{\gamma_{1}C_{F}/(2\pi\beta_{0})}\bigg{[}1+\frac{\alpha_{s}(\lambda/z^{2})C_{F}}{2\pi}\left(\gamma_{1}\ln\left(\frac{e^{2\gamma_{E}+1}}{4}\right)+d_{1}\right)\bigg{]}\,.$
(89)

We use the $\overline{MS}$ scale $\mu=2$ GeV and present the result of the
three matching procedures in Fig. 14. We use $\alpha_{s}(\mu^{2})=0.28$, a
value compatible with the usual $\overline{MS}$ value at that scale. The
relative arbitrariness of this choice is balanced by the possibility of
performing scale variations.

Figure 14: The gravitational form factor of the $H$ GPD at $t=-0.33$ GeV2 (cut
$t_{s}=3$) is an increasing function of the separation for
$z\in\\{a,...,6a\\}$ (blue points). We show three different attempts at
matching it to an $\overline{MS}$ scale of 2 GeV, using $\alpha_{s}(2\textrm{
GeV})=0.28$ according to the three matching formulas (87) - (89). The greyish
points represent an account of scale variation, which is considerable compared
to the statistical uncertainty in this well-constrained bin.

By construction, the first matching procedure does not present any divergence
except for $-z^{2}\rightarrow\infty$ and can be applied to as large $-z^{2}$
as desired. We observe some inconsistency at small $z$, followed by a
stabilization in the large $z$ regime. The second matching procedure, on the
other hand, cannot be applied once as $\lambda/z^{2}$ reaches the
$\Lambda_{QCD}$ divergence. With $\alpha_{s}(2$ GeV$)=0.28$ and three active
flavors, $\Lambda_{QCD}=165$ MeV, which allows us nonetheless to apply the
matching to all our values of $z^{2}$. However, the large values result in an
extremely large-scale fixing uncertainty, as one would expect.

The full leading-logarithm calculation exhibits the smallest scale fixing
uncertainties of the three matching procedures. Although
$\alpha_{s}(\lambda/z^{2})$ increases considerably at large $z$, since
$\gamma_{1}$ is negative, the ratio
$(\alpha_{s}(\lambda/z^{2})/\alpha_{s}(\mu^{2}))^{\gamma_{1}C_{F}/(2\pi\beta_{0})}$
decreases. On the other hand, the final factor of Eq. (89) increases at large
$z$, and both variations happen in quite similar proportions. The procedure
produces relatively constant results up to $z=0.3$ fm, after which it is no
longer constant in this bin. In principle, we do not have expectations that a
perturbative matching explains the $z^{2}$ dependence of the data when
$-1/z^{2}\ll 1$ GeV2. Therefore, this situation is not particularly
concerning. In fact, when applied to the $E$ GPD in the same bin, the full
leading-logarithmic matching still has the smallest scale fixing
uncertainties, and overall the best behavior as seen in Fig. 15. At this
stage, considering that different behaviors are observed for various bins, no
strong conclusion can be drawn, especially in view of the $z^{2}$ dependence
observed in the elastic form factors.

Figure 15: Same as Fig. 14, except for the $E$ GPD. The statistical
uncertainty is worse than the $H$ GPD.

Only a few bins of the generalized form factor $A_{2,0}$ are precise enough in
terms of statistics and systematic difference between the cut $t_{s}=3$ and 4
that scale variation in the leading-logarithmic matching produces a visible
difference. Overall, scale variation does not produce a significant
contribution to the uncertainty budget given the other uncertainties once the
$t$-dependence is reconstructed in the final fit. We will therefore stick to
the full leading-logarithmic matching with $\lambda=-1$ in the following. Even
if this is not a practical matter at the current stage of our calculations, we
foresee that applying perturbation theory at large $z$ is not good practice
when all uncertainties are under control. On the other hand, in the event that
our framework retains a strong sensitivity to the leading-twist operator even
at scales where perturbation theory becomes dubious, we have advocated that
evolution operators in $z^{2}$ should be computed on the lattice in Dutrieux
_et al._ (2024). However, this program requires a higher confidence in the
suppression of power corrections and other lattice systematic uncertainties
that we do not believe to have been achieved in the present calculation.

We present in Fig. 16 an overall picture of the matched moments in one bin in
$t$. The importance of the large $z$ data (and therefore large Ioffe time) to
constrain the higher moments is obvious.

Figure 16: The moments of the GPDs $H^{u-d}$ (left) and $E^{u-d}$ (right) at
$t=-0.33$ GeV2 with leading-logarithmic matching to $\mu=2$ GeV as a function
of $z/a$ (cut $t_{s}=4$). In the upper row, we show the skewness-independent
moments, and in the lower row the skewness-dependent moments.

###  6.2 $t$-dependent generalized form factors

The results of the fits of the $t$-dependence of moments of GPDs are presented
in Figs. 17 to 19 for the generalized form factors $A_{n,0}$ and $B_{n,0}$
that are independent of $\xi$, and in Fig. 20 for the $A_{n,2}$ and $B_{n,2}$
that are dependent on $\xi$. Since beyond the elastic form factor, the trend
in $z$ of the data is not discernible in a consistent fashion, we use a
constant fit on $z\in\\{a,...,6a\\}$ after having matched the data to
$\overline{MS}$ at 2 GeV using the leading-logarithmic matching derived in
section 2.3. We demonstrate in Appendix C that the correlated constant fit
gives very robust results on a diverse subset of separations, for instance
using only $z\in\\{a,...,4a\\}$. The colored bands are 68% confidence
intervals obtained by dropping the lower and upper 16% of the best-fit results
on the data with a cut in lattice time separation of three. This makes
practically no difference compared to a simple $1\sigma$ uncertainty band,
except for the most poorly determined moments $B_{4,0}$ and $A_{4,2}$,
$B_{4,2}$ where the non-linearity of the dipole fit becomes apparent and the
distributions of best fits are non-Gaussian. The uncolored band results from
the union of the result with a cut of three and a cut of four. It appears that
excited state uncertainty increases as we head towards larger moments,
unsurprisingly as the signal in the data weakens. The obtained results for the
dipole fit in terms of value at $t=0$ and mass are reported in Fig. 2.

Figure 17: The generalized form factors $A_{2,0}$ (left) and $B_{2,0}$ (right)
as a function of $t$. The data points are obtained for a time separation cut
of three in the analysis of correlation functions. The fit of those data
points, either by a dipole or a $z$ expansion, is represented by the uniformly
colored areas. When merging this uncertainty with the one using a time
separation cut of four, one obtains the wider uncertainties, which contain
information about excited state contamination. Figure 18: The generalized form
factors $A_{3,0}$ (left) and $B_{3,0}$ (right) (see caption of Fig. 17)
Figure 19: The generalized form factors $A_{4,0}$ (left) and $B_{4,0}$ (right)
(see caption of Fig. 17)

Figure 20: The skewness-dependent generalized form factors $A_{3,2}$ (top
left), $B_{3,2}$ (top right), $A_{4,2}$ (bottom left) and $B_{4,2}$ (bottom
right) (see caption of Fig. 17).

###  6.3 Radial distributions

One of the main physical motivations behind the study of GPDs is their
characterization of radial distributions of partonic properties inside hadrons
(momentum, energy, pressure, etc.). Such distributions are obtained by a
Fourier transform of the $t$-dependence of the GPD, usually at zero skewness
Burkardt (2003). The impact parameter distribution of an unpolarized quark in
an unpolarized proton is given by

$I(x,\vec{b}_{\perp})=\int\frac{\mathrm{d}^{2}\vec{\Delta}_{\perp}}{(2\pi)^{2}}\,e^{-i\vec{b}_{\perp}\cdot\vec{\Delta}_{\perp}}H(x,\xi=0,t=-\vec{\Delta}^{2}_{\perp})\,,$
(90)

whereas the GPD $E$ allows one to characterize the unpolarized quark
distribution inside a transversely polarized proton,

$I^{T}(x,\vec{b}_{\perp})=\int\frac{\mathrm{d}^{2}\vec{\Delta}_{\perp}}{(2\pi)^{2}}\,e^{-i\vec{b}_{\perp}\cdot\vec{\Delta}_{\perp}}\bigg{[}H(x,\xi=0,t=-\vec{\Delta}^{2}_{\perp})+i\frac{\Delta_{y}}{2m}E(x,\xi=0,t=-\vec{\Delta}^{2}_{\perp})\bigg{]}\,.$
(91)

Using the model-dependent extrapolation to any $t$ provided by the previous
fits, we can extract the first moments of the isovector component of these
distributions,

$\displaystyle I^{u-d}_{n+1}(\vec{b}_{\perp})$
$\displaystyle=\int\frac{\mathrm{d}^{2}\vec{\Delta}_{\perp}}{(2\pi)^{2}}\,e^{-i\vec{b}_{\perp}\cdot\vec{\Delta}_{\perp}}A^{u-d}_{n+1,0}(-\vec{\Delta}^{2}_{\perp})\,,$
(92) $\displaystyle I^{u-d,T}_{n+1}(\vec{b}_{\perp})$
$\displaystyle=\int\frac{\mathrm{d}^{2}\vec{\Delta}_{\perp}}{(2\pi)^{2}}\,e^{-i\vec{b}_{\perp}\cdot\vec{\Delta}_{\perp}}\bigg{[}A^{u-d}_{n+1,0}(-\vec{\Delta}^{2}_{\perp})+i\frac{\Delta_{y}}{2m}B^{u-d}_{n+1,0}(-\vec{\Delta}^{2}_{\perp})\bigg{]}\,.$
(93)

In the following we use the dipole fits as the $t$-dependence of the $z$
expansion is too unconstrained at large $t$ to give rise to sensible radial
distributions. The extraction presented here is therefore strongly model-
dependent. Even then, the $t$-dependence of $A_{4,0}$ and $B_{4,0}$ is very
poorly determined and therefore we only present results for $n\in\\{0,1,2\\}$.

Figure 21: The unpolarized IPD moments as a function of $|\vec{b}_{\perp}|$,
rescaled for visibility. The errorband are obtained from the dipole fits of
the data merging both the cuts of 3 and 4 in Euclidean time.

The unpolarized moments (92) are spherically symmetric in $|\vec{b}_{\perp}|$.
We represent them in Fig. 21. As can already be inferred by the increasing
value of the fitted dipole mass (see the summary table of Fig. 2), when the
order of the skewness-independent moments increases, the radial extent of the
distribution shrinks. This can be largely understood from a kinematic point of
view, as higher-order Mellin moments are increasingly dominated by the large
$x$ domain of the distribution. The center of the proton with respect to which
our radial coordinate system is defined is the barycenter of the partonic
longitudinal momentum. When the longitudinal momentum fraction of the active
parton $x$ approaches one, the active parton is by necessity close to the
center of the proton, so the impact parameter distribution at large $x$ is
very narrow.

The transverse impact parameter distribution moments (93) are not spherically
symmetric. We represent them, along with the unpolarized moments, in Fig. 22.
We only use the central value of the dipole fits. A subtle feature of the
transverse plots, which can also be observed in Bhattacharya  _et al._
(2023b), is the fact that moments with $n$ even (extracted from the imaginary
part of the Ioffe-time distributions) seem to show a lesser distortion in the
transverse impact parameter distribution than the odd moments. In Bhattacharya
_et al._ (2023b), where the $u+d$ component is extracted as well, it also
seems to hold for $n=3$. Here again, the distribution shrinks at large $n$.
One can notice that perturbation theory predicts the GPD to be independent of
$t$ when $x\rightarrow 1$ Yuan (2004):

$H^{q}(x,\xi,t)\sim\frac{(1-x)^{3}}{(1-\xi^{2})^{2}}\,,\ \ \
E^{q}(x,\xi,t)\sim\frac{(1-x)^{5}}{(1-\xi^{2})^{2}}\,.$ (95)

This is equivalent to the large $n$ GPD moments becoming independent of $t$,
or the impact parameter distribution becoming a simple Dirac peak at
$\vec{b}_{\perp}=0$.

Figure 22: The moments $n=0,1,2$ (respectively first, second and third row) of
the unpolarized $u-d$ quark combination in an unpolarized proton (left) and in
a transversely polarized proton (right).

###  6.4 $D$-term

For the isovector $D$-term, Eq. (45) invites us to fit:

$\textrm{Im}\,\mathcal{A}_{5}(\nu,\xi,t,z^{2})=\xi\nu
C_{2}(t,z^{2})+\mathcal{O}(\nu^{3})+\mathcal{O}(\Lambda_{\rm
QCD}^{2}z^{2},tz^{2})\,.$ (96)

As we have noticed before, there is a trade-off between the value of $\xi$ and
$\nu$ as larger $\xi$ requires sacrificing some of the $P_{3}$ momentum to
introduce a momentum transfer. As a result, the range in $\nu\xi$ available to
extract the $D$-term does not exceed 1.2 when $z=6a$, and we can only achieve
sensitivity to the $C_{2}$ moment of the $D$-term. The situation in one of the
bins in $t$ is depicted in Fig. 23. Because of the absence of any discernable
trend in $t$, we do not attempt an intrabin $t$ correction in the fashion we
used for the rest of the GPD.

Figure 23: The values of $\mathcal{A}_{5}(\nu,\xi,t\approx-1.36$
GeV${}^{2},z=6a)$ with a cut $t_{s}=3$ inside a single bin for the real (left)
and imaginary (right) parts, along with the results of fits by the functional
form of Eq. (96). We present the data as a function of $\xi\nu$ instead of the
separate $\xi$ and $\nu$.

In general, $\mathcal{A}_{5}$ is compatible with 0 and small in magnitude.
Some kinematics $(\vec{p}_{f},\vec{p}_{i})$ show noticeably non-zero values in
the real part, a problem we have already identified when discussing the
contribution of $\mathcal{A}_{5}$ in the local matrix element. Therefore,
there is a clear sign of potential lattice discretization errors, which we
have no means of subtracting from the data within this analysis. We also note
that Ref. Radyushkin (2024) highlights the potential presence of “kinematic”
higher-twists in $\mathcal{A}_{5}$. When using the full statistics available
in some bins, it is possible to fit the imaginary part with the functional
form of Eq. (96) and obtain an extraction for $C_{2}$ only a few standard
deviations away from 0. The imaginary part of $\mathcal{A}_{5}(\nu,\xi)$ is
compatible with being an odd function of the product $\nu\xi$, which mirrors
the fact that the $D$-term is an odd function of $x/\xi$. The matching of
$C_{2}$ then proceeds exactly similarly to the one of $A_{2,0}$ and $B_{2,0}$
(45).

The final result on $C_{2}$ is depicted in Fig. 24 matched to 2 GeV using the
data $z=a$ up to $6a$. The dipole and the $z$-expansion fit produce very
different $t$-dependence and extrapolation in the limit where $t\rightarrow
0$, due to the weakness of the constraint imposed by the data. The small
magnitude of $C_{2}$ is consistent with the results of other lattice
calculations from local matrix elements, e.g. Hagler  _et al._ (2008); Bali
_et al._ (2019); Alexandrou  _et al._ (2020a). In the large $N_{c}$
approximation, one can note that the $u-d$ contribution to the $D$-term is
suppressed by a factor $1/N_{c}$ compared to the $u+d$ flavor combination
Goeke  _et al._ (2001).

Figure 24: The isovector $D$-term form factor $\int d\alpha\,\alpha
D(\alpha)$.

##  7 Conclusions

In this work, we present the calculation of off-forward nucleon matrix
elements with extensive coverage in initial and final momenta. These matrix
elements scan a wide range of $t$ and $\xi$ relevant for understanding the
three-dimensional tomography of the nucleon, through a factorization
relationship to GPDs analogous to those for DVCS and DVMP. This lattice
calculation was made possible by utilizing the technique of distillation to
explicitly project onto the momentum states for all operators in the
correlation functions. In addition to significantly improving the signal, this
procedure provides a natural method for an efficient calculation by allowing
many momenta to be calculated, while recycling the components for each
operator.

With the current statistical precision, control of excited state contamination
is critical for an accurate analysis. We have compared results from two
different methods for extracting matrix elements and from different cuts on
the data. This step is critical for creating realistic error estimates for the
matrix elements. The data with different cuts in Euclidean time can lead to
significant deviations in the final results, as was seen in the extraction of
generalized form factors. In future studies, we will use the distillation to
implement distinct nucleon operators for each momentum. This allows for a GEVP
analysis to explicitly remove contributions from excited states. The excited
state contamination will only become more difficult to control in the chiral
limit, so this full control of excited state contamination will be crucial for
physical extrapolations.

By using this Lorentz-invariant approach, the $\nu$ dependence of this data
can be analyzed to study the Mellin moments of the GPDs. In this study, we
have utilized relatively small momenta in the direction of the Wilson line,
which scans the low $\nu$ region that is dominated by the lowest few moments.
We extract the moments of $H$ and $E$ and study the $t$ dependence and, for
the first time, the $\xi$ dependence. The $\xi^{2}$ coefficients of the
moments $A_{n,2}$ and $B_{n,2}$ are larger than the $\xi^{0}$ coefficients
$A_{n,0}$ and $B_{n,0}$. We also extracted radial distributions at zero
skewness by performing the Fourier transform of the $t$ dependence.

The use of non-local separations in the matrix elements up to $z=0.6$ fm gives
rise to challenges for the proper evaluation of the systematic uncertainty
caused by the perturbative matching, as well as the contamination by higher
twist contributions, which we do not have the possibility to properly treat in
this study. In view of some clear signal of lattice discretization errors,
which we cannot easily differentiate from higher twist contributions, we see
the realization of a continuum limit of this study as an important further
milestone.

This calculation represents a step on the road to a systematically controlled
calculation of nucleon tomography through GPDs. Nucleon tomography requires
significantly more matrix elements, to scan the full three-dimensional space
of GPDs than determinations of PDFs, particularly for a full study of
systematic errors. We have demonstrated a method capable of tackling the
problem of kinematical coverage and for studying the significant impact of
excited state contamination.

##  8 Acknowledgments

This project was supported by the U.S. Department of Energy, Office of
Science, Contract No. DE-AC05-06OR23177, under which Jefferson Science
Associates, LLC operates Jefferson Lab. KO and HD were supported in part by
U.S. DOE Grant #DE-FG02-04ER41302. CM is supported in part by U.S. DOE ECA
#DE-SC0023047. AR acknowledges support by U.S. DOE Grant #DE-FG02-97ER41028.
The work of HD and DR was conducted in part under the Laboratory Directed
Research and Development Program (LDRD 2412) at Thomas Jefferson National
Accelerator Facility for the U.S. Department of Energy. This work has
benefited from the collaboration enabled by the Quark-Gluon Tomography (QGT)
Topical Collaboration, U.S. DOE Award #DE-SC0023646. This research was funded,
in part, by l’Agence Nationale de la Recherche (ANR), project
ANR-23-CE31-0019. For the purpose of open access, the authors have applied a
CC-BY public copyright licence to any Author Accepted Manuscript (AAM) version
arising from this submission. Computations for this work were carried out in
part on facilities of the USQCD Collaboration, which are funded by the Office
of Science of the U.S. Department of Energy. This work was performed in part
using computing facilities at William & Mary which were provided by
contributions from the National Science Foundation (MRI grant PHY-1626177),
and the Commonwealth of Virginia Equipment Trust Fund. This work used the
Extreme Science and Engineering Discovery Environment (XSEDE), which is
supported by National Science Foundation grant number ACI-1548562.
Specifically, it used the Bridges system, which is supported by NSF award
number ACI-1445606, at the Pittsburgh Supercomputing Center (PSC) Towns  _et
al._ (2014); Nystrom  _et al._ (2015). In addition, this work used resources
at NERSC, a DOE Office of Science User Facility supported by the Office of
Science of the U.S. Department of Energy under Contract #DE-AC02-05CH11231, as
well as resources of the Oak Ridge Leadership Computing Facility at the Oak
Ridge National Laboratory, which is supported by the Office of Science of the
U.S. Department of Energy under Contract No. #DE-AC05-00OR22725. The software
codes Chroma Edwards and Joo (2005), QUDA Clark  _et al._ (2010); Babich  _et
al._ (2010), QPhiX Joó  _et al._ (2016), and Redstar Chen  _et al._ (2023)
were used in our work. The authors acknowledge support from the U.S.
Department of Energy, Office of Science, Office of Advanced Scientific
Computing Research and Office of Nuclear Physics, Scientific Discovery through
Advanced Computing (SciDAC) program, and of the U.S. Department of Energy
Exascale Computing Project (ECP). The authors also acknowledge the Texas
Advanced Computing Center (TACC) at The University of Texas at Austin for
providing HPC resources, like Frontera computing system Stanzione  _et al._
(2020) that has contributed to the research results reported within this
paper. The authors acknowledge William & Mary Research Computing for providing
computational resources and/or technical support that have contributed to the
results reported within this paper.

We acknowledge the EuroHPC Joint Undertaking for awarding this project access
to the EuroHPC supercomputer LUMI, hosted by CSC (Finland) and the LUMI
consortium through a EuroHPC Extreme Scale access call. This work also
benefited from access to the Jean Zay supercomputer at the Institute for
Development and Resources in Intensive Scientific Computing (IDRIS) in Orsay,
France under project A0080511504.

##  Appendix A Full set of Lorentz amplitudes

We present here the full set of Lorentz amplitudes $\mathcal{A}_{k}$ for
$k\in\\{1,...,8\\}$ following the decomposition of Eq. (10) for two pairs of
momentum: $\vec{p}_{f}=(2,1,0)$ and $\vec{p}_{i}=(0,1,2)$ (lattice units) on
Fig. 25, and $\vec{p}_{f}=(1,-1,0)$ and $\vec{p}_{i}=(-1,1,0)$ on Fig. 27. The
data is presented as a function of $z$ in lattice units. We only show positive
$z$ as we average over $\pm z$ using the symmetries discussed in section 2.4.
We remind the reader that this means that each data point exists therefore at
a different “scale” of the order of $-1/z^{2}$. The amplitudes
$\mathcal{A}_{2,3,8}$ which do not enter the construction of light-cone GPDs
are represented together in the last plots, whereas the other five get their
individual plot. Let us comment on some interesting features of both figures.