Generalised parton distributions from the off-forward Compton amplitude in lattice QCD

We determine the properties of generalised parton distributions (GPDs) from a lattice QCD calculation of the off-forward Compton amplitude (OFCA). By extending the Feynman-Hellmann relation to second-order matrix elements at off-forward kinematics, this amplitude can be calculated from lattice propagators computed in the presence of a background field. Using an operator product expansion, we show that the deeply-virtual part of the OFCA can be parameterised in terms of the low-order Mellin moments of the GPDs. We apply this formalism to a numerical investigation for zero-skewness kinematics at two values of the soft momentum transfer, $t = -1.1, -2.2 \;\text{GeV}^2$, and a pion mass of $m_{\pi}\approx 470\;\text{MeV}$. The form factors of the lowest two moments of the nucleon GPDs are determined, including the first lattice QCD determination of the $n=4$ moments. Hence we demonstrate the viability of this method to calculate the OFCA from first principles, and thereby provide novel constraint on the $x$- and $t$-dependence of GPDs.


I. INTRODUCTION
Since the 1990s, generalised parton distributions (GPDs) have been recognised as crucial observables in understanding hadron structure [1][2][3].They encode the spatial distribution of quarks and gluons in a fastmoving hadron [4].Moreover, their Mellin moments contain information about the spin and orbital angular momentum of hadron constituents [2], which would resolve the decades-old 'proton spin puzzle' [5,6].Finally, more recent research has explored the relationship between GPDs and 'mechanical' properties: pressure, energy and force distributions within hadrons [7,8].
GPDs can be measured from off-forward Compton scattering processes, such as deeply virtual Compton scattering (DVCS), which have been carried out at HERA [9][10][11][12][13], COMPASS [14], JLab [15][16][17][18], and are planned to be carried out in the future at the electronion collider [19].However, due to the high dimensionality of GPDs, they are difficult to extract directly from experiment, and global fits require assumptions about their functional form [20,21]. Therefore, a stronger theoretical understanding of GPD behaviour would allow for more precise experimental determinations.
More recently, there have also been major efforts to reconstruct the full x-dependence of parton distributions in lattice QCD, using the pseudo- [39] and quasidistribution [40] methods-see Refs.[41,42] for reviews.This includes recent calculations of quasi-GPDs [43][44][45].These methods aim to extract the light-cone distributions directly, whereas in the present work we are interested in the Compton scattering amplitude, from which GPDs may be accessed experimentally.
In this paper, we determine properties of GPDs from a calculation of the off-forward Compton amplitude (OFCA) in lattice QCD.The OFCA is defined as (1) and describes the process of γ * (q)N (P ) → γ * (q )N (P ), with q µ = 0 = q µ .Here, j µ is the hadronic vector The Feynman diagram for off-forward γ * (q)N (P ) → γ * (q )N (P ) scattering.
current, and we limit ourselves to the case where the scattered hadron is a nucleon.
Besides GPDs, this amplitude gives access to a range of interesting physical quantities, including generalised polarisabilities [46][47][48][49] and the subtraction function [50,51].In the high energy region (|q 2 | and/or |q 2 | Λ 2 QCD ), it is dominated by contributions from GPDs.By calculating the Compton amplitude, we overcome the issues of power-divergent renormalisation that the leading-twist matrix elements suffer from [52,53].Moreover, with the correction for lattice systematics, our calculation contains the same higher-twist contributions as the physical amplitude, which are of interest beyond their connection to leading-twist GPDs [54].Therefore, the present calculation bears many similarities to the hadronic tensor approach, which aims to access the forward structure functions from the direct calculation of four-point functions [55,56].
The method presented here to calculate the OFCA is an extension of Feynman-Hellmann methods used previously to determine the forward Compton amplitude [57,58].Two-point correlators calculated in the presence of a weakly-coupled background field field can be expanded in powers of the coupling, with their secondorder contribution in terms of four-point functions.As such, Feynman-Hellmann methods are a feasible alternative to the direct calculation of four-point functions.
The numerical results presented here are the first lattice QCD determination of the off-forward Compton amplitude.This calculation is performed at the SU(3) flavour symmetric point and larger-than-physical pion mass [59], for two values of the soft momentum transfer, t = −1.1,−2.2 GeV 2 , with zero-skewness kinematics.In this preliminary work, we assume leading-twist dominance, since our hard scale is in the perturbative region: Q2 ≈ 6 − 7 GeV 2 .As such, we fit Mellin moments of the OFCA, and interpret these as the moments of GPDs.
The structure of this paper follows: in section II we review key properties of the OFCA; in section III we derive the Feynman-Hellmann relation that allows us to determine the OFCA; in section IV we use an operator product expansion to parameterise the scalar amplitudes of the OFCA in terms of GPD moments; in section V we outline the details of our numerical calculation; and finally in section VI we present our results.
From these, we can form at most four linearly independent scalar variables: two scaling variables, and the soft and hard momentum transfers, respectively, In terms of these scalars, the usual skewness variable [60] is ξ = ϑ/ω, and hence ϑ = 0, ω = 0 implies that ξ = 0.In terms of the conventional deeply virtual Compton scattering (DVCS) kinematics, where q 2 = 0, we have that ϑ 1 and ω ξ −1 for large −q 2 .
As a consequence of the Ward identities of the OFCA, q µ T µν = 0 = q ν T µν , contributions to the Compton amplitude that are proportional to q ν or q µ are not linearly independent.Hence we can write the OFCA as where the gauge projector is and Tµν is the OFCA with no q ν or q µ terms.We will choose a basis for the tensor decomposition of Tµν , since all other terms are entirely determined by the Ward identities.In our chosen basis, the OFCA (before gauge projection) is where we have introduced the Dirac bilinears In Eq. ( 4), there are nine K, five unpolarised (H and E) and four polarised ( H and Ẽ) amplitudes, which gives 18 in total.
The basis in Eq. ( 4) is chosen to match onto the highenergy limit, which we will derive in section IV.While this does introduce kinematic singularities into our basis, these are not relevant to the leading-twist contribution or our numerical calculation.
The amplitudes of Eq. ( 4) also reduce in the forward (t → 0) limit to the more well-known functions of the forward Compton amplitude: where F 1,2 are the Compton structure functions [58] and Img 1,2 = 2πg 1,2 , for g 1,2 the spin-dependent, deep inelastic structure functions [65].On the other hand, the K amplitudes vanish in the forward limit.

Dispersion Relation
As in the forward case, we can use the analytic features of the amplitudes in Eq. ( 4) to write out a dispersion relation.For instance, following Refs.[47,66], H 1 and E 1 satisfy subtracted dispersion relations: where we have introduced and similarly for The subtraction function in Eq. ( 6) is a generalisation of the forward Compton amplitude subtraction function [67]: , which has been studied elsewhere [50,51].The amplitudes H 2,3 and E 2 require no subtraction in their dispersion relations [47,66].
The forward limit of H 1 is where F 1 is the deep inelastic scattering structure function [47].However, unlike the forward case, there is no optical theorem to relate ImH 1,2 to an inclusive cross section.Instead, these amplitudes can be measured directly by exclusive processes such as DVCS.

Generalised Parton Distributions
At high energies ( Q2 Λ 2 QCD ), the amplitudes of Eq. ( 4) are dominated by convolutions of GPDs [2,68]: where G is a GPD.Or, in the Euclidean region, |ω| < 1, Formally, GPDs are defined by the off-forward matrix element of a light-cone operator.For a light-like vector n µ such that n • P = 1 (and hence ξ = −n • ∆/2) and taking light-cone gauge n • U = 0, we have [2,69] dλ 2π e iλx P | ψq (−λn/2)/ nψ q (λn/2)|P = H q (x, ϑ/ω, t)ū(P )γ µ n µ u(P ) where H q and E q are the unpolarised twist-two GPDs for a quark of flavour q.It is not possible to directly calculate the quantity in Eq. ( 7) on the lattice, due to the Euclidean signature of spacetime.Instead, we can relate GPDs to a basis of leadingtwist local operators.These local operators are where See appendix A for the symmetrisation convention of the Lorentz indices.
By Taylor expanding Eq. ( 7), one can relate the GFFs from Eq. ( 9) to the GPDs H and E: recalling that ξ = ϑ/ω in the scalars defined at the start of this section.These equations are the famous 'polynomiality' of GPDs [60].
A proof-of-principle determination of GPD moments is the ultimate aim of this paper.Specifically, we will calculate the linear combination of zero-skewness moments, Equivalent expressions for polarised GPDs, H and Ẽ, are given in appendix A.

III. FEYNMAN-HELLMANN RELATION
In this section, we will show how to calculate the offforward Compton amplitude from Feynman-Hellmann methods in lattice QCD.Feynman-Hellmann methods are a subset of background field methods, in which a two-point function is calculated in the presence of a weakly-coupled field or current.This induces perturbations to the two-point function, thereby giving access to observables that may be difficult to calculate with a direct n-point function.
For the Feynman-Hellmann derivation presented here, we expand the perturbed propagator by means of a Dyson expansion [70].This is used to approximate a derivative of the propagator, similar to Refs.[71][72][73][74], and hence extract the OFCA for off-forward kinematics.
This differs from our previous proof of the forward Feynman-Hellmann relation [58], where we expressed the perturbed correlators as G λ (τ ) A λ e −E λ τ , and related the derivatives of the perturbed energy, E λ , to the Compton amplitude.While it is still possible to derive a Feynman-Hellmann relation for the OFCA in terms of derivatives of perturbed energies [75], such a proof is made difficult by the fact that degeneracies in the unperturbed spectrum cause there to be two low-lying perturbed energies.Similar considerations are needed for nucleon electromagnetic form factors from Feynman-Hellmann [76].By contrast, the Dyson expansion and correlator derivative formalism presented below circumvents this difficulty.
We introduce two spatially-oscillating background fields to the QCD lagrangian density: where j 3 (x) = Z V ψq (x)iγ 3 ψ q (x), and Z V is the lattice renormalisation constant for a local vector current.Therefore, the perturbed Hamiltonian is where Simulating with the perturbed Lagrangian in Eq. ( 11) leads to a modified lattice two-point propagator: where λ = (λ 1 , λ 2 ), and Γ is the spin-parity projector.
+ O(λ 3 ).Inserting two complete sets of states and taking χ(τ ) = e −H FH τ χ(0), Eq. ( 13) becomes Note that states and energies without a λ subscript are unperturbed.We can expand the time evolution operator, e −H FH τ , with a Dyson series: and hence Eq. ( 14) becomes Note that we have dropped the spin structure for brevity, but will reintroduce it in the final result.
From Eq. ( 15), we see that the O(λ 2 ) terms of the two-point propagator contain four-point functions (see Figure 2).In particular, the λ 2 i term has both currents inserting momentum ±q i , and hence provides access to the forward Compton amplitude.Only the mixed, second-order term, proportional to λ 1 λ 2 , will have different incoming/outgoing momenta, and therefore offforward kinematics.
To isolate the mixed second-order term, we define the combination of nucleon propagators, Having established how to isolate the second-order, offforward contribution to the perturbed propagator, we are now interested in how to ensure ground state saturation at the source and sink.
As detailed in Ref. [58], provided that none of the intermediate states are lower energy than E N , we have Using this result, Eq. 16 becomes neglecting O(λ 4 ) corrections and where Unlike a direct four-point function approach, ground state saturation at the source is ensured by a judicious choice of kinematics, not by large Euclidean time separations-see appendix C for a full calculation.To summarise, these kinematic restrictions require that the current insertion momenta, q 1 and q 2 , and sink momentum, p , are chosen such that vents the intermediate states from going on-shell, • and |p | = |p +q 1 −q 2 |, which keeps the incoming and outgoing states energy degenerate.
After these restrictions are imposed, Eq. ( 17) can be written, up to O(λ 4 ) corrections, as where T 33 is the µ = ν = 3 component of the OFCA for a single quark flavour with unit charge, and p = p + q 1 − q 2 .The term C is constant in both λ and τ ; it is made up of contributions for which the source is not the ground state (see Eq. (C4)).Therefore, by fitting R λ (τ, p) in τ and λ, we can isolate the OFCA.

IV. THE OFF-FORWARD COMPTON AMPLITUDE
Given the method to calculate the OFCA presented in the previous section, we now show how to parameterise the invariant amplitudes of the tensor decomposition, Eq. ( 4), in terms of GPDs.
The suitable tool for a perturbative expansion of the OFCA in the Euclidean region is the operator product expansion (OPE), which is an expansion about points in coordinate space and momentum space (z µ = 0 and ω = 0, respectively) that are accessible in a spacetime with Euclidean signature [77].There exist in the literature several OPEs of the OFCA [78][79][80][81].However, as these largely focus on the spin-zero case and/or significantly pre-date GPDs, in this section we give our own OPE.

Operator Product Expansion
Formally, the leading-twist contribution to the coordinate-space current product is given by the 'handbag' contributions [77,82]: where S µρνκ = g µρ g νκ + g µκ g νρ − g µν g ρκ , and the operators are defined in Eqs. ( 8) and (A3).
To obtain the leading-twist OFCA, Eq. ( 1), we must take the off-forward matrix element of Eq. ( 19) and Fourier transform it.Details of this calculation are presented in appendix B.
The final result is for the symmetric in µ ↔ ν contribution to the OFCA defined in Eq. ( 1), while for the anti-symmetric contribution, where we have used the bilinear definitions given in Eq. ( 5).Recall that the usual skewness variable is ξ = ϑ/ω in our chosen scalars.One can verify, by taking the Sudakov decomposition and DVCS kinematics ω ξ −1 , ϑ 1, that Eqs. ( 20) and ( 21) recover the standard twist-two DVCS amplitude [69].
We can now use the OPE results, Eqs. ( 20) and ( 21), to interpret the high energy limit of each of the scalar amplitudes in the tensor decomposition, Eq. ( 4): • The scalar amplitudes either vanish at leadingtwist, or can be parameterised in terms of convolutions of GPDs.For instance: See Eq. (B2) for a full list.
• We have off-forward equivalents of the Callan-Gross relation [86]: In the forward case, Feynman-Hellmann methods have recently been used to determine powersuppressed Callan-Gross breaking terms [87].
• The moments of polarised scalar amplitudes have the following relation at leading-twist and similarly for the replacement H → Ẽ.In the forward limit, this reduces to a relation between the spin-dependent structure functions [88].
• The K scalar amplitudes vanish at leading-twist, but do contribute at twist-three in terms of transverse GPDs [64,89,90].

Parameterisation of the lattice calculation
From the Feynman-Hellmann relation, Eq. ( 18), we calculate s,s tr Γu(P , s )T 33 ū(P, s) s tr[Γu(P , s)ū(P , s)] ≡ R(ω, t, Q2 ).(22) Therefore, to get a parameterisation that can be compared to the lattice, we use the tensor decomposition of section II and the OPE with the following additional conditions: • We choose the µ = ν = 3 component of our Compton amplitude.
The zero-skewness condition removes the tensor structures with scalar amplitudes K 3,4,6,8,9 and Ẽ2 .Further, by calculating the µ = ν = 3 component and taking a spin trace, the tensor structures associated with the polarised amplitudes H and Ẽ are made irrelevant.
Finally, since we take Q2 ∼ 7 GeV 2 , we will consider the remaining amplitudes, K 1,2,5,7 to be suppressed, since they have no leading-twist contribution.
Although a more complete study of the Q2 -dependence is essential, for this exploratory work we will neglect the Q2 suppressed K 1,2,5,7 amplitudes, keeping only the H 1,2,3 and E 1,2 amplitudes.Therefore, Eq. ( 22) is with P µν as defined in Eq. ( 3), and using Euclidean conventions now to match the lattice.Next, we subtract off the ω = 0 contribution: which is equivalent to replacing H 1 → H 1 and E 1 → E 1 in Eq. (23).
As in our previous study of the forward Compton amplitude, we find anomalous asymptotic behaviour of the S 1 subtraction function.A method for controlling this behaviour has been presented in the forward case, where the anomalous behaviour of S 1 is found to have minimal effect on the ω-dependence [95].An extension to the OFCA is a goal of future work.
We then take only the leading-twist contributions to the amplitudes, a full list of which is given in Eq. (B2).
Imposing the off-forward Callan-Gross relation reduces the number of linearly independent amplitudes in Eq. ( 23) from five to two.The final form is then where E N = m 2 N + p 2 is the sink energy, and For a first approximation of extracting the GPD moments, we will calculate R(ω, t, Q2 )/K 33 ( P3 , q3 , P • q, Q2 ).
Since our lattice calculations are in frames that are roughly near the rest frame (i.e.E N ≈ m N ), we can approximately treat the combination of GFFs in Eq. ( 25) as a Lorentz scalar: A determination of the A and B GFFs independently, rather than the linear combination defined in Eq. ( 27), is desirable.To this end, note that we can also use the spin-parity projector, which would give linearly independent combinations of the A and B form factors compared with Eq. ( 25), in a manner analogous to the separation of F 1 and F 2 electromagnetic form factors. Hence a separation of the A and B form factors by varying the spin-parity projector is a goal of future work.

V. SIMULATION DETAILS
For this calculation, we use the same gauge ensembles as Ref. [58].Note, in particular, that we are at the SU(3) flavour symmetric point, κ l = κ s , with a largerthan-physical pion mass, m π = 466(13) MeV, and a lattice spacing of a = 0.074(2) fm.See Table I for a summary of the gauge configurations.

Feynman-Hellmann Implementation
The Feynman-Hellmann implementation is almost identical to our previous study of the forward Compton amplitude [58].In practice, the objects we calculate are perturbed quark propagators, given by where M is the usual fermion matrix.
For our case, where we choose to calculate the µ = ν = 3 component of the OFCA, the operators are Then, the usual formulae for hadrons in terms of quark propagators apply, except with one or more of these propagators replaced with a perturbed propagator.
The Feynman-Hellmann perturbation is applied to the connected contributions only.While it is possible to perturb the disconnected contributions, this would be much more computationally expensive [96,97].
The determination of the ratio in Eq. ( 16) requires four separate sets of correlators at each magnitude of λ.We calculate two magnitudes of λ = 0.0125, 0.025, chosen based on λ-tuning tests carried out in the forward case [58,98]

Kinematics
We calculate two sets of correlators on the same gauge configurations (Table II).
Table II.Current insertion momenta, q1,2, and derived kinematics for two sets of correlators.To fit GPD moments, we need multiple ω values.However, we are restricted by the conditions of the Feynman-Hellmann relation to a frame for which our sink momentum, p , and our momenta from the current insertions, q 1 and q 2 , must obey: which limits the number of ω values that are accessible for each q 1,2 pair.
For each set of correlators, the ω value is determined by the value of the sink momentum, p : The explicit values of ω for our kinematics are shown in Table III.Moreover, since our amplitude is invariant under the exchanges ∆ µ → −∆ µ , ω → −ω, we average over ±p , ±(p − q 1 + q 2 ) to increase our statistics.
Table III.ω values for the two sets of correlators.Note that |ω| > 1 values are omitted.

VI. RESULTS AND DISCUSSION
To demonstrate what can be accomplished with the method outlined in the preceding sections, we determine the first two even moments of the nucleon GPD.
First, we fit the combination of correlators, R λ (τ, p ) from Eq. ( 16) to the function f (τ ) = c 1 τ +c 2 , where τ is Euclidean time.From the Feynman-Hellmann relation, Eq. ( 18), the slope, c 1 , is proportional to the OFCA, while c 2 is a superfluous parameter.In fitting this linear function, we apply a consistent fit window in Euclidean  time for all sink momenta.The χ 2 /dof for these fits are reported in Table IV and it is found that χ 2 /dof ∼ 1 for all the momenta, which demonstrates that the data is largely well described by a linear fit.An example of the Euclidean time fits for set #1 is given in Figure 3.After the fits in Euclidean time have been performed, we next investigate the behaviour of the ratio, R λ (p ), as a function of the Feynman-Hellmann coupling, λ.From the Feynman-Hellmann relation, Eq. ( 18), the λ 2 contribution to this ratio is proportional to the OFCA, and the next-to-leading contribution is O(λ 4 ), which is suppressed for our calculations at λ ∼ 10 −2 .
Therefore, to test the effects of the λ 4 contributions, we compare the quadratic coefficient of the ratio as extracted with a purely quadratic fit function, f (λ) = bλ 2 , to that extracted with the function g(λ) = bλ 2 + cλ 4 .We find that the quartic coefficient, c, is consistent with . Plot of R, as defined in Eq. ( 25), divided by the kinematic factor, K33, from Eq. ( 26).The red curve is a parameterisation of results from Ref. [58].The blue and green curves are from the moment fits.
zero, and that the quadratic coefficients, b, calculated using the two fit functions agree within errors.However, since the quartic fit determines two parameters from two λ values, it is not a reliable estimate of the higher order contaminations.Therefore, to further examine the effect of these contaminations, we calculate the quotient (λ 2 1 R λ2 )/(λ 2 2 R λ1 ), which is 1 for perfectly quadratic results.In Table IV, we can see that, although the central value of this quotient is close to 1 for all momenta, not all are within errors of 1.This indicates a 2−4% contamination from higher order terms, which is negligible compared to our overall errors.
Hence for this preliminary study, we find it sufficient to use the purely quadratic fit function, f (λ) = bλ 2 .In Figure 4, we plot the normalised ratio, R λ /λ 2 , as a function of λ, and compare this to the quadratic coefficient from the fit.We observe that the data is reasonably well described by a purely quadratic fit.
Table IV.Parameters demonstrating the quality of fits in Euclidean time and the Feynman-Hellmann parameter, λ, for the up quark results from set # 1.
Using the Feynman-Hellmann relation, Eq. ( 18), we can now interpret the quadratic coefficient as the offforward Compton amplitude.Then, by varying the sink momentum, we can calculate the amplitude at multiple values of the scaling variable, ω.The results for the up quark in the nucleon are shown in Figure 5.
The forward t = 0 curve in this plot is a fit to the Q 2 = 7.13 GeV 2 results from Ref. [58].As that study also used the Feynman-Hellmann method and the same gauge configurations as the present calculation, we can compare it to our off-forward, t = 0, results to determine the t-dependence of the OFCA.

Moment Fitting
Using the results of our OPE in section IV, we can interpret the moments of the OFCA as GPD moments, defined in Eq. ( 27).Hence, using Eq. ( 25), a fit in ω to the function yields the first J even GPD moments at fixed t and Q2 values.At leading-twist, these moments are Unlike the forward case, there is no optical theorem connecting the OFCA to the scattering cross section, and therefore no requirement for the scalar amplitudes to be positive definite.However, our moments, defined in Eq. ( 27), are dominated by A n,0 (t), the moments of the zero-skewness GPD H(x, t), which is typically treated as positive in model-dependent parameterisations (for instance, Refs.[99][100][101][102]), while the E(x, t) GPD is suppressed by t/8m 2 N in our moments.Therefore, it is reasonable for this proof-of concept calculation to treat the underlying distribution, H(x, t) + (t/8m 2 N )E(x, t), as strictly positive on the main x ∈ [−1, 1], and thus its moments as monotonically decreasing for fixed t: Future work will aim at a more extensive treatment of the conditions on moments, such as incorporating model-independent positivity constraints on GPDs [103][104][105] and on the Compton amplitude [106].
To fit these moments, we use a Markov chain Monte Carlo method [107,108].In contrast to a least squares fit, this method allows us to efficiently sample prior distributions that reflect physically-motivated constraints [58].The upper two plots use monotonic priors distributions, while the lower two plots use uniform positive distributions.J is the number of moments fit in the parameterisation, Eq. ( 29).
In Figure 6, we compare the up quark, t = −1.1 GeV 2 moments fit using monotonically decreasing priors, as in Eq. ( 30), to those fit with uniform positive priors, M n (t) ∈ [0, 100].Since we truncate the series of moments at a finite order, Figure 6 also compares the values of the first two moments fit at different orders of truncation, J, as in Eq. ( 29).We observe that, for both the monotonic moments, the value of J has little effect on the leading moment, M 2 .Moreover, the values of M 2 , as extracted with the monotonic and uniform moments, are highly consistent.
On the other hand, the value of M 4 differs significantly depending on whether uniform or monotonic priors are used.For the uniformly sampled moments, the M 4 distributions are heavily skewed towards zero, and do not converge with J.By contrast, the monotonically sampled moments, M 4 , do not depend greatly on the order of truncation for J > 2, and the distributions appear only slightly skewed towards zero.The higher moments require larger values of ω and more precise data to constrain them.Therefore, the inconsistencies in the M 4 results likely reflect the fact that we have a limited number of larger ω values, which have significant errors.Moreover, these inconsistencies may reflect that the monotonicity condition is too severe for small moments.Investigating these issues is a goal of future studies.
For this preliminary study, we choose to fit the first four even moments, n = 2, 4, 6, 8, using monotonic con- Figure 7.The t-dependence of the first two even moments, M q n (t), defined in Eq. ( 27), for up and down quarks.The t = 0 points are from a fit to results in Ref. [58].
ditions, and report the first two even moments.For consistency, we only fit the first four moments of the forward results as well.
We present results for the t-dependence of the leading moments in Figure 7.The values of the n = 2 GPD moments are statistically consistent with moments from three-point calculations at a comparable pion mass [27].However, the n = 4 moments have never been determined from three-point methods, and therefore the results presented here are a first look at the t behaviour of such moments.

Comment on systematics
As the present numerical results are exploratory, a detailed assessment of systematic uncertainties remains an objective of future work.A list of the most salient systematics and proposals to control them is given below.
1. To better isolate the leading-twist contribution, a range of Q2 values must be calculated, and the constant, leading-twist moments fit from this, as in Ref. [58].
2. The two data sets (#1 and #2) have different q3 , which means that the O(a) Ward identity violating terms, induced by discretisation, will differ between the two data sets.Hence it is preferable to use the conserved vector current, for which exact Ward identities are known [109,110].
3. The OPE performed in section IV is a continuum relation, and therefore a continuum extrapolation, similar to that in Refs.[111][112][113], is desirable.
4. Finally, there are all the usual lattice systematics: non-physical quark masses, finite volume, and excited state contamination, which must be accounted for.

VII. SUMMARY AND CONCLUSIONS
This study has presented a novel means to determine the off-forward Compton amplitude (OFCA) using lattice QCD, and thereby calculate the properties of generalised parton distributions (GPDs).We derived a Feynman-Hellmann relation to calculate the OFCA.In our parameterisation of the OFCA, we presented new results and collected old ones, which lay the groundwork for comprehensive calculations of GPDs from the OFCA.Finally, the nucleon moments presented here are the first determination of n = 4 GPD moments.
We are now in a position to realise the full potential of this method.A more detailed investigation of the systematics is a priority, including calculations at different lattice spacings and with the conserved vector current.Currently, such tests are being conducted for the forward Compton amplitude [95].Similarly, future work will be aimed at calculating a greater spread of Q2 and t, which will provide physical insights and allow us to more accurately determine the leading-twist contribution.Furthermore, we aim to separate out the H and E scalar amplitudes-equivalently the A and B generalised form factors.
Taking these steps would provide us with a wealth of physical information.For instance, we could investigate the non-perturbative features of the OFCA, including the off-forward subtraction function and generalised polarisabilites.Moreover, we could investigate GPD properties, such as their scaling behaviour, and higher-twist contributions to the Compton amplitude.Finally, this method allows us to constrain GPDs, by calculating their moments, fitting models, and other methods to extract parton distributions from the Euclidean Compton amplitude directly [114].
software library [116]) 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 ( For symmetrisation and anti-symmetrisation of a rank-2 tensor, we use the notation The general expression for a fully symmetrised rank-n tensor used in this paper is where S n is the group of permutations of the numbers 1, 2, ..., n, and σ is an element of S n .Here, we denote the i th component of some group element, σ ∈ S n , as σ(i).
Next, we use the identity to evaluate the integrals over the Fourier conjugates.After applying these steps, we arrive at Eqs. ( 20) and (21).
Next, after inserting a complete set of states, Eq.C1 becomes dτ 2 e −(E X (p X )−E N (p ))τ1 e (E X (p X )−E Y (p))τ2 Because of our choice of perturbing potential, the only values the source momentum can take are p = p + nq 1 + mq 2 for m, n ∈ Z at order O(λ m+n ).As we stated before, we choose our kinematics so that |p| ≤ |p + nq 1 + mq 2 |.Therefore, for any state in the nucleon spectrum X and any momentum q = p + nq 1 + mq 2 , we must have E X (q) ≥ E N (p ).
This ensures two things: (1) that the exponentials in Eqs.C3 and C4 are decaying, and (2) that if E Y (p) = E N (p ), then Y = N , and hence we have ground state saturation of the source.

(C5)
The exponentially decaying terms will be heavily suppressed for τ a compared to the purely linear in τ terms and the constant.Therefore, we will neglect these.For the moment we neglect the term that is constant in τ ; however, we will consider this in our fit to the lattice data.
Noting that the OFCA for a single quark flavour and unit charge can be expressed as T 33 (p , q, q ) = N (p )| Ô(p , q)|N (p) + N (p )| Ô(p , −q )|N (p) , equation C9 becomes where we have used the fact that A λ N (p ) = 1 + O(λ) at most, but once again odd powers of λ vanish as a consequence of our combination of propagators, Eq. 16.

Figure 2 .
Figure 2. Illustration of the expansion of the perturbed propagator, Eq. (15), where we have suppressed the tower of states at source and sink.Euclidean time increases left to right.

Figure 4 .
Figure 4. Plot of λ-dependence of the combination of correlators, R /λ 2 , as defined in Eq. (16), after fitting in Euclidean time.The shaded bands are the quadratic coefficient for each momentum.Momenta correspond to those in Figure 3.

Figure 6 .
Figure 6.Density distributions for the first two up quark moments at t = −1.10GeV 2 .The upper two plots use monotonic priors distributions, while the lower two plots use uniform positive distributions.J is the number of moments fit in the parameterisation, Eq. (29).

Table I .
Details of the gauge ensembles used in this work.
supported by the Australian Commonwealth Government) and the Phoenix HPC service (University of Adelaide).AHG is supported by an Australian Government Research Training Program (RTP) Scholarship.RH is supported by STFC through grant ST/P000630/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 grant DP190100297.