CSSM/QCDSF/UKQCD Collaborations

# Reconstructing generalised parton distributions from the lattice off-forward
Compton amplitude

A. Hannaford-Gunn  CSSM, Department of Physics, The University of Adelaide,
Adelaide SA 5005, Australia  K. U. Can  CSSM, Department of Physics, The
University of Adelaide, Adelaide SA 5005, Australia  J. A. Crawford  CSSM,
Department of Physics, The University of Adelaide, Adelaide SA 5005, Australia
R. Horsley  School of Physics and Astronomy, University of Edinburgh,
Edinburgh EH9 3JZ, UK  P. E. L. Rakow  Theoretical Physics Division,
Department of Mathematical Sciences, University of Liverpool, Liverpool L69
3BX, UK  G. Schierholz  Deutsches Elektronen-Synchrotron DESY, Notkestr. 85,
22607 Hamburg, Germany  H. Stüben  Regionales Rechenzentrum, Universität
Hamburg, 20146 Hamburg, Germany  R. D. Young  CSSM, Department of Physics, The
University of Adelaide, Adelaide SA 5005, Australia  J. M. Zanotti  CSSM,
Department of Physics, The University of Adelaide, Adelaide SA 5005, Australia

###### Abstract

We present a determination of the structure functions of the off-forward
Compton amplitude $\mathcal{H}_{1}$ and $\mathcal{E}_{1}$ from the Feynman-
Hellmann method in lattice QCD. At leading twist, these structure functions
give access to the generalised parton distributions (GPDs) $H$ and $E$,
respectively. This calculation is performed for an unphysical pion mass of
$m_{\pi}=412\;\text{MeV}$ and four values of the soft momentum transfer,
$t\approx 0,-0.3,-0.6,-1.1\;\text{GeV}^{2}$, all at a hard momentum scale of
$\bar{Q}^{2}\approx 5\;\text{GeV}^{2}$. Using these results, we test various
methods to determine properties of the real-time scattering amplitudes and
GPDs: (1) we fit their Mellin moments, and (2) we use a simple GPD ansatz to
reconstruct the entire distribution. Our final results show promising
agreement with phenomenology and other lattice results, and highlight specific
systematics in need of control.

††preprint: ADP-24-08/T1247††preprint: DESY-24-065††preprint: Liverpool LTH
1370

##  I Introduction

Generalised parton distributions (GPDs) are among the most important hadron
structure observables, giving access to the spatial probability distributions
of quarks Burkardt (2000), as well as the spin Ji (1997), force and pressure
distributions within a hadron Polyakov (2003); Burkert  _et al._ (2018).

While GPDs can in principle be determined from hard exclusive scattering
processes, in practice such experimental determinations face certain
difficulties Kumerički  _et al._ (2016); Guidal  _et al._ (2013); Bertone  _et
al._ (2021). In this context, first principles calculations of GPDs from
lattice QCD can be extremely useful in both guiding and comparing to
experiment. There have already been a handful of studies exploring the ability
of lattice results to aid experimental GPD fits Dutrieux  _et al._ (2022); Guo
_et al._ (2022, 2023); Riberdy  _et al._ (2024); Mamo and Zahed (2024).

While lattice calculations cannot directly determine GPDs, they can be used to
reconstruct them. Historically, this has been limited to determinations of the
first three Mellin moments of GPDs, calculated from local twist-two operators
Hägler  _et al._ (2003); Göckeler  _et al._ (2004, 2005, 2007); Hägler  _et
al._ (2008); Brömmel  _et al._ (2008); Bratt  _et al._ (2010); Alexandrou  _et
al._ (2011). Recent advances since 2015 have opened new avenues for
reconstructing GPDs from non-local operators, notably the quasi- Ji (2013) and
pseudo-distribution Radyushkin (2017) approaches, which have yielded promising
initial determinations of GPDs Chen  _et al._ (2020); Lin (2021); Alexandrou
_et al._ (2020); Bhattacharya  _et al._ (2023); Lin (2023).

In this paper, we apply the Feynman-Hellmann method to determine the off-
forward Compton amplitude (OFCA) and thereby access GPDs. A major distinction
between this method and those mentioned above is that we calculate a lattice
version of the scattering amplitude from which GPDs are determined
experimentally. As such, we have the potential to calculate phenomenologically
interesting properties that are inaccessible to other methods.

So far, the _forward_ Compton amplitude has been the focus of the Feynman-
Hellmann calculations Chambers  _et al._ (2017); Can  _et al._ (2020),
including studies of the scaling behaviour and higher-twist structures
Hannaford-Gunn  _et al._ (2020); Batelaan  _et al._ (2023) and the subtraction
function Hannaford-Gunn  _et al._ (2022a). Similarly, a calculation of the
off-forward Compton amplitude could determine the $\bar{Q}^{2}$ scaling of
this amplitude, of which there have only been limited experimental studies
Camacho Muñoz  _et al._ (2006); Defurne  _et al._ (2015, 2017); Georges
(2018). For the OFCA the $\bar{Q}^{2}$ scaling behaviour could be extremely
useful in isolating the leading-twist contribution, as most deeply virtual
Compton scatting experiments have a relatively modest hard scale of
$Q^{2}\approx 1-12\;\text{GeV}^{2}$ and contain additional $|t|/Q^{2}$
corrections Braun  _et al._ (2012); Guo  _et al._ (2021). Moreover, this
method can determine the off-forward subtraction function, which is a key
input for experimental determinations of the proton pressure distribution
Burkert  _et al._ (2018); Kumerički (2019).

In this work we focus on determining the subtracted structure functions of the
OFCA for a single hard scale, $\bar{Q}^{2}\approx 5\;\text{GeV}^{2}$. We build
upon our previous paper on the OFCA Hannaford-Gunn  _et al._ (2022b), where we
developed the formalism for this calculation and presented exploratory
numerical results. Here, we present major improvements on this numerical
calculation: we separately determine the $\mathcal{H}_{1}$ and
$\mathcal{E}_{1}$ structure functions, which can be uniquely related to twist-
two GPDs, over a wider range of kinematics. These improvements allow us to
focus on modelling and determining GPD properties: (1) we determine the Mellin
moments of $\mathcal{H}_{1}$ and $\mathcal{E}_{1}$ with a largely model-
independent fit, and (2) fit our Compton amplitude results using a
phenomenologically motivated GPD ansatz.

##  II Background

Our purpose in this paper is to determine the off-forward Compton amplitude
(OFCA) and its structure functions from lattice QCD. The notation of this
paper follows that of Ref. Hannaford-Gunn  _et al._ (2022b), in which can be
found derivations and full expressions for analytic results used here.

The OFCA is defined as

$\displaystyle T_{\mu\nu}\equiv i\int d^{4}ze^{\frac{i}{2}(q+q^{\prime})\cdot
z}\langle P^{\prime}|T\\{j_{\mu}(z/2)j_{\nu}(-z/2)\\}|P\rangle,$ (1)

where $j_{\mu}$ is the electromagnetic current operator. See Fig. 1 for the
diagram of this process.

Figure 1: Diagram of off-forward
$\gamma^{*}(q)N(P)\to\gamma^{*}(q^{\prime})N(P^{\prime})$ scattering.

In terms of kinematics, we use the basis of momentum vectors

$\displaystyle\bar{P}=\frac{1}{2}(P+P^{\prime}),\quad\bar{q}=\frac{1}{2}(q+q^{\prime}),\quad\Delta=P^{\prime}-P,$
(2)

which gives us two scaling variables and two momentum scales:

$\displaystyle\bar{\omega}=\frac{2\bar{P}\cdot\bar{q}}{-\bar{q}^{2}},\quad\xi=\frac{\Delta\cdot\bar{q}}{2\bar{P}\cdot\bar{q}},\quad\bar{Q}^{2}=-\bar{q}^{2},\quad
t=\Delta^{2}.$ (3)

For this work, we focus on zero-skewness kinematics, $\xi=0$, which is
enforced by choosing $\bar{q}\cdot\Delta=0$.

The general OFCA is parameterised by 18 linearly independent off-forward
structure functions Perrottet (1973); Tarrach (1975)

$\displaystyle{T}_{\mu\nu}$
$\displaystyle=-\frac{1}{2\bar{P}\cdot\bar{q}}\Big{(}h\cdot\bar{q}\mathcal{H}_{1}+e\cdot\bar{q}\mathcal{E}_{1}\Big{)}\left(g_{\mu\nu}-\frac{q^{\prime}_{\mu}q_{\nu}}{q\cdot
q^{\prime}}\right)+\dots$ (4)

where we have defined the Dirac bilinears

$h_{\mu}=\bar{u}(P^{\prime})\gamma_{\mu}u(P)\text{, and
}e_{\mu}=\bar{u}(P^{\prime})\frac{i\sigma_{\mu\kappa}\Delta^{\kappa}}{2m_{N}}u(P).$
(5)

See Ref. Hannaford-Gunn  _et al._ (2022b) for a full parameterisation of the
Compton amplitude; the tensor decomposition with all zero-skewness structures
is given in Eq. (43). For this work, we are interested in the helicity-
conserving $\mathcal{H}_{1}$ and helicity-flipping $\mathcal{E}_{1}$ structure
functions, which respectively give access to the $H$ and $E$ GPDs at leading-
twist. Our structure functions $\mathcal{H}_{1}$ and $\mathcal{E}_{1}$ are
comparable to the $\mathcal{H}$ and $\mathcal{E}$ Compton form factors used in
experimental analyses Belitsky  _et al._ (2002, 2014); Braun  _et al._ (2014).

The Euclidean Compton amplitude we calculate in lattice QCD can only be
related to the Minkowski amplitude if we use the kinematics $|\bar{\omega}|<1$
Can  _et al._ (2020). However, the real-time scattering amplitude at zero-
skewness has the kinematics $|\bar{\omega}|>1$.

As such, we use a dispersion relation to connect our lattice results
${\mathcal{H}}_{1}(\bar{\omega},t,\bar{Q}^{2})$ to the real-time structure
function, $\mathcal{H}_{1}(x,t,\bar{Q}^{2})$:

$\displaystyle\begin{split}\overline{\mathcal{H}}_{1}(\bar{\omega},t,\bar{Q}^{2})&=\frac{2\bar{\omega}^{2}}{\pi}\int_{0}^{1}dx\frac{x\text{Im}\mathcal{H}_{1}(x,t,\bar{Q}^{2})}{1-x^{2}\bar{\omega}^{2}},\end{split}$
(6)

where
$\overline{\mathcal{H}}_{1}(\bar{\omega},t,\bar{Q}^{2})={\mathcal{H}}_{1}(\bar{\omega},t,\bar{Q}^{2})-{\mathcal{H}}_{1}(0,t,\bar{Q}^{2})$,
the subtracted structure function. An analogous result can be derived for the
$\mathcal{E}_{1}$ structure function.

At $t=0$, the forward limit, $\mathcal{H}_{1}$ reduces to the forward Compton
structure function, $\mathcal{F}_{1}$, and one may use the optical theorem to
write Eq. (6) as

$\displaystyle\overline{\mathcal{F}}_{1}(\omega,{Q}^{2})={4\omega^{2}}\int_{0}^{1}dx\frac{xF_{1}(x,Q^{2})}{1-x^{2}\omega^{2}},$
(7)

where $F_{1}$ is the deep inelastic scattering structure function.

In the leading-twist approximation ($\bar{Q}^{2}\gg\Lambda_{\text{QCD}}^{2}$),
Eq. (6) becomes Hannaford-Gunn  _et al._ (2022b)

$\displaystyle\begin{split}\overline{\mathcal{H}}_{1}(\bar{\omega},t)&=2\bar{\omega}^{2}\int_{-1}^{1}dx\frac{x{H}(x,t)}{1-x^{2}\bar{\omega}^{2}},\end{split}$
(8)

where $H(x,t)$ is the twist-two helicity-conserving GPD. An analogous result
exists for the replacements $\mathcal{H}_{1}\to\mathcal{E}_{1}$ and
${H}\to{E}$.

From Eqs. (6) and (8) we see that it is in principle possible to reconstruct
the real-time scattering amplitude and the GPDs from our Euclidean Compton
amplitude. However, such a reconstruction is hampered by the fact that Eqs.
(6) and (8) have the form of a Fredholm integral equation of the first kind,
which is an ill-conditioned inverse problem known to have numerically unstable
solutions Horsley  _et al._ (2020).

As such, in this work we employ two strategies to access the structure
functions—similar techniques have been explored for the pseudo-distribution
method Karpie  _et al._ (2019); Joó  _et al._ (2019). First, we expand Eq. (6)
about $\bar{\omega}=0$, which gives us a power series in its Mellin moments.
By varying $\bar{\omega}$, we can then determine a finite set of the moments.
Such an extraction of structure function moments has been performed for the
forward Compton amplitude Chambers  _et al._ (2017); Can  _et al._ (2020);
Batelaan  _et al._ (2023) and in an exploratory calculation of the OFCA
Hannaford-Gunn  _et al._ (2022b).

While such a fit on the level of moments gives us a largely model-independent
determination, it is difficult to reconstruct the full structure function/GPD
from a limited number of moments. As such, for our second fit method we use a
phenomenologically-motivated GPD ansatz. In this work, we take

$\displaystyle H(x,t)\propto x^{-\alpha_{0}-\alpha^{\prime}t}(1-x)^{\beta}.$
(9)

Similar parameterisation have been widely used in both the phenomenological
studies of GPDs Goeke  _et al._ (2001); Belitsky  _et al._ (2004); Diehl  _et
al._ (2005); Guidal  _et al._ (2005); Jenkovszky (2006); Schoeffel (2007);
Kumerički  _et al._ (2008); Diehl and Kugler (2008); Gonzalez-Hernandez  _et
al._ (2013); Kumerički and Müller (2010); Kroll (2015) and fits to lattice
results from other methods Guo  _et al._ (2022, 2023). Our implementation of
this model-dependent fit follows previous work trialled in the forward case
Horsley  _et al._ (2020); Can  _et al._ (2023).

For both fitting approaches, we emphasise that the purpose is not to perform a
perturbative matching and extract a leading-twist GPD. Instead, our long-term
goal is a determination of the moments and real-time structure function
_including their power corrections and higher-twist effects_. Such corrections
provide us with useful information about the scale dependence of the
amplitude, which has not been studied experimentally in great detail for off-
forward Compton scattering. However, for the present work we perform
calculations for a single $\bar{Q}^{2}$ value to focus on determining the
structure functions and their moments.

Table 1:  Details of the gauge ensembles used in this work. $N_{f}$ |  $c_{SW}$ |  $\beta$ |  $\kappa_{l},\kappa_{s}$ |  $N_{L}^{3}\times N_{T}$ |  $a$ |  $m_{\pi}$ |  $Z_{V}$ |  $N_{\text{cfg}}$
---|---|---|---|---|---|---|---|---
|  |  |  |  |  [fm] |  [MeV] |  |
$2+1$ |  2.48 |  5.65 |  0.122005 |  $48^{3}\times 96$ |  0.068 |  $412$ |  0.871 |  537

##  III Determination of the off-forward Compton amplitude

To calculate the Compton amplitude in lattice QCD, we use the Feynman-Hellmann
(FH) method, which has proven to be a powerful tool to compute matrix elements
with one and two operator insertions, and a useful alternative to the direct
computation of three- and four-point functions.

This is implemented on the level of quark propagators, which are perturbed by
two background fields:

$\displaystyle
S_{(\lambda_{1},\lambda_{2})}(x_{n}-x_{m})=\big{[}M-\lambda_{1}\mathcal{O}_{1}-\lambda_{2}\mathcal{O}_{2}\big{]}^{-1}_{n,m},$
(10)

where $M$ is the Wilson fermion matrix and $\lambda_{1,2}$ are the FH
couplings.

The perturbing operators are

$\displaystyle\big{[}\mathcal{O}_{i}\big{]}_{n,m}=\delta_{n,m}(e^{i\mathbf{q}_{i}\cdot\mathbf{z}_{n}}+e^{-i\mathbf{q}_{i}\cdot\mathbf{z}_{n}})\bm{\gamma}\cdot\mathbf{\hat{e}}_{k},\;i=1,2,$
(11)

where ${\mathbf{q}}_{1,2}$ are the inserted momentum, and
$\hat{\mathbf{e}}_{k}$ picks the direction of the vector current.

From these perturbed two-point quark propagators, we then construct a
perturbed nucleon correlator:

$\begin{split}&\mathcal{G}^{\Gamma}_{\bm{\lambda}}(\tau,\mathbf{p}^{\prime})=\sum_{\mathbf{x}}e^{-i\mathbf{p}^{\prime}\cdot\mathbf{x}}\Gamma_{\beta\alpha}\prescript{}{\bm{\lambda}}{\langle}\Omega|\chi_{\alpha}(\mathbf{x},\tau)\chi^{\dagger}_{\beta}(0)|\Omega\rangle_{\bm{\lambda}},\end{split}$
(12)

where $\mathbf{p}^{\prime}$ is the sink momentum, $\tau$ is the Euclidean
time, $\Gamma$ is a spin-parity projector,
$\bm{\lambda}=(\lambda_{1},\lambda_{2})$ and $|\Omega\rangle_{\bm{\lambda}}$
the perturbed vacuum.

These perturbed nucleon propagators can be related to the OFCA by the Feynman-
Hellmann like relation Hannaford-Gunn  _et al._ (2022b)

$\displaystyle\begin{split}&R^{\Gamma}_{\lambda}(\mathbf{p}^{\prime},\tau)\overset{\tau\gg
a}{\simeq}\frac{\tau\lambda^{2}}{2E_{N}(\mathbf{p}^{\prime})}\mathcal{R}^{\Gamma}_{\mu\mu},\end{split}$
(13)

where $R^{\Gamma}_{\lambda}$ is the following combination of perturbed and
unperturbed nucleon correlators

$R^{\Gamma}_{\lambda}=\frac{\mathcal{G}^{\Gamma}_{(\lambda,\lambda)}+\mathcal{G}^{\Gamma}_{(-\lambda,-\lambda)}-\mathcal{G}^{\Gamma}_{(\lambda,-\lambda)}-\mathcal{G}^{\Gamma}_{(-\lambda,\lambda)}}{4\mathcal{G}^{\Gamma_{\text{unpol}}}_{(0,0)}},$
(14)

at some sink momentum $\mathbf{p}^{\prime}$ and sink time $\tau$, and
$\lambda$ is the magnitude of the Feynman-Hellmann coupling.

The quantity $\mathcal{R}^{\Gamma}_{\mu\nu}$ is defined as

$\displaystyle\mathcal{R}^{\Gamma}_{\mu\nu}=\frac{\sum_{s,s^{\prime}}\text{tr}\big{[}\Gamma
u(P^{\prime},s^{\prime})T_{\mu\nu}\bar{u}(P,s)\big{]}}{\sum_{s}{\text{tr}[\Gamma_{\text{unpol}}u(P^{\prime},s)\bar{u}(P^{\prime},s)]}},$
(15)

and $T_{\mu\nu}$ is the OFCA.

The spin-parity projectors we use are defined in Euclidean space as

$\Gamma_{\text{unpol}}=\frac{1}{2}(\mathbb{I}+\gamma_{4}),\quad\Gamma_{\text{pol}-\hat{\mathbf{e}}}=-\frac{i}{2}(\mathbb{I}+\gamma_{4})\hat{\mathbf{e}}\cdot\bm{\gamma}\gamma_{5}.$
(16)

We insert the perturbed quark propagators either for the doubly-represented
($u$) quarks, or the singly-represented ($d$) quarks individually. This means
that in practice we calculate

$\displaystyle T_{\mu\nu}^{ff}=i\int d^{4}ze^{i\bar{q}\cdot z}\langle
P^{\prime}|T\\{j^{f}_{\mu}(z/2)j^{f}_{\nu}(-z/2)\\}|P\rangle,$ (17)

where $j^{f}_{\mu}=Z_{V}\bar{\psi}_{f}\gamma_{\mu}\psi_{f}$, the local vector
current with flavour $f=u$ or $d$, and renormalisation factor $Z_{V}$. For
most analytic expressions we suppress flavour indices; however, we include
them for our numerical results.

The kinematics of the Compton amplitude are then completely derived from our
sink momentum, $\mathbf{p}^{\prime}$, and our two inserted momenta,
$\mathbf{q}_{1,2}$. For instance our momentum vectors, Eq. (2), become

$\displaystyle\bar{\mathbf{p}}=\frac{1}{2}(\mathbf{p}^{\prime}+\mathbf{p}),\quad\bar{\mathbf{q}}=\frac{1}{2}(\mathbf{q}_{1}+\mathbf{q}_{2}),\quad\mathbf{\Delta}=\mathbf{q}_{1}-\mathbf{q}_{2},$
(18)

where $\mathbf{p}=\mathbf{p}^{\prime}-\mathbf{\Delta}$ is the source momentum.
From these, the momentum scalars in Eq. (3) become

$\displaystyle\begin{split}\bar{\omega}=\frac{2\bar{\mathbf{p}}\cdot\bar{\mathbf{q}}}{\bar{\mathbf{q}}^{2}},\quad\bar{Q}^{2}=\bar{\mathbf{q}}^{2},\quad
t=-\mathbf{\Delta}^{2}.\end{split}$ (19)

To ensure the skewness variable is zero, we keep
$\mathbf{\Delta}\cdot\bar{\mathbf{q}}=0$, or equivalently
$|{\mathbf{q}}_{1}|=|{\mathbf{q}}_{2}|$.

The key kinematic choice made in this work in contrast to Ref. Hannaford-Gunn
_et al._ (2022b) is to keep
$\mathbf{\hat{e}}_{k}\propto{\mathbf{q}}_{1}-{\mathbf{q}}_{2}\equiv\mathbf{\Delta}$,
where $\mathbf{\hat{e}}_{k}$ picks the current direction in Eq. (11). Note
that $\mathbf{\hat{e}}_{k}$ is always chosen to be a purely spatial vector.
This means our vector current, $j_{\mu}$, is collinear with the soft momentum
transfer, $\Delta_{\mu}$.

Table 2: Current insertion momenta, $\mathbf{q}_{1,2}$, and derived kinematics for the four sets of correlators. $\frac{L}{2\pi}\mathbf{q}_{1}$, $\frac{L}{2\pi}\mathbf{q}_{2}$ |  $\frac{L}{2\pi}\mathbf{\Delta}$ |  $\frac{L}{2\pi}\bar{\mathbf{q}}$ |  $t$
$[\text{GeV}^{2}]$ |  $\bar{Q}^{2}$ $[\text{GeV}^{2}]$ |  $N_{\text{meas}}$
---|---|---|---|---|---
$(5,3,0)$ |  — |  — |  $0$ |  $4.86$ |  $1605$
$(4,3,3)$ $(3,4,3)$ |  $(1,-1,0)$ |  $(\frac{7}{2},\frac{7}{2},3)$ |  $-0.29$ |  $4.79$ |  $1031$
$(5,3,1)$ $(5,3,-1)$ |  $(0,0,2)$ |  $(5,3,0)$ |  $-0.57$ |  $4.86$ |  $1072$
$(4,2,4)$ $(2,4,4)$ |  $(2,-2,0)$ |  $(3,3,4)$ |  $-1.14$ |  $4.86$ |  $1031$

As shown in Appendix A, the choice
$\mathbf{\hat{e}}_{k}\propto\mathbf{\Delta}$ eliminates all structure
functions except $\mathcal{H}_{1}$ and $\mathcal{E}_{1}$ from the $\mu=\nu=k$
OFCA:

$\displaystyle\begin{split}{T}_{kk}&=-\frac{1}{2\bar{P}\cdot\bar{q}}\Big{(}h\cdot\bar{q}\mathcal{H}_{1}+e\cdot\bar{q}\mathcal{E}_{1}\Big{)},\end{split}$
(20)

where $h_{\mu}$ and $e_{\mu}$ are the bilinears defined in Eq. (5); recall
that $k$ is always a spatial direction (i.e. it has no temporal component).
This is a major improvement on our previous work Hannaford-Gunn  _et al._
(2022b), where only a linear combination of the $\mathcal{H}_{1,2,3}$ and
$\mathcal{E}_{1,2}$ was accessible, and we were required to make leading-twist
approximations to isolate these.

Inserting Eq. (20) into the ratio of spin-parity traces given in Eq. (15), we
obtain

$\mathcal{R}^{\Gamma}_{kk}=N^{\mathcal{H}}_{\Gamma}\mathcal{H}_{1}+N^{\mathcal{E}}_{\Gamma}\mathcal{E}_{1},$
(21)

where the $N$ factors are from the bilinear coefficients of Eq. (20) inserted
into the traces of Eq. (15).

Therefore, by using the two spin-parity projectors of Eq. (16), we determine
$\mathcal{H}_{1}$ and $\mathcal{E}_{1}$ from the matrix equation:

$\left(\begin{array}[]{c}\mathcal{R}^{\text{unpol}}_{kk}\\\
\mathcal{R}^{\text{pol}}_{kk}\end{array}\right)=\left(\begin{array}[]{cc}N^{\mathcal{H}}_{\text{unpol}}&N^{\mathcal{E}}_{\text{unpol}}\\\
N^{\mathcal{H}}_{\text{pol}}&N^{\mathcal{E}}_{\text{pol}}\end{array}\right)\left(\begin{array}[]{c}\mathcal{H}_{1}\\\
\mathcal{E}_{1}\end{array}\right).$ (22)

The matrix of $N$ factors has a determinant $\sim 1$ for all our kinematics,
making the inversion practical in all cases.

Figure 2: Top: the ratio defined in Eq. (15) for unpolarised and $y$-direction
polarised spin-parity projectors with the $\bar{\omega}=0$ term subtracted.
Bottom: The resulting subtracted off-forward structure functions,
$\overline{\mathcal{H}}_{1}$ and $\overline{\mathcal{E}}_{1}$. All results are
for $t=-0.57\;\text{GeV}^{2}$ and $uu-dd$ quarks.

### Calculation details

We perform our calculation of the off-forward Compton amplitude on a single
set of gauge fields generated by QCDSF/UKQCD Bietenholz  _et al._ (2011) with
$2+1$ quark flavours at the SU(3) flavour symmetric point, which yields an
unphysical pion mass of $m_{\pi}=412\;\text{MeV}$. See Table 1 for further
details.

We determine four sets of perturbed propagators, corresponding to four sets of
kinematics for the Compton amplitude with soft momentum transfers
$t=0,-0.29,-0.57,-1.14\;\text{GeV}^{2}$ all with a hard momentum transfer of
$\bar{Q}^{2}\approx 5\;\text{GeV}^{2}$ (see Table 2 for details). Furthermore,
we determine the propagators for two magnitudes of the Feynman-Hellmann
coupling $\lambda$: for $t=0,-0.29\;\text{GeV}^{2}$ we use
$\lambda=(0.0125,0.025)$, and for $t=-0.57,-1.14\;\text{GeV}^{2}$ we use
$\lambda=(0.00625,0.0125)$.

Our determination of the Compton amplitude from these correlators is similar
to that from Ref. Hannaford-Gunn  _et al._ (2022b): we fit our correlator
ratio in Eq. (13) as a linear function in Euclidean time, using a similar
weighted averaging method as in Ref. Beane  _et al._ (2021). Subsequently, we
fit this as a quadratic in $\lambda$—see Appendix B for details.

For each set of $\mathbf{q}_{1,2}$ we require a new set of inversions—see Eq.
(19). However, given a pair of $\mathbf{q}_{1,2}$, the $\bar{\omega}$ variable
can be expressed as

$\displaystyle\bar{\omega}=\frac{4\mathbf{p}^{\prime}\cdot(\mathbf{q}_{1}+\mathbf{q}_{2})}{(\mathbf{q}_{1}+\mathbf{q}_{2})^{2}}$
(23)

for $\mathbf{\Delta}\cdot\bar{\mathbf{q}}=0$. As such, by varying the sink
momentum, $\mathbf{p}^{\prime}$, we can obtain results at multiple values of
$\bar{\omega}$ for a single set of $t$ and $\bar{Q}^{2}$ values. See Appendix
C, Tab. 6 for all $\bar{\omega}$ values used in this work.

In Fig. 2, we plot the $\bar{\omega}$ dependence of the
$t=-0.57\;\text{GeV}^{2}$ results, and illustrate the determination of the
$\mathcal{H}_{1}$ and $\mathcal{E}_{1}$ structure functions from the spin-
parity traced quantity $\mathcal{R}$, given in Eq. (15). We observe good
signal for $\mathcal{R}$ for both spin-parity projectors, as well as the
subtracted Compton structure functions $\overline{\mathcal{H}}_{1}$ and
$\overline{\mathcal{E}}_{1}$. Moreover, these quantities all have the
$\bar{\omega}^{2}$ polynomial behaviour expected from Eq. (6). Results for the
$uu$ quark $\overline{\mathcal{H}}_{1}$ structure function across all $t$
values are presented in Fig. 3; see Fig. 11 in Appendix D for $dd$ quark
results.

##  IV Mellin Moment Fit

Our first fit strategy is to determine the Mellin moments of the real-time
scattering amplitudes. To this end we Taylor expand the dispersion relation
Eq. (6) about $\bar{\omega}=0$:

$\displaystyle\begin{split}\overline{\mathcal{H}}_{1}(\bar{\omega},t,\bar{Q}^{2})&=2\sum_{n=1}^{\infty}\bar{\omega}^{2n}\mathcal{A}_{2n,0}(t,\bar{Q}^{2}),\\\
\overline{\mathcal{E}}_{1}(\bar{\omega},t,\bar{Q}^{2})&=2\sum_{n=1}^{\infty}\bar{\omega}^{2n}\mathcal{B}_{2n,0}(t,\bar{Q}^{2}),\end{split}$
(24)

where $\mathcal{A}_{2n,0}$ and $\mathcal{B}_{2n,0}$ are the $n^{\text{th}}$
Mellin moments of the zero-skewness structure functions:

$\displaystyle\begin{split}\mathcal{A}_{2n,0}(t,\bar{Q}^{2})&=\frac{2}{\pi}\int_{0}^{1}dxx^{n-1}\text{Im}\mathcal{H}_{1}(x,t,\bar{Q}^{2}),\\\
\mathcal{B}_{2n,0}(t,\bar{Q}^{2})&=\frac{2}{\pi}\int_{0}^{1}dxx^{n-1}\text{Im}\mathcal{E}_{1}(x,t,\bar{Q}^{2}).\end{split}$
(25)

In a twist expansion, these moments can be expressed as Hannaford-Gunn  _et
al._ (2022b)

$\displaystyle\begin{split}\mathcal{A}_{2n,0}(t,\bar{Q}^{2})&=\mathcal{C}\left(\alpha_{S}(\bar{Q})\right)A_{2n,0}(t)+\text{higher
twist},\\\
\mathcal{B}_{2n,0}(t,\bar{Q}^{2})&=\mathcal{C}\left(\alpha_{S}(\bar{Q})\right)B_{2n,0}(t)+\text{higher
twist},\end{split}$ (26)

where $\mathcal{C}$ is the Wilson coefficient and $A_{n,0}$ and $B_{n,0}$ are
the generalised form factors (GFFs), given by the moments of the GPDs $H(x,t)$
and $E(x,t)$, respectively.

We emphasise that Eq. (26) is included to illustrate the connection between
the structure function moments ($\mathcal{A}_{n,0}$ and $\mathcal{B}_{n,0}$)
and the leading-twist GFFs (${A}_{n,0}$ and ${B}_{n,0}$). Although possible,
in this work we do not attempt a perturbative matching—our structure function
moments contain all the relevant non-leading-twist contributions.

### Fit details

From Eq. (24), we can fit the leading $N_{\text{max}}$ moments of the
subtracted structure functions to a power series in $\bar{\omega}$:

$\displaystyle\begin{split}&f(\bar{\omega})=2\sum_{n=1}^{N_{\text{max}}}\bar{\omega}^{2n}M_{2n},\end{split}$
(27)

where $M_{2n}$ are the moments of either $\overline{\mathcal{H}}_{1}$ or
$\overline{\mathcal{E}}_{1}$ as per Eq. (25).

For the $t=-0.29\;\text{GeV}^{2}$ results we note we can not access the
$\bar{\omega}=0$ point; see Appendix C for an explanation. As such, we fit the
$\bar{\omega}=0$ value simultaneously with our moments to the _unsubtracted_
structure functions:

$\displaystyle\begin{split}&f(\bar{\omega})=f(0)+2\sum_{n=1}^{N_{\text{max}}}\bar{\omega}^{2n}M_{2n}.\end{split}$
(28)

Figure 3: The subtracted Compton structure function
$\overline{\mathcal{H}}_{1}(\bar{\omega},t)$ for all $t$ values; $uu$ results
only. Note that our fits only use points for which
$|\mathbf{p}^{\prime}|<1\;\text{GeV}$; the shaded points are those with a sink
momentum greater than this. Curves correspond to all the fits performed in
this work: the model-independent fit, Eq. (27) (‘moment’); and the model-
dependent fits, Eq. (39), using the three sets of priors in Tab. 4 (‘wide’,
‘medium’ and ‘thin’). Note: the upper limit on the $y$-axis is held fixed
between the panels to demonstrate the change in magnitude with $-t$.

This allows us to determine the leading moments of each of the Compton
structure functions $\mathcal{H}_{1}$ and $\mathcal{E}_{1}$, which at leading-
twist are the GFFs $A_{2n,0}$ and $B_{2n,0}$, respectively.

We use the Bayesian Markov chain Monte Carlo (MCMC) package PyMC Salvatier
_et al._ (2016); Hoffman and Gelman (2014) to perform this fit. This allows us
to sample the model parameters from prior distributions that reflect physical
constraints.

For the $t=0$ kinematics our results are simply the forward Compton structure
function, $\mathcal{F}_{1}$, which is positive definite and hence has
monotonically decreasing Mellin moments Can  _et al._ (2020):

$\displaystyle 0\leq\mathcal{A}_{n+1,0}(t=0)\leq\mathcal{A}_{n,0}(t=0),$ (29)

where in this section we suppress the $\bar{Q}^{2}$ argument of the moments
for convenience.

Hence for the $n^{\text{th}}$ moment, we use a uniform prior distribution in
the range $\mathcal{A}_{n,0}(0)\in[0,\mathcal{A}_{n-2,0}(0)]$. We choose
$\mathcal{A}_{2,0}(0)\in[0,1]$ for the prior on the leading moment.

For the off-forward results, we no longer have monotonicity, so we use
positivity constraints on the GPDs Pobylitsa (2002), which at $\xi=0$ are:

$\displaystyle|H(x,t)|\leq
q(x),\quad|E(x,t)|\leq\frac{2m_{N}}{\sqrt{-t}}q(x),$ (30)

where $q(x)$ is the leading-twist parton distribution function. From these, it
is simple to determine the bounds on the moments at $\xi=0$:

$\displaystyle\big{|}A_{2n,0}(t)\big{|}\leq
a_{2n},\quad\big{|}B_{2n,0}(t)\big{|}\leq\frac{2m_{N}}{\sqrt{-t}}a_{2n},$ (31)

where $a_{n}$ is the $n^{\text{th}}$ parton distribution function moment.

Although Eq. (31) is derived for leading-twist GPDs we use it nonetheless,
noting that we are at a reasonably large $\bar{Q}^{2}$ and that these bounds
are not overly strict. As such, we adapt Eq. (31) to the prior:

$\displaystyle\big{|}\mathcal{A}_{2n,0}(t)\big{|}\leq\mathcal{A}_{2n,0}(0),\quad\big{|}\mathcal{B}_{2n,0}(t)\big{|}\leq\frac{2m_{N}}{\sqrt{-t}}\mathcal{A}_{2n,0}(0).$
(32)

For the $\mathcal{A}_{n,0}(0)$ bounds we use the mean plus one standard
deviation of the moments calculated from the $t=0$ results.

This fit is performed individually for each $t$ value across the
$\bar{\omega}$ values given in Tab. 6 of Appendix C. Note that in our fit we
only use $\bar{\omega}$ values for which the sink momentum is

$\displaystyle|\mathbf{p}^{\prime}|<1\;\text{GeV},$ (33)

as these are the points for which (1) we can better insure ground state
isolation, and (2) $\mathcal{O}(ap_{\mu})$ discretisation artefacts are
expected to be negligible for these points. We discuss these systematics
further in the next section.

### Results

As we can see in Fig. 3, the moment fit (labelled ‘moment’) describes the
$\bar{\omega}$-dependence of $\overline{\mathcal{H}}_{1}$ data well, with most
points being consistent within a standard deviation of the fits. Moreover, the
fits are even consistent with some of the $|\mathbf{p}^{\prime}|\geq
1\;\text{GeV}$ points that were not included in the fits. We note that for
points with $|\mathbf{p}^{\prime}|<1\;\text{GeV}$, the
$t=0,-0.29,-0.57\;\text{GeV}^{2}$ results do not have $\bar{\omega}$ points
greater than $\bar{\omega}\approx 0.8$, while for $t=-1.14\;\text{GeV}^{2}$,
there are no such points beyond $\bar{\omega}\approx 0.6$. This limits our
ability to constrain the higher moments in the present work. Moreover, since
the $t=-0.29\;\text{GeV}^{2}$ results require us to fit the $\bar{\omega}=0$
subtraction function, instead of determining it directly, these fits are
generally of a poorer quality. Note in Fig. 3 for the
$t=-0.29\;\text{GeV}^{2}$ we subtract off the fitted $\bar{\omega}=0$ point.

See Appendix E for the posterior distributions from the Bayesian Markov chain
Monte Carlo fits. Note that while using the fit function in Eq. (27), it is
necessary to choose the number of moments to fit, $N_{\text{max}}$. To make
this choice, we compare the effect of varying this parameter in Fig. 13. For
$N_{\text{max}}\geq 3$ the order of truncation has negligible impact on the
leading moments, and therefore we take $N_{\text{max}}=7$.

Figure 4: The $n=2$ off-forward Compton amplitude moments $\mathcal{A}_{2n,0}$
(top) and $\mathcal{B}_{2n,0}$ (bottom) determined from a moment fit. These
are the moments of $\mathcal{H}_{1}$ and $\mathcal{E}_{1}$ as in Eq. (25). The
Compton amplitude results are for the $uu-dd$ quarks. We compare this to
$A_{2,0}$ and $B_{2,0}$ GFFs for $u-d$ quark combination, calculated from the
local twist-two lattice operators on the same set of gauge configurations.
Note the twist-two operators are renormalised using
$\text{RI}^{\prime}\text{/MOM}$., while the OFCA amplitude moments contain all
power-corrections and higher-twist contributions.

In addition to determining our OFCA moments, we also determine the generalised
form factors $A_{n,0}$ and $B_{n,0}$ for $n=1,2$ using the local twist-two
operators on the same set of gauge configurations. The matrix elements of
these local operators are computed using standard three-point function
methodsGöckeler  _et al._ (2004)—see Appendix F for the details. As per Eq.
(26), the off-forward structure function moments $\mathcal{A}_{2n,0}$ and
$\mathcal{B}_{2n,0}$ correspond to the GFFs $A_{2n,0}$ and $B_{2n,0}$,
respectively, up to higher-twist, power-corrections and the Wilson
coefficient. As such, the $A_{2,0}$ and $B_{2,0}$ GFFs determined from local
operators are a useful point of comparison.

Finally, we fit $n=2$ moments as a function of $t$, using the dipole
parameterisation:

$\displaystyle G(t)=\frac{G(0)}{\left(1-t/m_{\text{dip}}^{2}\right)^{2}}.$
(34)

Given the large uncertainties on our points, we only use this simple
parameterisation and do not test the effects of different parameterisations
for form factor fits. See Tab. 3 for the parameters of our dipole fits for
$uu-dd$ quarks.

Table 3: Summary of parameters from dipole fit. All results for $uu-dd$ quarks. | $G(0)$ | $m_{\text{dip}}$
---|---|---
$\mathcal{A}_{2,0}$ | 0.226(59) | 1.8(1.1)
$\mathcal{B}_{2,0}$ | 0.50(26) | 1.8(1.3)

In Fig. 4 we plot the OFCA moments $\mathcal{A}_{2,0}$ and $\mathcal{B}_{2,0}$
as functions of the soft momentum transfer $t$, including a comparison to
${A}_{2,0}$ and ${B}_{2,0}$ GFFs. We observe good agreement between the
helicity-conserving moments $\mathcal{A}_{2,0}$ and the GFF $A_{2,0}$ across
the range of $t$ values. Similarly, there is reasonable agreement between the
helicity-flipping moments $\mathcal{B}_{2,0}$ and $B_{2,0}$.

However, we emphasise that, even with complete control of lattice systematics,
our structure function moments should be distinct from the leading-twist GFFs,
and as such we do not attempt a strong comparison between these results. Nor
do we attempt a separation of the leading-twist contributions and the power-
corrections; such a determination has been performed on our forward Compton
amplitude results where more $Q^{2}$ values are available Batelaan  _et al._
(2023). An equivalent study of the $\bar{Q}^{2}$ dependence of the off-forward
Compton amplitude could provide useful information on the non-leading-twist
contributions to these moments. It is, however, encouraging that there is
reasonable agreement between the $t$ dependence of the Compton amplitude
moments and that of the local twist-two operators.

Finally, we note that the parameters from our dipole fits broadly agree with
other fits to generalised form factors calculated from local twist-two
operators at similar pion masses Hägler  _et al._ (2008). Moreover, we
determine the $u-d$ quark angular momentum from the Ji sum rule Ji (1997):

$\displaystyle\langle J_{u-d}^{3}\rangle\approx\frac{1}{2}[\mathcal{A}^{uu-
dd}_{2,0}(0,\bar{Q}^{2})+\mathcal{B}^{uu-dd}_{2,0}(0,\bar{Q}^{2})]=0.36(16),$
(35)

neglecting the non-leading-twist corrections to our moments. Again, this
agrees with determinations from local twist-two operators at similar pion
masses Hägler  _et al._ (2008), although our errors are very large mostly
owing to the statistical uncertainties on the
$\mathcal{B}_{2,0}(0,\bar{Q}^{2})$ dipole fit.

Despite the size of these uncertaintes, this calculation provides an
alternative means of determining the Ji sum rule. Moreover, determinations of
the OFCA with multiple $\bar{Q}^{2}$ values would allow us to analyse the hard
scale dependence of this quantity, which is not achievable from other methods.

##  V Model Fit

In the previous section, we determined the Mellin moments of the real-time
off-forward structure functions from our Euclidean OFCA. While this
determination is largely model-independent, it is difficult to reconstruct the
complete real-time structure functions (and hence GPDs) from a limited set of
Mellin moments.

As such, for our second fit strategy we use the phenomenological
parameterisation of a GPD (or off-forward structure function),

$\displaystyle H(x,t)=Cx^{-\alpha(t)}(1-x)^{\beta},$ (36)

with $\alpha(t)=\alpha_{0}+\alpha^{\prime}t$, where $\alpha^{\prime}$ is the
Regge slope parameter.

Note that this parameterisation is normalised by the factor

$\displaystyle
C=A\int_{0}^{1}dxx^{-\alpha_{0}}(1-x)^{\beta}=A\frac{\Gamma(3-\alpha_{0}+\beta)}{\Gamma(2-\alpha_{0})\Gamma(\beta+1)},$
(37)

which ensures that $\mathcal{A}_{2,0}(t=0,\bar{Q}^{2})=A$. This gives us a
total of four parameters in our model: $A,\alpha_{0},\alpha^{\prime}$ and
$\beta$. We then perform a global fit with this parameterisation to our
Compton amplitude results for all $t$ values.

The model in Eq. (36) and similar Regge-inspired parameterisations of GPDs
have been used widely to determine GPD properties from various experimental
processes Goeke  _et al._ (2001); Belitsky  _et al._ (2004); Diehl  _et al._
(2005); Guidal  _et al._ (2005); Jenkovszky (2006); Schoeffel (2007);
Kumerički  _et al._ (2008); Diehl and Kugler (2008); Gonzalez-Hernandez  _et
al._ (2013); Kumerički and Müller (2010); Kroll (2015), as well as in fits to
other lattice results Guo  _et al._ (2022, 2023).

Inserting the ansatz in Eq. (36) into our dispersion relation, Eq. (8), we
obtain

$\displaystyle\begin{split}&\overline{\mathcal{H}}_{1}(\bar{\omega},t)=2C\bar{\omega}^{2}\frac{\Gamma(2-\alpha(t))\Gamma(\beta+1)}{\Gamma(3+\beta-\alpha(t))}\\\
&\times\prescript{}{3}{F}_{2}\bigg{[}{1,(2-\alpha(t))/2,(3-\alpha(t))/2\atop(3+\beta-\alpha(t))/2,(4+\beta-\alpha(t))/2};\bar{\omega}^{2}\bigg{]},\end{split}$
(38)

where $\prescript{}{3}{F}_{2}$ is a generalised hypergeometric function.

Figure 5: Posterior distributions of our three model parameters for $u$
quarks; note that $\alpha_{0}$ is not a model parameter but its posterior can
be reconstructed with Eq. (42).

Equation (38) can be expressed as the sum of moments:

$\displaystyle\overline{\mathcal{H}}_{1}(\bar{\omega},t)=2C\sum_{n=1}^{\infty}\bar{\omega}^{2n}\frac{\Gamma(2n-\alpha(t))\Gamma(\beta+1)}{\Gamma(1+2n-\alpha(t)+\beta)},$
(39)

with the $n^{\text{th}}$ moment as

$\displaystyle\mathcal{A}_{n,0}(t)=C\frac{\Gamma(n-\alpha(t))\Gamma(\beta+1)}{\Gamma(1+n-\alpha(t)+\beta)},$
(40)

which is similar to Regge-inspired models of elastic form factors Di Vecchia
and Drago (1969).

To simplify the implementation of the fit, we use Eq. (39) as our fit
function, truncating at a very high order, $n=50$, which ensures even marginal
effects from the higher moments are negligible.

We note that this model is best justified for valence quarks, though our
results include sea quark contributions (i.e. $H(-x,t)$) that we take to be
suppressed. The assumption that our distributions are dominated by valence
quark contributions allows us to make the further constraint on our parameters
that

$\displaystyle\int_{0}^{1}dxH^{q}(x,t=0)=F^{q}_{1}(|t|=0)=N_{q},$ (41)

for $N_{q}$ the number of valence quarks of flavour $q$. While Eq. (41) is
strictly only true for leading-twist, it should nonetheless be a good
approximation for $\bar{Q}^{2}\approx 5\;\text{GeV}^{2}$.

Figure 6: The Compton amplitude moments $\mathcal{A}_{n,0}$, determined from
the model-dependent fit, Eq. (39) with the three sets of priors (‘wide’,
‘medium’ and ‘thin’). We compare these to the moments from the direct moment
fit, Eq. (27) (‘moment’); all results for $u$ quarks. In addition we compare
our Compton amplitude moments to the GFFs calculated with local twist-two
operators on the same gauge configurations (‘$A_{n,0}$ 3-pt’). We emphasise
again that the Compton amplitude moments, $\mathcal{A}_{n,0}$ are distinct
from the leading-twist GFFs $A_{n,0}$, although they can be related via Eq.
(26).

Equation (41) gives us

$\displaystyle\alpha_{0}=\frac{N_{q}-A(2+\beta)}{N_{q}-A},$ (42)

which we use to remove the parameter $\alpha_{0}$ from our fits.

### Fit details

We fit the model in Eq. (39) simultaneously to all our soft momentum transfer
values: $t=0,-0.29,-0.57,-1.14\;\text{GeV}^{2}$. Note that we only fit the
$\mathcal{H}_{1}$ structure function as the $\mathcal{E}_{1}$ results are
typically poorer quality and lack the $t=0$ Compton amplitude.

The fit is again performed with Bayesian MCMC. However, in contrast to the
moment fits, there are no model-independent priors for these parameters. As
such, to test the dependence of this fit on our priors we vary the width of
the priors around approximate phenomenologically expected values:
$\alpha^{\prime}=0.9\;\text{GeV}^{-2}$ Goeke  _et al._ (2001); Belitsky  _et
al._ (2004); Diehl  _et al._ (2005) and $\beta=3$, while keeping $A\in[0,1]$
for all fits. We fit to our $uu$ and $dd$ structure functions separately, with
three parameters for each flavour. All prior distributions are uniform. See
Tab. 4 for a summary of the three priors.

As with the Mellin moment fits, we only use $\bar{\omega}$ values for which
the sink momentum is $|\mathbf{p}^{\prime}|<1\;\text{GeV}$.

Table 4: Priors for the three parameters. | $A$ |  $\alpha^{\prime}$ [$\text{GeV}^{-2}$] | $\beta$
---|---|---|---
wide | $[0,1]$ | $[0.0,2.5]$ | $[0,8]$
medium | $[0,1]$ | $[0.2,1.7]$ | $[1,6]$
thin | $[0,1]$ | $[0.4,1.4]$ | $[2,4]$

### Results

In Fig. 3, we plot the fits to the $uu$ quark Compton structure function
$\overline{\mathcal{H}}_{1}$ for all $t$ values; see Fig. 11 for $dd$ quark
results. Among the three sets of priors for model-dependent fits, we see good
agreement in $\bar{\omega}$-space apart from the regions that are not
constrained (i.e. where there are no $\bar{\omega}$ values in the fit).
Comparing the model-independent moment fit and our model fits, discrepancies
are apparent for $\bar{\omega}\gtrsim 0.3$.

Table 5: The values of the Regge slope parameter $\alpha^{\prime}$ for the two flavours and three sets of model priors. $\alpha^{\prime}$ $[\text{GeV}^{-2}]$ | $u$ | $d$
---|---|---
wide | $1.32(57)$ | $1.12(63)$
medium | $0.92(36)$ | $0.92(41)$
thin | $0.85(26)$ | $0.89(28)$

In Fig. 5, we plot the posterior distributions of the model parameters; see
Fig. 10 of Appendix D for $d$ quark results. We similarly observe good
agreement for the $A$ parameter among the three sets of priors, reasonable
agreement in $\alpha^{\prime}$, and a discrepancy among the three fits for the
$\beta$ parameters.

The discrepancies between the model fits for large $\bar{\omega}$ and the
model-dependence of $\beta$ parameter are related. As discussed in the
previous section, higher moments are constrained by the large $\bar{\omega}$
values of our Compton amplitude. In terms of our model parameterisation,
higher moments are more dependent on the $x\to 1$ behaviour of the GPD, which
is determine by $\beta$. As such, a unique determination of $\beta$ requires
accurate and precise determinations of the Compton amplitude for large
$\bar{\omega}$, which we currently do not have. We emphasise that the
determination of higher-quality large $\bar{\omega}$ values is simply a matter
of systematic improvement, not a fundamental limitation.

For the ‘medium’ and ‘thin’ priors, we find $\alpha^{\prime}$ values of
$\approx 0.9\;\text{GeV}^{-2}$ for both $u$ and $d$ quarks, while for the
‘wide’ priors the mean of the distributions is slightly larger, although still
consistent with the other two fits. This compares well with phenomenological
studies, which give the range $\alpha^{\prime}\approx
0.8-1.8\;\text{GeV}^{-2}$ for valence quarks found in both deeply virtual
Compton scattering analyses and fits to elastic form factors Goeke  _et al._
(2001); Belitsky  _et al._ (2004); Diehl  _et al._ (2005); Schoeffel (2007);
Gonzalez-Hernandez  _et al._ (2013); Kumerički and Müller (2010). Moreover,
the ‘thin’ and ‘medium’ results compare well with the value of
$\alpha^{\prime}=0.871(6)\;\text{GeV}^{-2}$ found in global fits to the light
hadron masses Jia  _et al._ (2017). See Tab. 5 for all $\alpha^{\prime}$
results.

Figure 7: The helicity-conserving GPD, $H(x,t)$, for $u$ quarks, determined
from fitting the model parameterisation Eq. (39) to our Compton amplitude
results. This uses the ‘medium’ priors given in Tab. 4.

In Fig. 6, we show the $u$ quark results for $\mathcal{A}_{n,0}$, the Mellin
moments of $\mathcal{H}_{1}$ determined from this model fit; see Fig. 12 of
Appendix D for $d$ quark results. We also include the moments determined from
the model-independent moment fit, Eq. (27). Moreover, we compare these to the
Dirac elastic form factor (i.e. $A_{1,0}$) and $A_{2,0}$ as discussed in the
previous section, both determined on the same set of gauge configurations from
twist-two local operators—see Appendix F for details of this calculation.
Again, we strongly emphasise, as per Eq. (26), our $\mathcal{A}_{n,0}$ moments
and the $A_{n,0}$ GFFs are not necessarily equivalent.

We see very strong agreement between our ‘medium’ and ‘thin’ fits and the
three-point results for $A_{1,0}$. While this agreement is enforced at $t=0$
by the condition Eq. (41), the agreement in the $t$-dependence is nonetheless
promising. However, the ‘wide’ parameters show a markedly steeper drop off in
$t$. For $\mathcal{A}_{2,0}$ we see a stronger agreement among the three
priors for the model fits than for $\mathcal{A}_{1,0}$, and reasonable
agreement of these with both the direct moment fits and the three-point
results.

Figure 8: The impact parameter space distribution $xH(x,{b})$ for $d$ quarks,
determined from the model parameterisation fit using the ‘medium’ priors given
in Tab. 4. Note colours correspond to $z$-axis value, and are simply to help
convey the shape.

For the $\mathcal{A}_{4,0}$ results, there is slightly less agreement among
the three sets of priors in the model fit and the direct moment fits are less
well-constrained compared to $\mathcal{A}_{2,0}$. As can be seen in Fig. 3,
the direct moment fits are highly unconstrained for large $\bar{\omega}$,
especially for $t=-0.29,-1.14\;\text{GeV}^{2}$, which explains the
irregularity in these results.

In Fig. 7, we plot the $x$-space GPD for the $u$ quark and ‘medium’ priors.
While this result is roughly the expected shape, we note that the GPD appears
to have a somewhat slow drop off for $x\to 1$, and the peak of $xH(x,0)$
should be closer to $0.2$ than $0.4$, although the position of this peak is
dependent on the hard scale and pion mass. Such features are likely a result
of the $\beta$ values skewing small, which is equivalent in moment space to
the higher moments not falling quickly enough. As discussed previously, more
and better quality results for large $\bar{\omega}$ are necessary to better
constrain the higher moments and the $\beta$ parameter.

Finally, in Fig. 8 we plot the impact parameter space distribution $xH(x,{b})$
Burkardt (2000) obtained via a Fourier transformation
($\mathbf{\Delta}\to\mathbf{b}$, where $b=|\mathbf{b}|$) for the $d$ quark,
using the mean values of the ‘medium’ priors. The distribution $H(x,{b})$ can
be interpreted as the probability density for a $d$ quark in a fast-moving
nucleon as a function of the momentum fraction $x$ and a distance $b$ from the
nucleon’s centre of mass. Therefore, weighted by $x$, the distribution
$xH(x,{b})$ is the momentum density. At this stage given our control of
lattice systematics, such a plot is simply a demonstration that our method can
be used to reconstruct the impact parameter space distributions.

In general, we note that the GPD ansatz, Eq. (36), is successful in describing
our data and reproducing physically expected properties of GPDs, despite its
simplicity. Future studies with higher quality data could test a range of
different GPD models.

### Systematics and Future Improvements

The results presented here indicate a great deal of potential for lattice QCD
calculations of the off-forward Compton amplitude, as a means to provide
unique and complementary physical information about GPDs and off-forward
scattering. However, we have also revealed and clarified areas in need of
improvement. The most significant of these is the need for more accurate and
precise determinations of points with large sink momenta. As discussed, these
points give us access to large $\bar{\omega}$ values and hence are crucial in
constraining higher moments.

To improve the quality of the signal, there exist numerous methods to better
isolate the ground state for correlators with large sink momentum Egerer  _et
al._ (2021); Bali  _et al._ (2016); Wu  _et al._ (2018), which are already
used widely in the calculation of quasi and Ioffe time distributions. Such
methods are capable of determining the ground state in correlators with sink
momenta as large as $|\mathbf{p}^{\prime}|\approx 3\;\text{GeV}$, three times
greater than the largest sink momentum used in our fits.

In addition, we have treated our $\mathcal{H}_{1}$ and $\mathcal{E}_{1}$ data
as if they were continuum results, without attempting to account for any
discretisation artefacts. We are currently expanding on previous work of a
lattice operator product expansion of the Compton amplitude Göckeler  _et al._
(2006) that would allow us to control such artefacts.

Finally, the lattice systematics that are not specific to our method—the
unphysical quark masses, lattice spacing and volume—must be accounted for
before a strong comparison can be made with phenomenology.

##  VI Summary and Conclusions

In this study, we present a significantly improved determination of the
structure functions of the off-forward Compton amplitude (OFCA). In
particular, we determine the structure functions $\mathcal{H}_{1}$ and
$\mathcal{E}_{1}$ independently for a wider range of kinematics than our
previous study Hannaford-Gunn  _et al._ (2022b). This separation allows us to
perform a more in-depth attempt at determining properties of the real-time
structure functions despite the inverse problem.

We also calculate the $n=1,2$ generalised form factors (GFFs) from the local
twist-two lattice operators to compare to our Compton amplitude results.
Although we do not perform a perturbative matching of our Compton amplitude
moments, we nonetheless note reasonable agreement between the Compton
amplitude moments and twist-two GFFs. Similarly, our determinations of the
Regge slope parameter $\alpha^{\prime}$ broadly agree with those from fits to
experiment. These agreements are promising, and show that our method is
capable of determining meaningful physical information.

Our analysis also clarifies key systematics in need of addressing; in
particular, our determinations of the structure functions for large
$\bar{\omega}$, which requires large sink momenta for our correlators, appear
to suffer from lattice artefacts. As discussed in the previous section,
addressing these systematics is simply a matter of using and/or building upon
existing techniques, making a precise and accurate determination of the off-
forward Compton amplitude completely achievable.

An improved lattice QCD determination of the OFCA would provide a valuable
comparison for studies in the quasi- and pseudo-distribution formalisms.
Moreover, our method is unique in being able to determine non-leading-twist
effects, as has been done for the forward Compton amplitude Hannaford-Gunn
_et al._ (2020); Batelaan  _et al._ (2023). Such effects could provide useful
phenomenological information, as most experimental studies of hard exclusive
processes have a modest hard scale of $Q^{2}\approx 1-12\;\text{GeV}^{2}$ and
contain additional $|t|/{Q}^{2}$ corrections Braun  _et al._ (2012); Guo  _et
al._ (2021).

Moreover, past work to determine the subtraction function of the forward
Compton amplitude Hannaford-Gunn  _et al._ (2022a) could be extended to the
OFCA subtraction function. This off-forward subtraction function is a key
input for determinations of the proton pressure distribution Burkert  _et al._
(2018); Kumerički (2019), and as such could significantly reduce the errors of
model-independent measurements of this quantity.

## Acknowledgements

The numerical configuration generation (using the BQCD lattice QCD program
Haar  _et al._ (2018))) and data analysis (using the Chroma software library
Edwards and Joó| (2005)) was carried out on the DiRAC Blue Gene Q and Extreme
Scaling (EPCC, Edinburgh, UK) and Data Intensive (Cambridge, UK) services, the
GCS supercomputers JUQUEEN and JUWELS (NIC, Jülich, Germany) and resources
provided by HLRN (The North-German Supercomputer Alliance), the NCI National
Facility in Canberra, Australia (supported by the Australian Commonwealth
Government) and the Phoenix HPC service (University of Adelaide). AHG and JAC
are supported by an Australian Government Research Training Program (RTP)
Scholarship. RH is supported by STFC through grants ST/T000600/1 and
ST/X000494/1. PELR is supported in part by the STFC under contract
ST/G00062X/1. GS is supported by DFG Grant No. SCHI 179/8-1. KUC, RDY and JMZ
are supported by the Australian Research Council grants DP190100297 and
DP220103098.

##  Appendix A Isolating $\mathcal{H}_{1}$ and $\mathcal{E}_{1}$

We start with the tensor decomposition from Ref. Hannaford-Gunn  _et al._
(2022b), removing the terms that vanish for $\xi=0$ Tarrach (1975):

$\begin{split}&\bar{T}_{\mu\nu}=\frac{1}{2\bar{P}\cdot\bar{q}}\bigg{\\{}-\Big{(}h\cdot\bar{q}\mathcal{H}_{1}+e\cdot\bar{q}\mathcal{E}_{1}\Big{)}g_{\mu\nu}+\frac{1}{\bar{P}\cdot\bar{q}}\Big{(}h\cdot\bar{q}\mathcal{H}_{2}+e\cdot\bar{q}\mathcal{E}_{2}\Big{)}\bar{P}_{\mu}\bar{P}_{\nu}+\mathcal{H}_{3}h_{\\{\mu}\bar{P}_{\nu\\}}\bigg{\\}}\\\
&+\frac{i}{2\bar{P}\cdot\bar{q}}\epsilon_{\mu\nu\rho\kappa}\bar{q}^{\rho}\bigg{\\{}\tilde{h}^{\kappa}\tilde{\mathcal{H}}_{1}+\tilde{e}^{\kappa}\tilde{\mathcal{E}}_{1}+\frac{1}{\bar{P}\cdot\bar{q}}\Big{[}\big{(}\bar{P}\cdot\bar{q}\tilde{h}^{\kappa}-\tilde{h}\cdot\bar{q}\bar{P}^{\kappa}\big{)}\tilde{\mathcal{H}}_{2}\Big{]}\bigg{\\}}+\Big{(}\bar{P}_{\mu}q^{\prime}_{\nu}+\bar{P}_{\nu}q_{\mu}\Big{)}\Big{(}h\cdot\bar{q}\mathcal{K}_{1}+e\cdot\bar{q}\mathcal{K}_{2}\Big{)}\\\
&+q_{\mu}q^{\prime}_{\nu}\big{(}h\cdot\bar{q}-e\cdot\bar{q}\big{)}\mathcal{K}_{5}+\Big{(}h_{\mu}q^{\prime}_{\nu}+h_{\nu}q_{\mu}\Big{)}\mathcal{K}_{7},\end{split}$
(43)

where we note $\bar{T}_{\mu\nu}$ is the Compton amplitude without gauge
projection; that is, it contains no terms with uncontracted $q^{\prime}_{\mu}$
and $q_{\nu}$. Further note that $h_{\mu}$ and $e_{\mu}$ are given in Eq. (5)
and $\tilde{h}_{\mu}=\bar{u}(P^{\prime})\gamma_{\mu}\gamma_{5}u(P)$ and
$\tilde{e}_{\mu}=\Delta_{\mu}\bar{u}(P^{\prime})\gamma_{5}u(P)/2m_{N}$.

The key kinematic choice for this work in contrast to Ref. Hannaford-Gunn  _et
al._ (2022b), is that we take $\mathbf{\hat{e}}_{k}\propto\mathbf{\Delta}$,
where $\mathbf{\hat{e}}_{k}$ is the vector that picks the direction of the
current in Eq. (11). As both $\mathbf{\bar{q}}$ and $\mathbf{\bar{p}}$ are
orthogonal to $\mathbf{\Delta}$, the choice
$\mathbf{\hat{e}}_{k}\propto\mathbf{\Delta}$ means that any terms in our
tensor decomposition with an uncontracted $\bar{q}_{\mu}$ or $\bar{P}_{\mu}$
vanish. Therefore, only tensor structures with $g_{\mu\nu}$ or
$\Delta_{\mu}\Delta_{\nu}$ survive. The former are associated with
$\mathcal{H}_{1}$ and $\mathcal{E}_{1}$ amplitude, while the latter are
suppressed.

As such, this kinematic choice minimises the effects of EM gauge dependent
terms (i.e. any terms with uncontracted $q^{\prime}_{\mu}$ and $q_{\nu}$) and
hence discretisation artefacts, as our local current does not satisfy the
continuum Ward identities Karsten and Smith (1981); Guerin (1987). Moreover,
it helps us isolate $\mathcal{H}_{1}$ and $\mathcal{E}_{1}$ instead of a
linear combination of other structure functions.

Explicitly, this choice means that the polarised structure functions,
$\tilde{\mathcal{H}}_{1,2}$ and $\tilde{\mathcal{E}}_{1,2}$, are attached to

$\Delta_{\\{\mu}\epsilon_{\nu\\}\sigma\rho\kappa}\Delta^{\sigma}\bar{q}^{\rho}\tilde{h}^{\kappa},$

after gauge projection. Since $\mu,\nu$ are the Compton amplitude’s indices,
with the $\mathbf{\hat{e}}_{k}\propto\mathbf{\Delta}$ condition we have
$\mu=\nu=\sigma$, and hence the above equation must vanish. This completely
removes all polarised amplitudes.

Further, the $\mathcal{K}_{1,2,5,7}$ amplitudes have no $g_{\mu\nu}$ tensor
structure. Therefore, after gauge projection, the only contribution that
survives is

$\frac{\Delta_{k}\Delta_{k}}{q\cdot
q^{\prime}}\mathcal{K}_{1,2,5,7}\sim\frac{-t}{\bar{Q}^{2}}\mathcal{K}_{1,2,5,7}.$

Hence these tensor structures, which are already non-leading-twist, receive an
additional kinematic suppression of $|t|/\bar{Q}^{2}$.

Therefore, up to highly suppressed terms containing the $\mathcal{K}$
amplitudes, the gauge-projected Compton amplitude is

$\begin{split}&{T}_{kk}=\frac{1}{2\bar{P}\cdot\bar{q}}\bigg{[}\Big{(}h\cdot\bar{q}\mathcal{H}_{1}+e\cdot\bar{q}\mathcal{E}_{1}\Big{)}\bigg{(}1+\frac{q^{\prime}_{k}q_{k}}{q\cdot
q^{\prime}}\bigg{)}+\frac{1}{\bar{P}\cdot\bar{q}}\Big{(}h\cdot\bar{q}\mathcal{H}_{2}+e\cdot\bar{q}\mathcal{E}_{2}\Big{)}\\\
&\times\bigg{(}\bar{P}_{k}\bar{P}_{k}-\frac{\bar{P}\cdot\bar{q}}{q\cdot
q^{\prime}}(q^{\prime}_{k}\bar{P}_{k}+\bar{P}_{k}q_{k})+\bigg{(}\frac{\bar{P}\cdot\bar{q}}{q\cdot
q^{\prime}}\bigg{)}^{2}q^{\prime}_{k}q_{k}\bigg{)}+\mathcal{H}_{3}\bigg{(}\bar{P}_{k}h_{k}-\frac{\bar{P}\cdot\bar{q}}{q\cdot
q^{\prime}}(q^{\prime}_{k}h_{k}+h_{k}q_{k})+\frac{\bar{P}\cdot\bar{q}h\cdot\bar{q}}{(q\cdot
q^{\prime})^{2}}q^{\prime}_{k}q_{k}\bigg{)}\bigg{]}.\end{split}$ (44)

The Dirac bilinear $h_{\mu}$ is orthogonal to $\Delta_{\mu}$, so that the
$h_{k}$ terms, which are proportional to $h\cdot\Delta$, must vanish.
Moreover, as previously explained, $\bar{P}_{k}=0=\bar{q}_{k}$. Therefore, Eq.
(44) becomes

$\begin{split}{T}_{kk}&=\frac{1}{2\bar{P}\cdot\bar{q}}\bigg{[}\Big{(}h\cdot\bar{q}\mathcal{H}_{1}+e\cdot\bar{q}\mathcal{E}_{1}\Big{)}\bigg{(}1-\frac{\Delta_{k}\Delta_{k}}{4q\cdot
q^{\prime}}\bigg{)}-\frac{\bar{P}\cdot\bar{q}}{4(q\cdot
q^{\prime})^{2}}\Big{(}h\cdot\bar{q}(\mathcal{H}_{2}+\mathcal{H}_{3})+e\cdot\bar{q}\mathcal{E}_{2}\Big{)}{\Delta_{k}\Delta_{k}}\bigg{]}.\end{split}$
(45)

In Ref. Hannaford-Gunn  _et al._ (2022b) it was shown that for large
$\bar{Q}^{2}$ the $\mathcal{H}$ and $\mathcal{E}$ structure functions satisfy
the off-forward Callan-Gross relation:

$\frac{\bar{\omega}}{2}\big{(}\mathcal{H}_{2}+\mathcal{H}_{3}\big{)}=\mathcal{H}_{1}+\Delta\mathcal{H}_{\text{CG}},\quad\frac{\bar{\omega}}{2}\mathcal{E}_{2}=\mathcal{E}_{1}+\Delta\mathcal{E}_{\text{CG}},$
(46)

where we have included $\Delta\mathcal{H}_{\text{CG}}$ and
$\Delta\mathcal{E}_{\text{CG}}$, the $\mathcal{O}(\alpha_{S})$ violations to
this Callan-Gross relation.

Further, note that $q\cdot q^{\prime}=-\bar{Q}^{2}+t/4$. Hence Eq. (45)
becomes

$\begin{split}{T}_{kk}&=\frac{1}{2\bar{P}\cdot\bar{q}}\bigg{\\{}\Big{(}h\cdot\bar{q}\mathcal{H}_{1}+e\cdot\bar{q}\mathcal{E}_{1}\Big{)}\bigg{(}1-\frac{\Delta_{k}\Delta_{k}}{4(\bar{Q}^{2}-t/4)}\frac{t}{4(\bar{Q}^{2}-t/4)}\bigg{)}\\\
&-\Big{(}h\cdot\bar{q}\Delta\mathcal{H}_{\text{CG}}+e\cdot\bar{q}\Delta\mathcal{E}_{\text{CG}}\Big{)}\frac{\Delta_{k}\Delta_{k}}{4(\bar{Q}^{2}-t/4)}\frac{\bar{Q}^{2}}{\bar{Q}^{2}-t/4}\bigg{\\}}.\end{split}$
(47)

Given that $\Delta\mathcal{H}_{\text{CG}}$ and $\Delta\mathcal{E}_{\text{CG}}$
are $\mathcal{O}(\alpha_{S})$, with the extra $|t|/\bar{Q}^{2}$ suppression,
they are at best $\bar{Q}^{-3}$. Therefore, up to $\bar{Q}^{-3}$ corrections,
the OFCA is

$\begin{split}{T}_{kk}&=\frac{1}{2\bar{P}\cdot\bar{q}}\Big{(}h\cdot\bar{q}\mathcal{H}_{1}+e\cdot\bar{q}\mathcal{E}_{1}\Big{)}.\end{split}$
(48)

This is a drastic improvement on Ref. Hannaford-Gunn  _et al._ (2022b), where
we truncated all terms that were not leading-order ($\bar{Q}^{-1}$ and
higher), and only isolated a linear combination of $\mathcal{H}_{1,2,3}$ and
$\mathcal{E}_{1,2}$. Here, we have either eliminated completely unwanted
tensor structures, or suppressed them by a further $|t|/\bar{Q}^{2}$, with
only a simple linear combination of $\mathcal{H}_{1}$ and $\mathcal{E}_{1}$
remaining.

##  Appendix B Determination of the Compton amplitude

As in Eq. (13), we determine the combination of perturbed correlators defined
in Eq. (14):

$\displaystyle
R_{\lambda}=\frac{\mathcal{G}_{(\lambda,\lambda)}+\mathcal{G}_{(-\lambda,-\lambda)}-\mathcal{G}_{(\lambda,-\lambda)}-\mathcal{G}_{(-\lambda,\lambda)}}{\mathcal{G}_{(0,0)}}.$

This combination isolates the $\lambda^{2}$ contribution up to
$\mathcal{O}(\lambda^{4})$ corrections.

First, we fit $R_{\lambda}$ as a function of the Euclidean time $\tau$.
Apropos Eq. (13), these correlators should have a linear $\tau$ dependence, so
we fit the function $f(\tau)=A\tau+B$, where the slope is proportional to the
Compton amplitude. This fit is performed using a weighted averaging method
similar to Ref. Beane  _et al._ (2021), where multiple Euclidean time windows
are fit and then averaged over with each window weighted by

$\tilde{w}^{i}=\frac{p(\delta A^{i})^{-2}}{\sum_{i^{\prime}}p(\delta
A^{i^{\prime}})^{-2}},$ (49)

where $A^{i}$ is the slope parameter from the $i^{\text{th}}$ fit window,
$\delta A^{i}$ is the statistical error and $p$ is the $p$-value determined by

$p=\tilde{\Gamma}(N_{\text{dof}}/2,\chi^{2}/2)/\tilde{\Gamma}(N_{\text{dof}}/2),$

where $\tilde{\Gamma}$ the regularised upper incomplete gamma function.

Figure 9: Left: the fits in $\tau$ to the quantity $R_{\lambda}$ defined in
Eq. (14). The light shaded band is the full extent of all fit windows, while
the darker shaded band is the fit window with the highest weight given by Eq.
(49). Right: the fits in $\lambda$ to $R_{\lambda}$ after the fits in
Euclidean time. The top panels are for the unpolarised spin-parity projector,
while the bottom are for the $y$-polarised projector. All results for
$t=-0.57\;\text{GeV}^{2}$ and $uu$ quark combination.

The minimum Euclidean time to include in the fits is chosen by eye, while the
largest is chosen using the unperturbed correlator’s noise-to-signal (we
choose a number of standard deviations, $N_{\sigma}$, and the largest time
slice is chosen as last time slice for which the unperturbed correlator is
more than $N_{\sigma}$ standard deviations from zero). Finally, to keep the
degrees of freedom greater than zero, we need a fit window of minimum length
$3a$.

After the Euclidean fits, we perform fits in the Feynman-Hellmann parameter
$\lambda$. The combination of correlators $R_{\lambda}$ must be a polynomial
in $\lambda^{2}$, where the $\lambda^{2}$ coefficient is proportional to the
Compton amplitude. Since we only determine two $\lambda$ values in the range
$[0.0625,0.025]$ for each of our inversions, we only perform a one parameter
fit: $f(\lambda)=b\lambda^{2}$. However, as tested in Ref. Hannaford-Gunn  _et
al._ (2022b) the higher order in $\lambda$ terms do not have a major impact on
our results, especially at our current level of precision.

The results of these fits for $t=-0.57\;\text{GeV}^{2}$ and $uu$ quarks are
presented in Fig. 9. We observe that in Euclidean time the results, while
still reasonably noisy, are well described by the linear fit function.
Similarly, in the Feynman-Hellmann parameter, the results are well-described
by the quadratic fit.

##  Appendix C Sink momenta

In our Compton amplitude determination, the sink momentum determines the value
of $\bar{\omega}$, which encodes the $x$-dependence of the GPD. In particular,
accurate and precise determinations of large $\bar{\omega}$ values are crucial
to determining higher moments and allowing a full reconstruction of the GPD.

It is convenient to define the dimensionless sink momentum:

$\mathbf{n}^{\prime}=\left(\frac{L}{2\pi}\right)\mathbf{p}^{\prime},$ (50)

where for our work the momentum interval is $2\pi/L\approx 380\;\text{MeV}$.

For the forward results we do not calculate the Compton amplitude for
$\mathbf{n}^{\prime 2}>9$, while for the off-forward Compton amplitude we do
not calculate beyond $\mathbf{n}^{\prime 2}>10$. Without methods to improve
the isolation of the ground state, beyond these sink momenta the signal is of
poor quality. As per our discussion in Section II, we do not calculate results
for $|\bar{\omega}|>1$.

Similarly, in all our fits we do not include results for which
$\mathbf{n}^{\prime 2}\geq 7$, which corresponds to $|\mathbf{p}^{\prime}|\geq
1\;\text{GeV}$, as ground state isolation is generally poorer and
$\mathcal{O}(ap_{\mu})$ artefacts are expected to be significant for these
sink momenta. This is discussed further in the text.

Table 6:  Dimensionless sink momenta, $\mathbf{n}^{\prime}$, and corresponding $\bar{\omega}$ values for the four different sets of Compton amplitude results. Momenta in italics have $|\mathbf{p}^{\prime}|\geq 1\;\text{GeV}$ and are hence excluded from our fits. | $t=0$
---
$\mathbf{n}^{\prime}$ | $\bar{\omega}$ | $\mathbf{n}^{\prime 2}$
$(0,0,0)$ | $0.0$ | $0$
$(-1,2,0)$ | $0.06$ | $5$
$(1,-1,0)$ | $0.12$ | $2$
$(0,1,0)$ | $0.18$ | $1$
$\mathit{(2,-2,0)}$ | $\mathit{0.24}$ | $\mathit{8}$
$(1,0,0)$ | $0.29$ | $1$
$(0,2,0)$ | $0.35$ | $4$
$(2,-1,0)$ | $0.41$ | $5$
$(1,1,0)$ | $0.47$ | $2$
$\mathit{(0,3,0)}$ | $\mathit{0.53}$ | $\mathit{9}$
$(2,0,0)$ | $0.59$ | $4$
$(1,2,0)$ | $0.65$ | $5$
$(2,1,0)$ | $0.76$ | $5$
$\mathit{(3,0,0)}$ | $\mathit{0.88}$ | $\mathit{9}$
$\mathit{(2,2,0)}$ | $\mathit{0.94}$ | $\mathit{8}$
$t=-0.29\;\text{GeV}^{2}$
---
$\mathbf{n}^{\prime}$ | $\bar{\omega}$ | $\mathbf{n}^{\prime 2}$
$(1,0,-1)$ | $0.03$ | $2$
$(0,-1,2)$ | $0.15$ | $5$
$(1,0,0)$ | $0.21$ | $1$
$\mathit{(2,1,-2)}$ | $\mathit{0.27}$ | $\mathit{9}$
$\mathit{(0,-1,3)}$ | $\mathit{0.33}$ | $\mathit{10}$
$(1,0,1)$ | $0.39$ | $2$
$(2,1,-1)$ | $0.45$ | $6$
$(1,0,2)$ | $0.57$ | $5$
$(2,1,0)$ | $0.63$ | $5$
$\mathit{(1,0,3)}$ | $\mathit{0.75}$ | $\mathit{10}$
$(2,1,1)$ | $0.81$ | $6$
$t=-0.57\;\text{GeV}^{2}$
---
$\mathbf{n}^{\prime}$ | $\bar{\omega}$ | $\mathbf{n}^{\prime 2}$
$(0,0,1)$ | $0.0$ | $1$
$(-1,2,1)$ | $0.06$ | $6$
$(1,-1,1)$ | $0.12$ | $3$
$(0,1,1)$ | $0.18$ | $2$
$\mathit{(2,-2,1)}$ | $\mathit{0.24}$ | $\mathit{9}$
$(1,0,1)$ | $0.29$ | $2$
$(0,2,1)$ | $0.35$ | $5$
$(2,-1,1)$ | $0.41$ | $6$
$(1,1,1)$ | $0.47$ | $3$
$\mathit{(0,3,1)}$ | $\mathit{0.53}$ | $\mathit{10}$
$(2,0,1)$ | $0.59$ | $5$
$(1,2,1)$ | $0.65$ | $6$
$(2,1,1)$ | $0.76$ | $6$
$\mathit{(3,0,1)}$ | $\mathit{0.88}$ | $\mathit{10}$
$\mathit{(2,2,1)}$ | $\mathit{0.94}$ | $\mathit{9}$
$t=-1.14\;\text{GeV}^{2}$
---
$\mathbf{n}^{\prime}$ | $\bar{\omega}$ | $\mathbf{n}^{\prime 2}$
$(1,-1,0)$ | $0.0$ | $2$
$(2,0,-1)$ | $0.12$ | $5$
$\mathit{(0,-2,2)}$ | $\mathit{0.12}$ | $\mathit{8}$
$(1,-1,1)$ | $0.24$ | $3$
$(2,0,0)$ | $0.35$ | $4$
$(1,-1,2)$ | $0.47$ | $6$
$(2,0,1)$ | $0.59$ | $5$
$\mathit{(3,1,0)}$ | $\mathit{0.71}$ | $\mathit{10}$
$\mathit{(2,0,2)}$ | $\mathit{0.82}$ | $\mathit{8}$

In Tab. 6 we present the dimensionless sink momentum and corresponding
$\bar{\omega}$ values for all results in this work. Note that we average over
equivalent kinematics: for the unpolarised projector, the value of
$\mathcal{R}_{\mu\mu}$ (see Eq. (15)) does not change with
$\bar{\omega}\to-\bar{\omega}$ or $\Delta\to-\Delta$, and hence we average
over these. For the polarised projector, there is a relative minus sign for
each of the changes $\bar{\omega}\to-\bar{\omega}$ or $\Delta\to-\Delta$ in
$\mathcal{R}_{\mu\mu}$, and hence we average over these accounting for the
minus sign.

As mentioned in the text, there is no $\bar{\omega}=0$ point for the
$t=-0.29\;\text{GeV}^{2}$ results. Recall from Eq. (19) that
$\bar{\omega}\propto\bar{\mathbf{p}}\cdot\bar{\mathbf{q}}$, where

$\displaystyle\bar{\mathbf{p}}=\frac{1}{2}({\mathbf{p}}^{\prime}+{\mathbf{p}})=\frac{1}{2}(2{\mathbf{p}}^{\prime}+{\mathbf{\Delta}}).$
(51)

Therefore, to access $\bar{\omega}=0$ we need to set
$\bar{\mathbf{p}}=\mathbf{0}$, which implies
$\mathbf{p}^{\prime}=-{\mathbf{\Delta}}/2$ from the above equation. However,
for the $t=-0.29\;\text{GeV}^{2}$ results,
$\mathbf{\Delta}=\frac{2\pi}{L}(1,-1,0)$, and hence $\bar{\omega}=0$ implies
that $\mathbf{p}^{\prime}=\frac{2\pi}{L}(-0.5,0.5,0)$, which is not accessible
with our discretised momentum.

##  Appendix D Additional $d$ quark results

In Fig. 11 we present all the fits performed in $\bar{\omega}$ space for the
$d$ quarks for the $\mathcal{H}_{1}$ structure function. We note that the $d$
quark results typically have a poorer signal-to-noise ratio than those for the
$uu$ quark combination.

In Fig. 10, we present the posterior distributions for the model fit of the
$d$ quark results. We note that the $\alpha^{\prime}$ posteriors have a less
well-defined Gaussian form when compared to the $u$ quark results, Fig. 5.
Further, the $\beta$ parameter is largely unconstrained, but in contrast to
the $u$ quark, the $d$ quark results for this parameter are not as strongly
skewed to the lower bound.

In contrast to the $u$ quark results (see Fig. 3), the $d$ quark model fits
appear to have less dependence on the prior distribution—i.e. the ‘wide’,
‘medium’ and ‘thin’ results show better agreement. This can also be seen in
Fig. 12, where we plot the results for the moments of the $d$ quarks.

Similar to the $u$ quark results, we note reasonable agreement between the
model fits, the direct moment fits and the leading-twist GFFs determine from
the local twist-two operators (labelled ‘$A_{n,0}$ 3-pt’). Moreover, the
$\mathcal{A}_{4,0}$ moment is not well-determined compared to
$\mathcal{A}_{2,0}$.

Figure 10: All results and labels as in Fig. 5 except for $d$ quarks. Figure
11: All results and labels as in Fig. 3 except for $d$ quarks. Figure 12: All
results and labels as in Fig. 6 except for $d$ quarks. Figure 13: Posterior
distributions for the moments of the off-forward structure functions
$\mathcal{H}_{1}$ (top) and $\mathcal{E}_{1}$ (bottom) using the fit function
in Eq. (27). Here $N_{\text{max}}$ is the number of moments included in the
fits, with $n=2N_{\text{max}}$ as the highest moment. The range of the uniform
priors is given in Eq. (32). All results for $t=-0.57\;\text{GeV}^{2}$ and
$uu$ quark combination only.

##  Appendix E Mellin Moment Fit

For the direct fits to the Mellin moments, using the fit function Eq. (27), we
use uniform prior distributions given by the constraints in Eq. (32).

The posteriors of these fits are presented in Fig. 13 for
$t=-0.57\;\text{GeV}^{2}$ and $u$ quarks. We observe that the order of
truncation, $N_{\text{max}}$ has negligible effect on the $n=2,4$ moments for
$N_{\text{max}}>2$. However, for the $\mathcal{A}_{n,0}$ moments, we note that
the $n=4,6$ the distributions are highly skewed towards the upper bound. This
suggests that the bounds in Eq. (32) are over constraining for our results.
This could be the result of (1) systematics such as the lack of large
$\bar{\omega}$ results or discretisation artefacts, or (2) that the GPD
positivity bound Eq. (30) is broken at $\bar{Q}^{2}\approx 5\;\text{GeV}^{2}$.
At our current precision and control of systematics, we can not draw strong
conclusions. However, this demonstrates a possible application for our data in
testing theoretical GPD constraints at moderately large $\bar{Q}^{2}$.

##  Appendix F Generalised form factors from local operators

Figure 14: Euclidean time dependence of the ratio in Eq. (56) for the operator
$\bar{\psi}\gamma_{4}\psi$ (left) and the operator $\mathcal{O}_{v_{2}b}$
(right). Both ratios have the unpolarised spin-parity projector and sink
momenta $\mathbf{p}^{\prime}=(1,0,0)$, which implies $\mathbf{\Delta}^{2}=0$.

At various points in this work, we have compared our Compton structure
function moments, $\mathcal{A}_{n,0}$ and $\mathcal{B}_{n,0}$, to the leading-
twist generalised form factors $A_{n,0}$ and $B_{n,0}$ (see Figs. 4, 6 and
12). We determine these generalised form factors using the standard technique
of calculating the matrix elements of twist-two local operators using three-
point functions.

For $A_{1,0}$ and $B_{1,0}$ (note we do not compare $B_{1,0}$ to our Compton
amplitude moments), we use all components of the local vector current:

$\displaystyle j_{\mu}=\bar{\psi}\gamma_{\mu}\psi,\quad\text{for
}\mu=1,2,3,4.$ (52)

For $A_{2,0}$ and $B_{2,0}$ we use the operators

$\mathcal{O}_{v_{2}a}=\frac{1}{2}\sum_{i=1}^{3}\left(\mathcal{O}_{4i}^{V}+\mathcal{O}_{i4}^{V}\right)\quad\text{and}\quad\mathcal{O}_{v_{2}b}=\mathcal{O}_{44}^{V}-\frac{1}{3}\sum_{i=1}^{3}\mathcal{O}_{ii}^{V},$
(53)

where

$\mathcal{O}_{\mu\nu}^{V}=\overline{\psi}\gamma_{\mu}i\overleftrightarrow{D}_{\nu}\psi.$
(54)

To determine the matrix elements of these operators, we compute the three-
point correlation function:

$\mathcal{G}^{\Gamma}_{3-\text{pt}}(\mathbf{p}^{\prime},\tau;\mathbf{\Delta},\tau_{\text{ins}})=\sum_{\mathbf{x}_{1},\mathbf{x}_{2}}e^{-i\mathbf{p}^{\prime}\cdot\mathbf{x}_{2}}e^{i\mathbf{\Delta}\cdot\mathbf{x}_{1}}\Gamma^{\alpha\beta}\langle\chi_{\beta}(\mathbf{x}_{2},\tau)\mathcal{O}(\mathbf{x}_{1},\tau_{\text{ins}}){\chi}^{\dagger}_{\alpha}(\mathbf{0},0)\rangle,$
(55)

The transfer 3-momentum is defined as
$\mathbf{\Delta}=\mathbf{p}^{\prime}-\mathbf{p}$. The local operator,
$\mathcal{O}$, is inserted at time slice $\tau_{\text{ins}}$, where
$\tau>\tau_{\text{ins}}>0$. We use the spin-parity projectors given in Eq.
(16), making use of all polarisation directions ($x,y$ and $z$).

To isolate ground state contribution we construct the ratio

${R}_{\text{3-pt}}=\frac{\mathcal{G}_{3-\text{pt}}^{\Gamma}(\mathbf{p}^{\prime},\tau;\mathbf{\Delta},\tau_{\text{ins}})}{\mathcal{G}_{2-\text{pt}}(\mathbf{p}^{\prime},\tau)}\left[\frac{\mathcal{G}_{2-\text{pt}}(\mathbf{p}^{\prime},\tau)\mathcal{G}_{2-\text{pt}}(\mathbf{p}^{\prime},\tau_{\text{ins}})\mathcal{G}_{2-\text{pt}}(\mathbf{p},\tau-\tau_{\text{ins}})}{\mathcal{G}_{2-\text{pt}}(\mathbf{p},\tau)\mathcal{G}_{2-\text{pt}}(\mathbf{p},\tau_{\text{ins}})\mathcal{G}_{2-\text{pt}}(\mathbf{p}^{\prime},\tau-\tau_{\text{ins}})}\right]^{\frac{1}{2}}.$
(56)

In order to control excited state contamination, we make use of a two-state
ansatz for the correlation functions similar to that in Ref. Ottnad (2021) for
our Euclidean time fits. See Fig. 14 for selected Euclidean time fits to
$R_{\text{3-pt}}$.

These calculations are performed on the same set of gauge configurations as
our OFCA calculation—see Tab. 1—using $N_{\text{meas}}=1074$, making their
statistics comparable to our Compton amplitude results (see Tab. 2). We
calculate the three-point correlators for 16 values of the soft momentum
transfer $t=-\mathbf{\Delta}^{2}$. The results are multiplicatively
renormalised in $\text{RI}^{\prime}\text{/MOM}$.

## References

  * Burkardt (2000) M. Burkardt, Phys. Rev. D 62, 071503 (2000), [Erratum: Phys.Rev.D 66, 119903 (2002)], arXiv:hep-ph/0005108 .
  * Ji (1997) X.-D. Ji, Phys. Rev. Lett. 78, 610 (1997), arXiv:hep-ph/9603249 .
  * Polyakov (2003) M. V. Polyakov, Phys. Lett. B 555, 57 (2003), arXiv:hep-ph/0210165 .
  * Burkert  _et al._ (2018) V. D. Burkert, L. Elouadrhiri, and F. X. Girod, Nature 557, 396 (2018).
  * Kumerički  _et al._ (2016) K. Kumerički, S. Liuti, and H. Moutarde, Eur. Phys. J. A 52, 157 (2016), arXiv:1602.02763 [hep-ph] .
  * Guidal  _et al._ (2013) M. Guidal, H. Moutarde, and M. Vanderhaeghen, Rept. Prog. Phys. 76, 066202 (2013), arXiv:1303.6600 [hep-ph] .
  * Bertone  _et al._ (2021) V. Bertone, H. Dutrieux, C. Mezrag, H. Moutarde, and P. Sznajder, Phys. Rev. D 103, 114019 (2021), arXiv:2104.03836 [hep-ph] .
  * Dutrieux  _et al._ (2022) H. Dutrieux, O. Grocholski, H. Moutarde, and P. Sznajder, Eur. Phys. J. C 82, 252 (2022), [Erratum: Eur.Phys.J.C 82, 389 (2022)], arXiv:2112.10528 [hep-ph] .
  * Guo  _et al._ (2022) Y. Guo, X. Ji, and K. Shiells, JHEP 09, 215 (2022), arXiv:2207.05768 [hep-ph] .
  * Guo  _et al._ (2023) Y. Guo, X. Ji, M. G. Santiago, K. Shiells, and J. Yang, JHEP 05, 150 (2023), arXiv:2302.07279 [hep-ph] .
  * Riberdy  _et al._ (2024) M. J. Riberdy, H. Dutrieux, C. Mezrag, and P. Sznajder, Eur. Phys. J. C 84, 201 (2024), arXiv:2306.01647 [hep-ph] .
  * Mamo and Zahed (2024) K. A. Mamo and I. Zahed, (2024), arXiv:2404.13245 [hep-ph] .
  * Hägler  _et al._ (2003) P. Hägler, J. W. Negele, D. B. Renner, W. Schroers, T. Lippert, and K. Schilling, Phys. Rev. D 68, 034505 (2003), arXiv:hep-lat/0304018 .
  * Göckeler  _et al._ (2004) M. Göckeler, R. Horsley, D. Pleiter, P. E. L. Rakow, A. Schäfer, G. Schierholz, and W. Schroers, Phys. Rev. Lett. 92, 042002 (2004), arXiv:hep-ph/0304249 .
  * Göckeler  _et al._ (2005) M. Göckeler, P. Hägler, R. Horsley, D. Pleiter, P. E. L. Rakow, A. Schäfer, G. Schierholz, and J. M. Zanotti, Phys. Lett. B 627, 113 (2005), arXiv:hep-lat/0507001 .
  * Göckeler  _et al._ (2007) M. Göckeler, P. Hägler, R. Horsley, Y. Nakamura, D. Pleiter, P. E. L. Rakow, A. Schäfer, G. Schierholz, H. Stüben, and J. M. Zanotti, Phys. Rev. Lett. 98, 222001 (2007), arXiv:hep-lat/0612032 .
  * Hägler  _et al._ (2008) P. Hägler, W. Schroers, J. Bratt, J. W. Negele, A. V. Pochinsky, R. G. Edwards, D. G. Richards, M. Engelhardt, G. T. Fleming, B. Musch, _et al._ , Phys. Rev. D 77, 094502 (2008), arXiv:0705.4295 [hep-lat] .
  * Brömmel  _et al._ (2008) D. Brömmel, M. Diehl, M. Göckeler, P. Hägler, R. Horsley, Y. Nakamura, D. Pleiter, P. E. L. Rakow, A. Schäfer, G. Schierholz, H. Stüben, and J. M. Zanotti, Phys. Rev. Lett. 101, 122001 (2008), arXiv:0708.2249 [hep-lat] .
  * Bratt  _et al._ (2010) J. D. Bratt, R. G. Edwards, M. Engelhardt, P. Hägler, H. W. Lin, M. F. Lin, H. B. Meyer, B. Musch, J. W. Negele, K. Orginos, _et al._ , Phys. Rev. D 82, 094502 (2010), arXiv:1001.3620 [hep-lat] .
  * Alexandrou  _et al._ (2011) C. Alexandrou, J. Carbonell, M. Constantinou, P. A. Harraud, P. Guichon, K. Jansen, C. Kallidonis, T. Korzec, and M. Papinutto, Phys. Rev. D 83, 114513 (2011), arXiv:1104.1600 [hep-lat] .
  * Ji (2013) X. Ji, Phys. Rev. Lett. 110, 262002 (2013), arXiv:1305.1539 [hep-ph] .
  * Radyushkin (2017) A. V. Radyushkin, Phys. Rev. D 96, 034025 (2017), arXiv:1705.01488 [hep-ph] .
  * Chen  _et al._ (2020) J.-W. Chen, H.-W. Lin, and J.-H. Zhang, Nucl. Phys. B 952, 114940 (2020), arXiv:1904.12376 [hep-lat] .
  * Lin (2021) H.-W. Lin, Phys. Rev. Lett. 127, 182001 (2021), arXiv:2008.12474 [hep-ph] .
  * Alexandrou  _et al._ (2020) C. Alexandrou, K. Cichy, M. Constantinou, K. Hadjiyiannakou, K. Jansen, A. Scapellato, and F. Steffens, Phys. Rev. Lett. 125, 262001 (2020), arXiv:2008.10573 [hep-lat] .
  * Bhattacharya  _et al._ (2023) S. Bhattacharya, K. Cichy, M. Constantinou, X. Gao, A. Metz, J. Miller, S. Mukherjee, P. Petreczky, F. Steffens, and Y. Zhao, Phys. Rev. D 108, 014507 (2023), arXiv:2305.11117 [hep-lat] .
  * Lin (2023) H.-W. Lin, Phys. Lett. B 846, 138181 (2023), arXiv:2310.10579 [hep-lat] .
  * Chambers  _et al._ (2017) A. J. Chambers, R. Horsley, Y. Nakamura, H. Perlt, P. E. L. Rakow, G. Schierholz, A. Schiller, K. Somfleth, R. D. Young, and J. M. Zanotti, Phys. Rev. Lett. 118, 242001 (2017), arXiv:1703.01153 [hep-lat] .
  * Can  _et al._ (2020) K. Can, A. Hannaford-Gunn, R. Horsley, Y. Nakamura, H. Perlt, P. Rakow, G. Schierholz, K. Somfleth, H. Stüben, R. Young, and J. M. Zanotti, Phys. Rev. D 102, 114505 (2020), arXiv:2007.01523 [hep-lat] .
  * Hannaford-Gunn  _et al._ (2020) A. Hannaford-Gunn, R. Horsley, Y. Nakamura, H. Perlt, P. E. L. Rakow, G. Schierholz, K. Somfleth, H. Stüben, R. D. Young, and J. M. Zanotti, PoS LATTICE2019, 278 (2020), arXiv:2001.05090 [hep-lat] .
  * Batelaan  _et al._ (2023) M. Batelaan, K. Can, A. Hannaford-Gunn, R. Horsley, Y. Nakamura, H. Perlt, P. Rakow, G. Schierholz, H. Stüben, R. Young, and J. Zanotti, Phys. Rev. D 107, 054503 (2023), arXiv:2209.04141 [hep-lat] .
  * Hannaford-Gunn  _et al._ (2022a) A. Hannaford-Gunn, E. Sankey, K. U. Can, R. Horsley, H. Perlt, P. E. L. Rakow, G. Schierholz, K. Somfleth, H. Stüben, R. D. Young, and J. M. Zanotti, PoS LATTICE2021, 028 (2022a), arXiv:2207.03040 [hep-lat] .
  * Camacho Muñoz  _et al._ (2006) C. Camacho Muñoz _et al._ (Jefferson Lab Hall A), Phys. Rev. Lett. 97, 262002 (2006), arXiv:nucl-ex/0607029 .
  * Defurne  _et al._ (2015) M. Defurne _et al._ (Jefferson Lab Hall A), Phys. Rev. C 92, 055202 (2015), arXiv:1504.05453 [nucl-ex] .
  * Defurne  _et al._ (2017) M. Defurne _et al._ (Jefferson Lab Hall A), Nature Commun. 8, 1408 (2017), arXiv:1703.09442 [hep-ex] .
  * Georges (2018) F. Georges,  _Deeply virtual Compton scattering at Jefferson Lab_ , Ph.D. thesis, Orsay, IPN (2018).
  * Braun  _et al._ (2012) V. M. Braun, A. N. Manashov, and B. Pirnay, Phys. Rev. Lett. 109, 242001 (2012), arXiv:1209.2559 [hep-ph] .
  * Guo  _et al._ (2021) Y. Guo, X. Ji, and K. Shiells, JHEP 12, 103 (2021), arXiv:2109.10373 [hep-ph] .
  * Kumerički (2019) K. Kumerički, Nature 570, E1 (2019).
  * Hannaford-Gunn  _et al._ (2022b) A. Hannaford-Gunn, K. U. Can, R. Horsley, Y. Nakamura, H. Perlt, P. E. L. Rakow, H. Stüben, G. Schierholz, R. D. Young, and J. M. Zanotti, Phys. Rev. D 105, 014502 (2022b), arXiv:2110.11532 [hep-lat] .
  * Perrottet (1973) M. Perrottet, Lett. Nuovo Cim. 7S2, 915 (1973).
  * Tarrach (1975) R. Tarrach, Nuovo Cim. A28, 409 (1975).
  * Belitsky  _et al._ (2002) A. V. Belitsky, D. Müller, and A. Kirchner, Nucl. Phys. B 629, 323 (2002), arXiv:hep-ph/0112108 .
  * Belitsky  _et al._ (2014) A. V. Belitsky, D. Müller, and Y. Ji, Nucl. Phys. B 878, 214 (2014), arXiv:1212.6674 [hep-ph] .
  * Braun  _et al._ (2014) V. M. Braun, A. N. Manashov, D. Müller, and B. M. Pirnay, Phys. Rev. D 89, 074022 (2014), arXiv:1401.7621 [hep-ph] .
  * Horsley  _et al._ (2020) R. Horsley, Y. Nakamura, H. Perlt, P. E. L. Rakow, G. Schierholz, K. Somfleth, R. D. Young, and J. M. Zanotti, PoS LATTICE2019, 137 (2020), arXiv:2001.05366 [hep-lat] .
  * Karpie  _et al._ (2019) J. Karpie, K. Orginos, A. Rothkopf, and S. Zafeiropoulos, JHEP 04, 057 (2019), arXiv:1901.05408 [hep-lat] .
  * Joó  _et al._ (2019) B. Joó, J. Karpie, K. Orginos, A. V. Radyushkin, D. G. Richards, R. S. Sufian, and S. Zafeiropoulos, Phys. Rev. D 100, 114512 (2019), arXiv:1909.08517 [hep-lat] .
  * Goeke  _et al._ (2001) K. Goeke, M. V. Polyakov, and M. Vanderhaeghen, Prog. Part. Nucl. Phys. 47, 401 (2001), arXiv:hep-ph/0106012 .
  * Belitsky  _et al._ (2004) A. V. Belitsky, X.-d. Ji, and F. Yuan, Phys. Rev. D 69, 074014 (2004), arXiv:hep-ph/0307383 .
  * Diehl  _et al._ (2005) M. Diehl, T. Feldmann, R. Jakob, and P. Kroll, Eur. Phys. J. C 39, 1 (2005), arXiv:hep-ph/0408173 .
  * Guidal  _et al._ (2005) M. Guidal, M. V. Polyakov, A. V. Radyushkin, and M. Vanderhaeghen, Phys. Rev. D 72, 054013 (2005), arXiv:hep-ph/0410251 .
  * Jenkovszky (2006) L. L. Jenkovszky, Phys. Rev. D 74, 114026 (2006), arXiv:hep-ph/0607340 .
  * Schoeffel (2007) L. Schoeffel, Phys. Lett. B 658, 33 (2007), arXiv:0706.3488 [hep-ph] .
  * Kumerički  _et al._ (2008) K. Kumerički, D. Müller, and K. Passek-Kumerički, Nucl. Phys. B 794, 244 (2008), arXiv:hep-ph/0703179 .
  * Diehl and Kugler (2008) M. Diehl and W. Kugler, Phys. Lett. B 660, 202 (2008), arXiv:0711.2184 [hep-ph] .
  * Gonzalez-Hernandez  _et al._ (2013) J. O. Gonzalez-Hernandez, S. Liuti, G. R. Goldstein, and K. Kathuria, Phys. Rev. C 88, 065206 (2013), arXiv:1206.1876 [hep-ph] .
  * Kumerički and Müller (2010) K. Kumerički and D. Müller, Nucl. Phys. B 841, 1 (2010), arXiv:0904.0458 [hep-ph] .
  * Kroll (2015) P. Kroll, EPJ Web Conf. 85, 01005 (2015), arXiv:1410.4450 [hep-ph] .
  * Can  _et al._ (2023) K. U. Can, M. Batelaan, A. Hannaford-Gunn, R. Horsley, Y. Nakamura, H. Perlt, P. E. L. Rakow, G. Schierholz, H. Stüben, R. D. Young, and J. M. Zanotti, in  _30th International Workshop on Deep-Inelastic Scattering and Related Subjects_ (2023) arXiv:2307.07904 [hep-lat] .
  * Bietenholz  _et al._ (2011) W. Bietenholz, V. Bornyakov, M. Göckeler, R. Horsley, W. G. Lockhart, Y. Nakamura, H. Perlt, D. Pleiter, P. E. L. Rakow, G. Schierholz, A. Schiller, T. Streuer, H. Stüben, F. Winter, and J. M. Zanotti, Phys. Rev. D 84, 054509 (2011), arXiv:1102.5300 [hep-lat] .
  * Beane  _et al._ (2021) S. Beane, W. Detmold, R. Horsley, M. Illa, M. Jafry, D. Murphy, Y. Nakamura, H. Perlt, P. Rakow, G. Schierholz, P. Shanahan, H. Stüben, M. Wagman, F. Winter, R. Young, and J. Zanotti, Phys. Rev. D 103, 054504 (2021), arXiv:2003.12130 [hep-lat] .
  * Salvatier  _et al._ (2016) J. Salvatier, T. V. Wiecki, and C. Fonnesbeck, PeerJ Computer Science 2:e55 (2016).
  * Hoffman and Gelman (2014) M. D. Hoffman and A. Gelman, Journal of Machine Learning Research 15, 1593 (2014).
  * Pobylitsa (2002) P. V. Pobylitsa, Phys. Rev. D 65, 114015 (2002), arXiv:hep-ph/0201030 .
  * Di Vecchia and Drago (1969) P. Di Vecchia and F. Drago, Lett. Nuovo Cim. 1S1, 917 (1969).
  * Jia  _et al._ (2017) D. Jia, C.-Q. Pang, and A. Hosaka, Int. J. Mod. Phys. A 32, 1750153 (2017), arXiv:1706.02788 [hep-ph] .
  * Egerer  _et al._ (2021) C. Egerer, R. G. Edwards, K. Orginos, and D. G. Richards, Phys. Rev. D 103, 034502 (2021), arXiv:2009.10691 [hep-lat] .
  * Bali  _et al._ (2016) G. S. Bali, B. Lang, B. U. Musch, and A. Schäfer, Phys. Rev. D 93, 094515 (2016), arXiv:1602.05525 [hep-lat] .
  * Wu  _et al._ (2018) J. J. Wu, W. Kamleh, D. B. Leinweber, R. D. Young, and J. M. Zanotti, J. Phys. G 45, 125102 (2018), arXiv:1807.09429 [hep-lat] .
  * Göckeler  _et al._ (2006) M. Göckeler, R. Horsley, H. Perlt, P. E. L. Rakow, G. Schierholz, and A. Schiller, PoS LAT2006, 119 (2006), arXiv:hep-lat/0610064 .
  * Haar  _et al._ (2018) T. R. Haar, Y. Nakamura, and H. Stüben, EPJ Web Conf. 175, 14011 (2018), arXiv:1711.03836 [hep-lat] .
  * Edwards and Joó| (2005) R. G. Edwards and B. Joó|, Nucl.Phys.Proc.Suppl. 140, 832 (2005), arXiv:hep-lat/0409003 [hep-lat] .
  * Karsten and Smith (1981) L. H. Karsten and J. Smith, Nuclear Physics B 183, 103 (1981).
  * Guerin (1987) F. Guerin, Nuclear Physics B 282, 495 (1987).
  * Ottnad (2021) K. Ottnad, Eur. Phys. J. A 57, 50 (2021), arXiv:2011.12471 [hep-lat] .