format: text

text:
# Deep HST Imaging of Sextans A. III. The Star Formation History1footnote 11footnote 1Based on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with proposal ID 7496.

## 1 Introduction

Sextans A (DDO 75) is a dwarf irregular galaxy near the edge of the Local Group. Its distance has recently been determined via a variety of methods to be $`d=1.32\pm 0.04`$ Mpc (Dolphin et al. 2003; hereafter paper II), which corresponds to a distance modulus of $`\mu _0=25.61\pm 0.07`$ and places it in close proximity to Sextans B, NGC 3109, and the Antlia dwarf.

It has long been known that Sextans A is an active star-forming galaxy, with blue supergiants being observed in the early photometric work of Sandage & Carlson (1982) and HII regions observed by Aparicio & Rodríguez-Ulloa (1992) and Hodge, Kennicutt, & Strobel (1994). Deeper CCD photometry by Hoessel, Schommer, & Danielson (1983) and Aparicio et al. (1987) indicated that the current star formation rate was significantly higher than it had been in the past.

The first attempt to quantitatively determine the star formation history of Sextans A was done by Dohm-Palmer et al. (1997), based on moderately deep WFPC2 observations. They computed star formation histories based on the main sequence (MS) and blue helium-burning (BHeB) sequence luminosity functions. They also used BHeB stars as age indicators (brighter means younger) to determine the spatially-resolved star formation history. Similar work was done by Van Dyk, Puche, & Wong (1988), who used wider-field but shallower data and found similar rates of recent star formation. All of this work was limited to the past Gyr, leaving the question of the remainder of Sextans A’s history open.

This is the third paper in a series based on new, deeper WFPC2 observations of Sextans A. Paper I (Dohm-Palmer et al., 2002) measured the recent star formation history (both global and spatially resolved) using the techniques of Dohm-Palmer et al. (1997); the geometry of the two sets of WFPC2 data allow for such a study to be carried out over most of the galaxy. In paper II, we studied the short-period variable star content and determined a new distance based on a variety of distance indicators.

In this paper, we seek to uncover the full star formation history of Sextans A through CMD modeling techniques (Dolphin 2002 and references therein). Although such work would have been possible with the earlier WFPC2 data, the present data reach two magnitudes deeper, significantly below the red clump region. In section 2, we describe the data and reductions. Section 3 contains a “manual” determination of the star formation history by looking at the various components separately. Finally, in section 4, we quantitatively measure the full star formation history using our CMD modeling techniques.

## 2 Data and Reduction

### 2.1 Observations

A detailed description of our observations is presented in paper I. We summarize here. As part of cycle 7 program GO-7496, we obtained 24 pairs of 1200s images with WFPC2. Eight pairs of images were in F555W and sixteen pairs were in F814W. The total integration times were thus 19200 seconds in F555W and 38400 seconds in F814W. One quarter of the images were taken at each of four pointings, each offset by 0.25 arcsec, or approximately 5.5 PC pixels. A 25th pair of images was taken in F656N, but will not be used for this paper.

### 2.2 Photometry

We retrieved the data from the STScI archive using on the fly calibration, which applies the best set of reference frames. Reduction and photometry was made with the HSTphot package (Dolphin, 2000a).

To achieve the greatest depth and photometry accuracy, we cosmic ray cleaned all images in each filter and then co-added all images at the same pointing. The result was eight deep images, corresponding to four pointings in each of two filters. The procedure was similar to that used in paper II, except that we coadded images whenever possible here. In paper II, we were looking for variable stars and thus kept images from different orbits separate even if the pointings were the same.

HSTphot was used to simultaneously photometer stars on all eight deep images, producing combined $`V`$ and $`I`$ magnitudes in addition to magnitudes calculated from the separate images. Calibration to the WFPC2 flight system was made using HSTphot-produced CTE corrections and calibrations (Dolphin, 2000b). For this work, we consider only the combined magnitudes. We measured a total of 42751 objects, of which 33340 passed our photometry cuts (observed in both filters, $`\chi 2.0`$, $`|`$sharpness$`|0.3`$, and S/N$`5`$). Figure 1 shows the CMD.

We made extensive artificial star tests over the color range $`0.5<VI<3.0`$ and the magnitude range $`18<V<30`$. Stars were distributed on the images and on the CMD as functions of the observed star distribution, though care was made to ensure that all parts of the CMD adequately sampled. In total, we made 1.36 million artificial star tests. We show the results of the artificial star tests in Figure 2.

## 3 Stellar Populations

Because the CMD of Sextans A contains stars of all ages, from recently-formed stars on the MS to ancient stars on the red giant branch (RGB), one can use the CMD to determine the star formation history of the galaxy. In nearby systems, one can simply make an unconstrained maximum likelihood solution for all unknowns: distance, extinction, SFR(t), and Z(t) using the techniques described by Dolphin (2002 and references therein). However, Dolphin (2002) has demonstrated that one needs photometry extending well below the horizontal branch (i.e., at least to $`M_V=+2.0`$) in order to constrain the oldest stellar populations. Unfortunately, this a luxury we do not have for galaxies at the distance of Sextans A. Thus we need to first examine the components of the CMD in turn before we can determine a meaningful star formation history. We will do so in order of importance on the CMD.

### 3.1 Young Stars (age $`<`$ 1 Gyr)

We first turn our attention to the young stars in Sextans A, which are seen in the CMD on the MS and BHeB sequence. Paper I presented a detailed derivation of the recent star formation history using these stars, so we will discuss them only briefly here.

Theoretical stellar evolution models predict that the position BHeB sequence is a strong function of metallicity. Using the Padova models (Girardi et al., 2000), we have created synthetic CMDs of young stars for metallicities Z=0.0004 ($`[\text{M/H}]=1.7`$) and Z=0.004 ($`[\text{M/H}]=0.7`$); these are shown and compared with the observed CMD in Figure 3. In preparing these figures, we adopted the distance ($`\mu _0=25.61`$) and extinction ($`A_V=0.12`$) determined in paper II and a power law IMF with a slope of $`\mathrm{\Gamma }=1.30`$ (compared with a Salpeter value of $`\mathrm{\Gamma }=1.35`$). Note that we intentionally do not plot isochrones over the observed CMD. The MS and BHeB sequence both have finite widths due to evolutionary effects which must be accounted for in modeling the CMD. Thus, simple line representations of these populations can be misleading.

Comparing the observed and synthetic CMDs, it is clear that the BHeB sequence in the $`[\text{M/H}]=0.7`$ CMD falls well above the observed data, while that in the $`[\text{M/H}]=1.7`$ sequence falls along the lower edge of the observed sequence. We note that both synthetic CMDs have sharper BHeB sequences than the observed data; we interpret this as evidence of a metallicity spread. Fitting the BHeB sequence position and spread, we find that the observed CMD is best fit by a current mean metallicity of $`[\text{M/H}]=1.45`$, with a metallicity spread of $`\pm 0.20`$ dex. Note that if the errors in our photometry were significantly underestimated, then the observations would be consistent with a smaller metallicity spread.

However, we also find that the synthetic MS is too narrow to fit the observed data, regardless of the adopted metallicity (or metallicity spread). There are several possible explanations for this. One is that the isochrones are grossly in error on the main sequence, a possibility whose detailed exploration is beyond the scope of our paper. A second possibility would be the presence of a greater metallicity spread. However we note that $`VI`$ colors of upper main sequence stars change by only $`0.03`$ magnitudes per dex of metallicity; thus even a huge metallicity spread would not explain the main sequence color spread we observe.

A more plausible option is the presence of differential extinction in the youngest populations, caused by the preferential proximity of young stars to dusty star-forming regions. Such an effect has been observed in the LMC by Harris, Zaritsky, & Thompson (1997), who find extinction values in excess of $`A_V=1.5`$ for OB stars but very few red giants with extinction values in excess of $`A_V=0.6`$, based on their $`UBVI`$ photometric survey. As Sextans A is also actively forming stars, one should expect to see such an effect here as well, albeit with lower overall extinction values due to its lower metallicity. This conclusion is strengthened by the photometry of Aparicio et al. (1987), who found an excess number of red points in their $`(UB)`$, $`(BV)`$ diagram that they interpreted as evidence of internal reddening.

Determination of the exact functional form of the extinction distribution is a more difficult task. We note that the Cepheid period-magnitude relations from paper II showed no sign of differential extinction, and thus assume no differential extinction at ages $`500`$ Myr. We also note that the width of the red helium-burning sequence fainter than $`V=21.25`$ ($`M_V=4.5`$) is explained solely by the small metallicity spread noted above; this implies little differential extinction at ages $`70`$ Myr. The MS shape is best fit with some differential extinction for ages younger than $`90`$ Myr and significantly higher amounts at ages younger than $`40`$ Myr.

To model the differential extinction as simply as possible (given the constraints above), we apply a flat distribution of differential extinction values to all stars younger than 100 Myr, with the maximum extinction increasing linearly from zero at an age of 100 Myr to 0.5 magnitudes in $`A_V`$ at an age of 40 Myr.

Incorporating all of these features (current metallicity, metallicity spread, and differential extinction) into the CMD produces the “young” synthetic CMD shown in Figure 4.

### 3.2 Old Stars (age $`>`$ 2 Gyr)

The dominant CMD feature that is not caused by young stars is the RGB; we now turn our attention to that. Although there are some younger stars falling in that part of the CMD, the RGB is predominantly composed of stars older than $`2`$ Gyr. The position of the RGB is a function of age and metallicity, in that a star that is older or more metal-rich will be redder. We show the range of possible metallicity values for 2, 5, 10, and 15 Gyr stars in Figure 5, calculated from the Girardi et al. (2000) theoretical isochrones.

The blue edge of the RGB is poorly-defined, as there are significant numbers of population I and AGB stars in that part of the CMD. However, the red edge is fairly sharp, and thus provides strong constraints for the metallicity of RGB stars. The red edge can be fit with 2 Gyr old stars of metallicity $`[\text{M/H}]0.9`$, 15 Gyr old stars of metallicity $`[\text{M/H}]1.3`$, or stars with intermediate ages and metallicities. In the previous section we determined a metallicity of $`1.45\pm 0.20`$ for the young stars. Thus, we would expect very few old stars to have metallicity higher than this. From Figure 5, we can conclude that the stars at the red edge of the RGB are likely old stars (ages $``$ 10 Gyr), with metallicity $`[\text{M/H}]1.3`$.

We note that the presence of old stars with metallicities similar to young stars implies that Sextans A reached its current metallicity very quickly. Such a feature is not unique to Sextans A; other metal-poor galaxies with extended star formation such as Leo I (Dolphin, 2002) and UGC 4438 (Dolphin et al., 2001) have shown little metal enrichment over most of their histories. This is not necessarily inconsistent with the bulk of their star formation occurring at later ages. In the cases of both Sextans A and UGC 4438, their present day metallicities (12 $`+`$ log(O/H) $``$ 7.5, Skillman et al. 1989; 1994) are not that much higher than that seen almost universally at high redshift (Pettini et al., 1999). A reasonable picture of Sextans A is thus that it formed from pre-enriched gas, and that most of the metals produced by its stars have escaped.

To account for the young and old stars, we will adopt a mean metallicity of $`[\text{M/H}]=1.45`$ with spread of $`\pm 0.20`$ dex for all ages. The synthetic CMD (still not containing stars between $`1`$ and $`2`$ Gyr) is shown in Figure 6.

The observed CMD (Figure 1) shows a feature that extends redwards from the RGB tip to $`(VI)3`$. Modeling this using the Girardi et al. (2000) isochrones would require a metallicity of $`[\text{M/H}]=0.6`$ at an age of 3 Gyr ago, which is inconsistent with our observations of the metallicity of younger stars and would require an RGB tip color of $`VI=1.95`$ (about half a magnitude redder than what is observed). We conclude that these stars are a combination of TP-AGB and Carbon stars, evolutionary phases that are observed at this color and absolute magnitude in other metal-poor galaxies such as IC 1613 (Albert, Demers, & Kunkel, 2000).

Finally, we note that there is no clearly-delineated horizontal branch in Sextans A, a result of the large number of young stars in that part of the CMD. Consequently a measurement of the number of stars $`>10`$ Gyr old can only be done via statistical modeling of the CMD.

### 3.3 Intermediate-age Stars

The final stars to be considered are the stars with ages between 1 and 2 Gyr. Although “intermediate-age” has many meanings, we adopt the practical definition of being older than the oldest BHeB stars but younger than the stars on the RGB. The synthetic CMD with just these stars is shown in Figure 7. It is clear from the figure that these stars contribute to the red clump, upper AGB, and lower MS; however disentangling them from the younger and older populations is not trivial because they do not occupy otherwise-empty parts of the CMD. The statistical decomposition of the CMD will be the topic of the next section.

## 4 Star Formation History

Armed with the components of the CMD that are produced by stars of a variety of ages, we set out to determine what combination of these components best fits the observed data. The simplest hypothesis is that the star formation rate has been constant over the history of the galaxy. A synthetic CMD computed with such a history is shown in Figure 8, adopting the metallicity and extinction from the previous section. From the difference between the observed and synthetic diagrams, it is clear that the synthetic diagram has too many old stars (from the oversubtraction of the RGB and lower red clump) and too few young stars (from the undersubtraction of the MS and BHeB sequence). Put differently, the star formation rate over the past Gyr or two is significantly higher than that $`>510`$ Gyr ago.

Because of the failure of the constant SFR model, a more sophisticated fit for the history is necessary. The statistical basis behind this procedure and tests of its accuracy were given by Dolphin (2002); we merely summarize here. We produce a large number of “single-population” synthetic CMDs, each containing only stars with limited ranges (0.1 dex) of age and metallicity. These CMDs are generated using a combination of theoretical isochrones, IMF, binary function, distance, and extinction distribution. Observational errors (from photon noise and blending) are simulated by using the completeness, photometric bias, and photometric scatter (all functions of color and magnitude) measured in the artificial star tests.

These synthetic CMDs, as well as simulated CMDs of foreground stars and bad points, are combined linearly to calculate the expected distribution of stars on the CMD for any star formation history. With synthetic and observed CMDs in hand, we calculate the likelihood that the observed data were produced by the star formation history used to produce that particular synthetic CMD; repeating this procedure for many histories to maximize the likelihood allows us to measure the history of Sextans A.

As with the Tolstoy et al. (1998) study of Leo A, we find that the CMD is better-fit with a shorter distance than with our preferred distance of $`\mu _0=25.61`$. We attribute this to uncertainties in the theoretical models of high-mass stars, and make our solution using the longer distance. To compensate for potential numerical instabilities caused by the isochrones, we ran three separate solutions – one in the normal way described by Dolphin (2002), one at low resolution (to emphasize the numbers of stars in each part of the CMD while deemphasizing positions of features), and one excluding the red clump region from the CMD (to reduce the reliance of the solution on that dominant feature). Each of the three solutions was consistent with the others, giving us confidence in the robustness of our results.

Sextans A appears to have had a significant event $`12.5`$ Gyr ago, when its star formation greatly increased. Star formation has continued since that event, with another increase of a factor of three $`0.060.1`$ Gyr ago. It is unclear what triggered the sudden increases in star formation. These events are frequently linked with interactions with other galaxies, and Sextans A has a very close companion, Sextans B. Interestingly, based on the HIPASS observation of a smooth HI distribution around Sextans A, Wilcots & Hunter (2002) have argued that Sextans A has not interacted with another galaxy in the last 6 – 10 Gyr. Since the global dynamics of the HI associated with the optical galaxy show large distortions (i.e., misalignment of optical and kinematic axes, Skillman et al. 1988), and the HIPASS observations are very low angular resolution, perhaps a recent interaction with Sextans B is not ruled out, and possibly even likely.

The star formation history at older ages is less certain. We believe the star formation rate to have been much lower than at present, but nonzero. However, because the old stars are dominated by the large numbers of younger stars, it is impossible to quantify the star formation rate or give detailed properties of the old population.

Confirming what we found in the previous section, we measured a roughly constant (to within 0.1 dex) metallicity of $`[\text{M/H}]=1.4`$ over the measurable history of the galaxy. The CMD was best fit with an rms spread of 0.15 dex, a spread present at all ages. We note that the recovered metallicity is consistent with that of $`4\%`$ solar determined by Skillman, Kennicutt, & Hodge (1989) from HII region spectroscopy.

We show the star formation history in Table 1 and Figure 9, normalized to a lifetime average of 1. A Monte Carlo test was used to determine the statistical uncertainties; systematic errors were estimated by comparing our three solutions. Conversion of our relative star formation rates into absolute rates depends heavily on the assumed IMF for low-mass stars. Assuming the IMF of Kroupa (2001), we calculate a lifetime average star formation rate of $`600\pm 300M_{}Myr^1kpc^2`$ for this field, and an average rate of $`1600\pm 500M_{}Myr^1kpc^2`$ during the past 2.5 Gyr.

The best synthetic CMD is shown in Figure 10; the improvement over Figure 8 is readily apparent but we were still unable to perfectly fit the CMD. This was expected, given the uncertainties in the input physics in theoretical models. Nevertheless, we see that most parts of the CMD were fit very well; no CMD bin had an error of more than 5 $`\sigma `$. Furthermore, the discrepancies are errors in position rather than errors in number of objects in a particular phase of stellar evolution. Compared with the observed CMD, the theoretical MS is narrower, the RGB is broader, the red supergiants are bluer, and the red clump is fainter. None of these discrepancies is sufficiently large to give us reason to believe the measured star formation history is incorrect.

## 5 Recent Star Formation History Revisited

Paper I presented recent (age $`<0.7`$ Gyr) star formation history calculations based on the MS and BHeB sequence luminosity functions. The $`V`$ magnitudes along each sequence were converted into star masses, which could be transformed into star formation rates based on assumed lifetimes of those phases and an assumed initial mass function. Summarizing those results briefly, from the BHeB stars it was determined that Sextans A had a fairly constant star formation rate (within a factor of 2) over the past 0.7 Gyr, followed by a rise (factor of 2 - 3) during the period between 50 and 100 Gyr ago, and followed by another steep rise in the last 25 Myr. The star formation histories constructed from the MS stars were roughly consistent with this, only spanning a smaller time interval ($``$ 200 Myr) and showing a larger dependence on which stellar evolution models were used. In paper I we emphasized that the method of directly converting the MS luminosity function to a recent star formation history is inherently uncertain because, unlike the BHeB stars, a position in the CMD corresponds to a large range of ages and there is significant luminosity evolution on the MS itself. In paper I we suggested that the star formation history could be derived more accurately from the MS star using statistical methods, and thus we make that comparison here.

In our present best model we find a recent star formation history which is consistent with that found in paper I, with the notable exception that we do not see evidence of the steep increase at 25 Myr. This feature had appeared to be quite firm as it was found in both the MS and BHeB analysis (using the Padova models, although it was not as strong and there was less agreement when the Geneva models were used). On the other hand, the analysis of the BHeB stars is much more secure for ages greater than 25 Myr (for a given star formation rate there are more stars, and the evolution of these stars is less affected by uncertainties such as the role of stellar winds and rotation).

To understand the cause of the discrepancy, we have made similar fits to those in paper I by adopting their distance, metallicity, and CMD sections and removing binary stars and differential extinction from our CMD models. We find that we can reproduce the BHeB star formation history, indicating that the method used for that calculation in paper I is robust. However, we are unable to reproduce the MS star formation history, a result we believe to stem from the treatment of evolution on the main sequence. Since the method used in the present paper accounts for evolution, we believe the history presented here to be more accurate.

We also note that, if we model the CMD at the distance, extinction, metallicity, and binarity used in this paper but still using the restricted CMD sections (MS and BHeB stars only), we find a young star formation history consistent with that shown in Table 1, confirming the results presented here. To clarify, our star formation history does show an elevated star formation rate during the last $``$ 60 Myr relative to the almost constant star formation rate from 100 - 500 Myr ago (as found in paper I), it is only that there is no increase of another factor of 2 to 4 in the last 25 Myr. Our modeling here supports the use of the BHeB stars for reconstructing recent star formation histories, but suggests a limitation that these star formation histories be limited to parts of the CMD where the contamination of the BHeB sequence from reddened MS stars is minimal.

## 6 Summary

We have presented deep WFPC2 photometry of the dwarf irregular galaxy Sextans A, and determined its history of star formation through comparisons with synthetic CMDs. We have made these comparisons both manually and using a numerical algorithm and find consistent results.

Sextans A began forming stars at a high rate sometime between 1 and 2.5 Gyr ago and has continuously formed stars since then. The star formation rate over the past 0.06 Gyr is approximately three times the average over the past 2.5 Gyr, although we do not have the time resolution to determine whether similar 0.1 Gyr “bursts” happened more than 1 Gyr ago. We also found a small but nonzero number of stars older than 2.5 Gyr, but the present data do not permit an accurate determination of their characteristics.

We found very little evidence of significant chemical enrichment within Sextans A. Based on our CMD fitting, it appears that the galaxy reached a metallicity of $`[\text{M/H}]1.4`$ more than 10 Gyr ago, and still has that metallicity today. This metallicity is consistent with literature values from nebular spectroscopy. An rms scatter of $`0.15`$ at all ages dex helps fit the observed CMD more accurately.

Our determination of the recent star formation history (age $``$ 0.7 Gyr) agrees with that of paper I for all but the last 25 Myr. In paper I we found a steep rise in the star formation rate that is not recovered here. Since the present analysis demands consistency of both the MS and BHeB stars and takes into account MS luminosity evolution, we prefer the solution presented here. This indicates that star formation histories determined solely from BHeB stars are probably quite unreliable at magnitudes where reddened MS stars contaminate the BHeB sequence.

Support for this work was provided by NASA through grant number GO-07496 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555. EDS is grateful for partial support from NASA LTSARP grant No. NAG5-9221.


text:
# Moments of Nucleon Generalized Parton Distributions in Lattice QCD

## I Introduction

Light-cone correlation functions play a special role in the experimental exploration of the quark and gluon structure of hadrons. Asymptotic freedom allows quantitative separation of the reaction mechanism from the structure of the probed hadron at high energy, so that spin-independent scattering experiments unambiguously measure matrix elements of the light-cone operator

$$𝒪(x)=\frac{d\lambda }{4\pi }e^{i\lambda x}\overline{\psi }(\frac{\lambda }{2}n)\overline{)}n𝒫e^{ig_{\lambda /2}^{\lambda /2}𝑑\alpha nA(\alpha n)}\psi (\frac{\lambda }{2}n),$$
(1)
where $`n`$ is a light-cone vector, and $`𝒫`$ denotes a path-ordering of the gauge fields in the exponential. Since these matrix elements are singled out by their experimental accessibility, it is essential to use all our tools of analytical methods and lattice field theory to explore and understand them as fully as possible.

Diagonal nucleon matrix elements, $`q(x)=P|𝒪(x)|P`$, measure the familiar quark distribution $`q(x)`$ specifying the probability of finding a quark carrying a fraction $`x`$ of the nucleon’s momentum in the light cone frame. Although light cone correlation functions cannot be calculated directly in lattice QCD, expansion of $`𝒪(x)`$ generates the tower of twist-two operators,

$$𝒪_q^{\{\mu _1\mu _2\mathrm{}\mu _n\}}=\overline{\psi }_q\gamma ^{\{\mu _1}\text{i}\stackrel{}{D}^{\mu _2}\mathrm{}\text{i}\stackrel{}{D}^{\mu _n\}}\psi _q,$$
(2)
with $`\stackrel{}{D}=\frac{1}{2}\left(\stackrel{}{D}\stackrel{}{D}\right)`$, and curly braces $`\{\mu _1\mathrm{}\mu _n\}`$ mean symmetrization of indices and subtraction of traces. The diagonal matrix elements $`P|𝒪_q^{\{\mu _1\mu _2\mathrm{}\mu _n\}}|P`$ can be calculated on the lattice and specify the $`(n1)^{th}`$ moments $`𝑑xx^{n1}q(x)`$. Note that expressions analogous to Eqns. (1) and (2) for spin-dependent observables differ only in their gamma matrix structure, but will not be considered in the present work.

Generalized parton distributions (GPDs), as introduced in Müller et al. (1994); Ji (1997); Radyushkin (1997), correspond to non-diagonal matrix elements $`P^{}|𝒪(x)|P`$. When expressed in terms of the relevant Lorentz invariants, $`P^{}|𝒪(x)|P`$ is specified in terms of two generalized parton distributions, $`H(x,\xi ,t)`$ and $`E(x,\xi ,t)`$, depending on three kinematical variables. In terms of the four-momentum transfer $`\mathrm{\Delta }=P^{}P`$, the invariant momentum transfer squared is $`t=\mathrm{\Delta }^2`$, the skewedness is $`\xi =n\mathrm{\Delta }/2`$, and $`x`$ denotes the momentum fraction. Since the dependence of the GPDs, $`H(x,\xi ,t)`$ and $`E(x,\xi ,t)`$, on three kinematical variables renders their physical interpretation more difficult than ordinary parton distributions, it is useful to recall several important physical properties. In the forward limit, i.e. $`\xi `$, $`t0`$, we recover the forward parton distribution function as $`H(x,0,0)=q(x)`$. In what is sometimes referred to as the local limit, integrating over the momentum fraction, $`x`$, yields the familiar electromagnetic form factors, $`𝑑xH(x,\xi ,t)=F_1(t)`$ and $`𝑑xE(x,\xi ,t)=F_2(t)`$. The first moment of the sum of $`H`$ and $`E`$ yields the total quark angular momentum $`\frac{1}{2}𝑑xx[H(x,\xi ,t)+E(x,\xi ,t)]=J_q`$. We note that both these results are independent of $`\xi `$. And finally, at skewedness $`\xi =0`$, the $`t`$ dependence specifies the transverse Fourier transform with respect to the impact parameter of the light cone wave function Burkardt (2000). For a possible interpretation of GPDs with non-zero $`\xi `$, see Refs. Ralston and Pire (2002); Diehl (2002).

On the lattice, instead of matrix elements of the light cone operator, Eq. (1), one again calculates non-diagonal matrix elements of the local operators, Eq. (2), yielding moments of the generalized parton distributions. Following the notation of Ref. Ji (1997), the non-diagonal matrix element $`P^{}|𝒪^{\{\mu _1\mathrm{}\mu _n\}}|P`$ may be expressed in terms of a set of generalized form factors (GFFs) $`A_{ni}(t)`$, $`B_{ni}(t)`$, and $`C_n(t)`$. The form factors $`A_{ni}(t)`$ and $`C_n(t)`$ multiplied by powers of $`\xi ^2`$ yield the moment $`H_n(\xi ,t)=𝑑xx^{n1}H(x,\xi ,t)`$ and the form factors $`B_{ni}(t)`$ and $`C_n(t)`$ multiplied by powers of $`\xi ^2`$ yield the moment $`E_n(\xi ,t)=𝑑xx^{n1}E(x,\xi ,t)`$.

The lowest three moments considered in this work are

$`P^{}|𝒪^{\mu _1}|P`$ $`=`$ $`\gamma ^{\mu _1}A_{10}(t)`$
$`+{\displaystyle \frac{\text{i}}{2m}}\sigma ^{\mu _1\alpha }\mathrm{\Delta }_\alpha B_{10}(t),`$
$`P^{}|𝒪^{\{\mu _1\mu _2\}}|P`$ $`=`$ $`\overline{P}^{\{\mu _1}\gamma ^{\mu _2\}}A_{20}(t)`$
$`+{\displaystyle \frac{\text{i}}{2m}}\overline{P}^{\{\mu _1}\sigma ^{\mu _2\}\alpha }\mathrm{\Delta }_\alpha B_{20}(t)`$
$`+{\displaystyle \frac{1}{m}}\mathrm{\Delta }^{\{\mu _1}\mathrm{\Delta }^{\mu _2\}}1C_2(t),`$
$`P^{}|𝒪^{\{\mu _1\mu _2\mu _3\}}|P`$ $`=`$ $`\overline{P}^{\{\mu _1}\overline{P}^{\mu _2}\gamma ^{\mu _3\}}A_{30}(t)`$ (3)
$`+{\displaystyle \frac{\text{i}}{2m}}\overline{P}^{\{\mu _1}\overline{P}^{\mu _2}\sigma ^{\mu _3\}\alpha }\mathrm{\Delta }_\alpha B_{30}(t)`$
$`+\mathrm{\Delta }^{\{\mu _1}\mathrm{\Delta }^{\mu _2}\gamma ^{\mu _3\}}A_{32}(t)`$
$`+{\displaystyle \frac{\text{i}}{2m}}\mathrm{\Delta }^{\{\mu _1}\mathrm{\Delta }^{\mu _2}\sigma ^{\mu _3\}\alpha }\mathrm{\Delta }_\alpha B_{32}(t),`$
where $`\overline{P}_\mu =(P_\mu +P_\mu ^{})/2`$ and $`\mathrm{\Gamma }=\overline{U}(P^{})\mathrm{\Gamma }U(P)`$. The GFFs, $`A_{ni}(t)`$, $`B_{ni}(t)`$, and $`C_n(t)`$, specify all the information about spin-independent generalized parton distributions that is known to be accessible on the lattice. The limits for $`H`$ and $`E`$ discussed above may be re-expressed in terms of the generalized form factors. The limit $`t0`$ of $`A_{n0}`$ is the familiar parton distribution moment, $`A_{n0}(0)=𝑑xx^{n1}q(x)`$. The electromagnetic form factors are given by $`A_{10}(t)=F_1(t)`$ and $`B_{10}(t)=F_2(t)`$ for the appropriate flavor combination. Finally, the total quark angular momentum is given by the sum of $`A_{20}+B_{20}`$ as $`t0`$, $`J_q=\frac{1}{2}[A_{20}(0)+B_{20}(0)]`$.

In the context of this brief review, we may now consider the challenges and opportunities in calculating moments of generalized parton distributions on the lattice, and compare them with the analogous issues for ordinary parton distributions. Experimentally, three decades of deep inelastic scattering experiments have provided impressive phenomenological determinations of parton distributions as a function of momentum fraction Lai et al. (1997); Glück et al. (1998); Martin et al. (2002). Hence, for parton distributions, the key issues are developing the lattice technology to the point of attaining quantitative agreement with moments of experimental parton distributions and using the lattice as a tool to obtain insight into how these distributions arise from QCD. The present status is that computational limitations restrict unquenched QCD calculations to the heavy quark regime in which the pion mass is heavier than roughly $`500`$MeV and naive linear extrapolation to the physical pion mass yields serious disagreement with experiment. For example, linear extrapolation of the quark momentum fraction exceeds experiment by the order of fifty percent Göckeler et al. (1996a); Best et al. (1997); Göckeler et al. (1997, 2001); Dolgov et al. (2002). Although there are strong indications that physical extrapolation to the chiral limit may introduce corrections of the required magnitude Detmold et al. (2001), there is presently no quantitative theory for the intermediate mass regime and we must exploit the emerging generation of computers to perform the requisite calculations sufficiently close to the chiral regime.

For generalized parton distributions, the situation is quite different. Since GPDs depend on three variables and experimental quantities involve convolutions, there is no prospect of measuring the full dependence on $`x`$, $`\xi `$, and $`t`$. Without additional input arising from first principles, extraction of GPDs from experiments like deeply virtual Compton scattering will necessarily be contaminated by uncontrolled assumptions. Hence, once computer power is sufficient to obtain quantitative agreement with moments of parton distributions, lattice calculations of moments of GPDs will become an essential tool to be used in conjunction with experiment to extract and understand the full dependence on $`x`$, $`\xi `$, and $`t`$. It is thus imperative to develop techniques to calculate these moments.

In addition, theorists may also obtain insight into how QCD works by studying the dependence of hadron structure on the quark mass. This study can begin immediately, addressing the behavior of hadrons in a world where the pion weighs more than $`500`$MeV. This heavy pion world is much closer to the non-relativistic quark model, and as we eventually lower the pion mass, we will learn how QCD evolves from the world of heavy quarks to the physical world of light quarks. Even in the heavy pion world, one can test contemporary assumptions, such as factorization of the $`t`$ dependence Belitsky et al. (2002).

On the computational side, it is essential to confront the additional challenges that arise for GPDs relative to the ordinary parton distributions. Already for forward parton distributions, the tower of operators, Eq. (2), involves operators that become increasingly subject to statistical noise as one progresses to higher and higher derivatives. For GPDs, however, we compound the noise of these operators with the additional noise from the finite momentum transfer, $`\mathrm{\Delta }`$. Hence, in this work, we address the problem of imposing the maximal statistical constraints a lattice calculation can provide on the form factors $`A_{ni}`$, $`B_{ni}`$, and $`C_n`$ appearing in Eq. (I). We note that we have at our disposal the choice of several alternative representations of the hypercubic group on the lattice corresponding to the same continuum operator $`𝒪^{\{\mu _1\mu _2\mathrm{}\mu _n\}}`$ as well as the choice of different kinematic variables corresponding to the same $`t`$. We will therefore use this freedom to construct an over-determined set of lattice observables corresponding to the continuum expressions in Eq. (I), and thereby significantly improve the measurement of the form factors. In section II, we will describe the details of the method. Section III will compare the results of our method with a conventional analysis and present results for the the angular momentum carried by quarks, $`J_q`$. The conclusions will be presented in the final section and an Appendix presents the necessary detailed expressions for matrix elements in terms of generalized form factors.

## II Lattice Calculation of Generalized Form Factors

We will calculate the matrix elements of $`P^{}|𝒪_q^{\{\mu _1\mu _2\mathrm{}\mu _n\}}|P`$ to extract the generalized form factors $`A_{ni}^q(t)`$, $`B_{ni}^q(t)`$, and $`C_n^q(t)`$ where, when relevant, we append a quark flavor label $`q`$. Since the calculation of disconnected diagrams raises yet another level of complexity, in this present work we will restrict our attention primarily to the flavor non-singlet combination $`ud`$, for which these diagrams cancel.

In the usual way, we calculate the non-diagonal matrix elements by the following ratio of three- and two-point functions,

$$R_𝒪(\tau ,P^{},P)=\frac{C_𝒪^{\text{3pt}}(\tau ,P^{},P)}{C^{\text{2pt}}(\tau _{\text{snk}},P^{})}\left[\frac{C^{\text{2pt}}(\tau _{\text{snk}}\tau +\tau _{\text{src}},P)C^{\text{2pt}}(\tau ,P^{})C^{\text{2pt}}(\tau _{\text{snk}},P^{})}{C^{\text{2pt}}(\tau _{\text{snk}}\tau +\tau _{\text{src}},P^{})C^{\text{2pt}}(\tau ,P)C^{\text{2pt}}(\tau _{\text{snk}},P)}\right]^{1/2}.$$
(4)
The factors relating this ratio to the physical continuum matrix element are given in Eq. (25). The correlation functions are given by

$`C^{\text{2pt}}(\tau ,P)`$ $`=`$ $`{\displaystyle \underset{j,k}{}}\left(\mathrm{\Gamma }_{\text{unpol}}\right)_{jk}\mathrm{\Omega }\left|N_k(\tau ,P)\overline{N}_j(\tau _{\text{src}},P)\right|\mathrm{\Omega },`$
$`C_𝒪^{\text{3pt}}(\tau ,P^{},P)`$ $`=`$ $`{\displaystyle \underset{j,k}{}}\left(\mathrm{\Gamma }_{\text{pol}}\right)_{jk}\mathrm{\Omega }\left|N_k(\tau _{\text{snk}},P^{})𝒪(\tau )\overline{N}_j(\tau _{\text{src}},P)\right|\mathrm{\Omega },`$ (5)
where $`|\mathrm{\Omega }`$ denotes the QCD vacuum state. The nucleon source, $`\overline{N}(\tau ,P)`$, and sink, $`N(\tau ,P)`$, create and annihilate states with the quantum numbers of the nucleon, and to maximize the overlap with the ground state, we used the smeared sources defined in Ref. Dolgov et al. (2002). The source is located at timeslice $`\tau _{\text{src}}`$, the operator $`𝒪`$ is inserted at timeslice $`\tau `$, and the sink is positioned at timeslice $`\tau _{\text{snk}}`$. Explicit expressions for the polarized and unpolarized projectors $`\mathrm{\Gamma }_{\text{pol/unpol}}`$ are given in Appendix A, Eqs. (26,27).

Inserting a complete set of states into Eq. (II) and using the time evolution operator yields

$`C^{\text{2pt}}(\tau ,P)`$ $`=`$ $`{\displaystyle \underset{l}{}}e^{E_l(P)(\tau \tau _{\text{src}})}\text{Tr}\left(\mathrm{\Gamma }_{\text{unpol}}\mathrm{\Omega }|N(\tau ,P)|ll|\overline{N}(\tau _{\text{src}},P)|\mathrm{\Omega }\right)`$ (6)
$`=`$ $`e^{E_0(P)(\tau \tau _{\text{src}})}{\displaystyle \frac{\left(Z(P)\overline{Z}(P)\right)^{1/2}}{E_0(P)}}\text{Tr}\left(\mathrm{\Gamma }_{\text{unpol}}U(P)\overline{U}(P)\right)`$
$`+\text{higher states},`$
$`C_𝒪^{\text{3pt}}(\tau ,P^{},P)`$ $`=`$ $`{\displaystyle \underset{k,l}{}}\text{Tr}\{\mathrm{\Gamma }_{\text{pol}}\mathrm{\Omega }|N(\tau _{\text{snk}},P^{})|kk|e^{E_k(P^{})\left(\tau _{\text{snk}}\tau \right)}`$ (7)
$`\times 𝒪|ll|e^{E_l(P)\left(\tau \tau _{\text{src}}\right)}|\overline{N}(\tau _{\text{src}},P)|\mathrm{\Omega }\}`$
$`=`$ $`e^{E_0(P)(\tau \tau _{\text{src}})E_0(P^{})(\tau _{\text{snk}}\tau )}{\displaystyle \frac{\left(Z(P)\overline{Z}(P^{})\right)^{1/2}}{E_0(P^{})E_0(P)}}\text{Tr}\left(\mathrm{\Gamma }_{\text{pol}}U(P^{})\overline{U}(P)\right)`$
$`\times P^{}\left|𝒪(\tau )\right|P+\text{higher states}.`$
The contributions from higher states in Eqns. (6) and (7) are suppressed by exponential prefactors when $`\tau _{\text{snk}}\tau `$ and $`\tau \tau _{\text{src}}`$ are significantly greater than the inverse of the excitation energy of the first excited state.

The ratio, Eq. (4), is constructed to exactly cancel all exponential and wave-function overlap factors. The two-point functions $`C^{\text{2pt}}(\tau _{\text{snk}},P^{})`$ and $`C^{\text{2pt}}(\tau _{\text{snk}},P)`$ decay exponentially for the full Euclidean distance between the source and sink and are thus particularly subject to statistical noise with finite statistics. In the worst case, they may even become negative, and these cases are excluded from the present work. We note that other possibilities besides Eq. 4 may be used to cancel the exponential and overlap factors, and this freedom will be explored in a subsequent work.

For sufficiently large time separations the ratio $`R(\tau ,P^{},P)`$ will exhibit a plateau yielding the desired lattice matrix element, and the plateau value, $`\overline{R}(P^{},P)`$, is obtained by averaging over an appropriate range of timeslices, $`\tau _{\text{min}}`$ to $`\tau _{\text{max}}`$,

$$\overline{R}_{\{\mu _1\mu _2\mathrm{}\mu _n\}}(P^{},P)=\frac{1}{\tau _{\text{max}}\tau _{\text{min}}}\underset{\tau =\tau _{\text{min}}}{\overset{\tau _{\text{max}}}{}}R_{\{\mu _1\mu _2\mathrm{}\mu _n\}}(\tau ,P^{},P).$$
(8)
To convert our lattice calculations to the continuum $`\overline{\text{MS}}`$-scheme, we use 1-loop perturbative matching at the scale $`\mu ^2=4\text{GeV}^2`$

$$𝒪_i^{\overline{\text{MS}}}(\mu )=\underset{j}{}Z_{ij}(\mu ,a)𝒪_j^{\text{lat}}(a),$$
(9)
so that the lattice matrix element is related to the continuum matrix element by

$$P^{}\left|𝒪_i^{\overline{\text{MS}}}\right|P=\sqrt{E(P^{})E(P)}\underset{j}{}Z_{ij}\overline{R}_j.$$
(10)
Note that the renormalization constant, $`Z_𝒪(\mu ,a)`$, depends only on the operator $`𝒪`$, but not on the external states, so that Eq. (10) is valid for any external momenta, $`P^{}`$ and $`P`$.

Finally, we write the Euclidean continuum relation between the renormalized matrix element of the generalized current $`P^{}|𝒪_{\{\mu _1\mu _2\mathrm{}\mu _n\}}^q|P`$ and the desired generalized form factors $`A_{ni}^q(t),B_{nj}^q(t)`$, and $`C_n^q(t)`$ in the following abbreviated notation:

$$P^{}\left|𝒪_{\{\mu _1\mu _2\mathrm{}\mu _n\}}^q\right|P=\underset{i}{}a_iA_{ni}^q+\underset{j}{}b_jB_{nj}^q+cC_n^q.$$
(11)
Full expressions for the kinematic factors $`\{a_i,b_j,c\}`$ are given in Appendix A for $`n=1,2,`$ and $`3`$, Eq. (28).

For a given $`n`$, we may evaluate the $`N_n^{GFF}`$ generalized form factors $`\{A_{ni}^q(t),B_{nj}^q(t),C_n^q(t)\}`$ as follows. We select $`NN_n^{GFF}`$ sets of operators $`𝒪_{\{\mu _1\mu _2\mathrm{}\mu _n\}}^q`$ and momenta $`\{P^{},P\}`$ such that Eq. (11) specifies $`N_n^{GFF}`$ linearly independent combinations of the form factors. Lattice matrix elements for these operators and momenta are calculated from the the ratios $`R(\tau ,P^{},P)`$, in Eq. (4) and matched to continuum operators via Eq. (9). If $`N=N_n^{GFF}`$, the GFFs are calculated by inverting Eq. (11), and if $`N>N_n^{GFF}`$, they are calculated by a least-squares fit to the overdetermined system Eq. (11).

Note, that in contrast to the case of forward parton distributions, where the moments correspond to a single number (denoted $`v_n`$ or $`x^{n1}`$), or electromagnetic form factors, where there are two form factors ($`F_1`$ and $`F_2`$), we have the complication of $`N_n^{GFF}`$ unknown generalized form factors to be determined. We therefore now discuss the strategy for selecting an appropriate set of operators and momenta for this task.

### II.1 Practical considerations

One practical concern is the numerical noise associated with momentum projection. The three-point function is subject to noise from the projection of the sink onto momentum $`P^{}`$ and of the operator onto momentum transfer $`\mathrm{\Delta }`$. In addition to being subject to the sink momentum projection, the two-point functions appearing in the lattice ratio also collect noise from the projection onto the source momentum $`P`$. Ideally, for each invariant momentum transfer $`t`$, one would like to select $`P`$ and $`P^{}`$ such as to minimize these errors, but each distinct choice of sink momentum $`P^{}`$ requires an expensive calculation of a new set of backward propagators, whereas changing $`\mathrm{\Delta }`$ for fixed $`P^{}`$ requires no new propagators. Hence, denoting the lowest momentum attained on a lattice with $`L_s`$ lattice points in a spatial dimension as $`p_l=2\pi /L_s`$. As a practical matter, we have used two sink momenta, $`\stackrel{}{P}^{}=(0,0,0)`$ and $`\stackrel{}{P}^{}=(p_l,0,0)`$. With these sinks, we generate a substantial set of momenta listed in Appendix A, Tab. 2, by using a range of values of $`\mathrm{\Delta }`$. However there are two special virtualities of note, the virtuality, $`t_{\text{vlow}}`$ and the Breit frame virtuality $`t_{\text{Breit}}`$ given by the following momentum combinations, or rotations thereof:

$`\stackrel{}{P}^{}=(0,0,0),\stackrel{}{P}=(p_l,0,0)`$ $``$ $`t_{\text{vlow}}=\left(m\sqrt{p_l^2+m^2}\right)^2p_l^2,`$ (12)
$`\stackrel{}{P}^{}=(p_l,0,0),\stackrel{}{P}=\stackrel{}{P}^{}`$ $``$ $`t_{\text{Breit}}=4p_l^2,`$ (13)
The virtuality $`t_{\text{vlow}}`$ is the lowest virtuality that can be placed on the lattice for the set of external momenta we use, for which no spatial momentum exceeds $`p_l`$, and would be expected to have minimal projection error. Since the time component does not reduce the virtuality in the Breit frame, the Breit frame provides the optimal means of providing a large virtuality with minimal momentum projection error. Hence, including both $`t_{\text{vlow}}`$ and $`t_{\text{Breit}}`$ in our full set of momentum selections ensures the presence of measurements at both ends of our $`t`$-range that have the minimal possible momentum-induced statistical error. In addition to these considerations, for certain different virtualities $`t`$, there are more available momentum combinations and thus more constraints than in the Breit frame. Hence, we obtain comparably small errors for most of the other virtualities as well.

Another practical issue in measuring generalized form factors is the presence of powers of the momentum transfer $`\mathrm{\Delta }`$ in the kinematical factors in Eq. (11). The four-momentum $`\mathrm{\Delta }`$ is — compared to the time component of external momenta in the cases we consider — a small number, so its presence in kinematical factors amplifies the effect of statistical errors in lattice matrix elements on measurements of the associated form factors. Consistent with the pattern observed in Eq. I, in the general case the factors $`a_i`$ multiplying $`A_{ni}^q(t)`$ contain $`i=\{0,2,\mathrm{}2[\frac{n1}{2}]\}`$ factors of $`\mathrm{\Delta }`$, the factors $`b_j`$ multiplying $`B_{nj}^q(t)`$ contain $`j=\{1,3,\mathrm{}1+2[\frac{n1}{2}]\}`$ factors of $`\mathrm{\Delta }`$, and the factor $`c`$, which is non-vanishing only for even $`n`$, has $`n`$ factors of $`\mathrm{\Delta }`$. Therefore, we expect $`A_{n0}^q(t)`$ to be the quantity that can be extracted most accurately from lattice calculations, $`B_{n0}^q(t)`$ should have slightly larger errors, $`C_n^q(t)`$ is the worst determined when it is non-vanishing, and the remaining GFFs, $`A_{ni}^q(t)`$ and $`B_{nj}^q(t)`$, should lie somewhere in between these extremes.

### II.2 Overdetermined set of lattice observables

To extract generalized form factors from lattice calculations, we now consider a fixed value of the virtuality, $`t`$, that can be achieved on the lattice, and abbreviate Eqns. (10) and (11) in the following schematic form:

$`𝒪_i^{cont}`$ $`=`$ $`{\displaystyle \underset{j}{}}a_{ij}_j,`$
$`𝒪_i^{cont}`$ $`=`$ $`\sqrt{E^{}E}{\displaystyle \underset{j}{}}Z_{ij}\overline{R}_j,`$ (14)
where we denote the generalized form factors generically by $`_j(t)`$ and $`j`$ runs over all the form factor labels for the $`n`$ under consideration. Thus, for n = 1, $`_1(t)=A_{10}^q(t)`$ and $`_2(t)=A_{20}^q(t)`$. Note that the familiar spin non-flip and spin-flip electromagnetic form factors are defined by these flavor form factors weighted with the electric charge, i.e. $`F_1(t)=2/3A_{10}^\text{u}(t)1/3A_{10}^\text{d}(t)`$, and $`F_2(t)=2/3B_{10}^\text{u}(t)1/3B_{10}^\text{d}(t)`$. Similarly, for $`n=2`$, corresponding to the operator $`𝒪_{\{\mu \nu \}}`$, the $`_i(t)`$ in Eq. (II.2) denote $`_1(t)=A_{20}^q(t)`$, $`_2(t)=B_{20}^q(t)`$, and $`_3(t)=C_2^q(t)`$.

Eliminating $`𝒪_i^{cont}`$ from Eqs. (II.2) yields the final form of the set of equations we will use for our analysis:

$`\overline{R}_i`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{E^{}E}}}{\displaystyle \underset{jk}{}}Z_{ij}^1a_{jk}_k`$ (15)
$``$ $`{\displaystyle \underset{j}{}}a_{ij}^{}_j.`$
In this set of equations, the label $`j`$ runs over the $`N_n^{GFF}`$ generalized form factors for the given $`n`$. The label $`i`$ denotes a specific choice of nucleon momenta $`P`$ and $`P^{}`$, and an index combination for a lattice H(4) representation of a continuum operator $`𝒪_{\{\mu _1\mu _2\mathrm{}\mu _n\}}^q`$. We consider $`N`$ distinct values of this label, $`i`$, corresponding to the same $`n`$ and $`t`$ and require $`N>N_n^{GFF}`$. Operationally, to obtain this collection of $`N`$ equations, we list all external lattice momenta, group them together in classes of identical virtual momentum transfer, and write down for each momentum combination the set of indices, $`\{i\}`$, of $`𝒪_{\{i\}}`$ that gives rise to non-trivial kinematic factors, $`a_{ij}`$. Note that sets of external momenta in which the spatial components differ only by cubic rotations still yield different constraints. Rotations changing components along the spin polarization axis in general change the kinematic factors and yield independent equations. Rotations that change components orthogonal to the spin polarization axis yield constraints that are physically equivalent, but provide additional statistical constraints in a lattice Monte Carlo calculation. Appendix A shows the explicit momentum combinations and the basis of diagonal index combinations we have chosen in this work.

Equation (15) is an overdetermined system of linear equations that in principle provides the maximal information on the generalized form factors that is attainable from the lattice. We will compare it in Sec. III.1 with the conventional technique of using a uniquely determined system of equations. In the case of measuring moments of parton distributions (see, e.g. Refs. Göckeler et al. (1996a); Best et al. (1997); Göckeler et al. (1997, 2001); Dolgov et al. (2002)) the conventional technique corresponds to selecting a single row of Eq. (15), and for electromagnetic form factors, it corresponds to taking linear combinations or selecting kinematic variables such that $`F_1`$ and $`F_2`$ appear separately in two equations.

It is important to note that Eq. (15) contains systematic, as well as statistical errors. The lattice action, the discrete approximations to the operators $`𝒪_{\{\mu _1\mu _2\mathrm{}\mu _n\}}^q`$, and the nucleon dispersion relation contain lattice artifacts, and since we use Wilson fermions, the calculation contains errors of $`𝒪(a)`$. In principle these lattice artifacts can by addressed by improved operators and extrapolation to the continuum limit, but these corrections are beyond the scope of the present work. In addition, the use of 1-loop perturbative renormalization introduces systematic errors in the analysis, which ultimately can be improved by higher order corrections or non-perturbative renormalization.

### II.3 Singular Value Decomposition

We now turn to the solution of the overdetermined set of equations, Eq. (15), which we write in matrix form as follows:

$$R^{}=A^{}.$$
(16)
Here, $``$ is the vector of the desired GFFs, $`R^{}`$ denotes the statistically measured lattice ratios, and $`A^{}`$ is the matrix of the coefficients $`a_{ij}^{}`$. Changing to a standard notation, we let $`n+1`$ denote the number of GFFs, which we previously called $`N^{GFF}`$, and continue to call the number of constraints $`N`$, so that $`A`$ is an $`N\times (n+1)`$ matrix, $``$ is a vector of length $`n+1`$, and $`R^{}`$ a vector of length $`N`$.

We solve this overdetermined problem by minimization of the $`\chi ^2`$ norm

$$\chi ^2=\underset{i=1}{\overset{N}{}}\left(\frac{_{j=1}^{n+1}A_{ij}^{}_jR_i^{}}{\sigma _i}\right)^2,$$
(17)
where $`\sigma _i`$ denotes the jackknife error of the lattice measurement of $`R_i^{}`$. In this case the jackknife procedure should be performed with an appropriate bin-size to eliminate autocorrelation effects. It is convenient to absorb the errors $`\sigma _i`$ by defining $`A_{ij}=A_{ij}^{}/\sigma _i`$ and $`R_i=R_i^{}/\sigma _i`$ so that we seek a vector $``$ that minimizes

$$\chi ^2=|AR|^2.$$
(18)
The method of choice to solve Eq. (18) in the presence of singularities or near singularities is Singular Value Decomposition (SVD) Teukolsky et al. (1992). It is based on the theorem that any $`N\times (n+1)`$ matrix may be decomposed as

$$A=U\text{diag}(w_1,\mathrm{},w_{n+1})V,$$

where $`U`$ is an $`N\times (n+1)`$ matrix, $`V`$ is $`(n+1)\times (n+1)`$, and all $`w_i`$, the singular values of $`A`$, are non-negative. If all the $`w_i`$ are non-zero, the solution that minimizes Eq. (18) is given by

$$F=V\text{diag}(1/w_1,\mathrm{},1/w_{n+1})\left(U^TR\right),$$
(19)
which provides the optimal vector of form factors, $``$, we can achieve with all the information available.

If none of the singular values is zero, the method is equivalent to standard least-squares minimization. If one or more of the singular values, $`w_i`$, are zero, the rank of the matrix $`A`$ is reduced by the number of zero singular values, implying that the associated directions in the solution space can not be explored. In the case of zero or nearly zero singular values, we may opt to avoid exploring the offending direction(s) and minimize Eq. (18) in the residual subspace by replacing the corresponding factors $`1/w_i`$ by zero. In this case, of course, only the components of $``$ corresponding to non-zero values of $`w_i`$ are determined.

As an example, in the forward case the singular values for all GFFs other than $`A_{n0}^q(0)`$ are zero and the result is the minimization in the subspace of $`A_{n0}^q(0)`$ alone. In this case, $`A_{n0}^q(0`$) is still determined simultaneously from several different lattice measurements and this procedure thereby automatically incorporates more physical constraints than the conventional method of fitting a single lattice observable.

The errors in $``$ are calculated using a jackknife analysis in the following way. The system (18) is solved simultaneously for the full set of lattice measurements $`R`$ to get the average solution vector $``$. Then, the appropriate subsamples $`R_l`$ are formed from which the subsample solutions $`_l`$ are calculated. The error vector for the solution, $`\sigma _{}`$, is formed from the $`_l`$s by the familiar jackknife formula (see e.g. Schroers (2001) and references therein for details).

In summary, the method introduced in this section provides a way to extract maximal information from a set of lattice calculations. It simultaneously includes all available data in a unified manner and corrects for correlations present in the underlying sample of gauge field configurations.

## III Results

In this section, we apply our procedure to the case $`n=2`$ using the operators in Eqs. (A.1,30) and thoroughly analyze the results. For this exploratory calculation, we use SESAM unimproved Wilson configurations on a $`16^3\times 32`$ lattice at $`\beta =5.6`$ with $`\kappa _{\text{sea}}=\kappa _{\text{val}}=0.1560`$. We note that this is the heaviest of the three quark masses used in earlier calculations of moments of quark distributionsDolgov et al. (2002), and using the lattice spacing $`a^1=2.01`$GeV from the chirally extrapolated nucleon mass, the pion mass $`(am)_\pi =0.446(3)`$ corresponds to $`m_\pi =896(6)`$MeV. All the computational details for calculating smeared sources and sinks, nonrelativistic sequential sources, the difference approximations to the operators $`𝒪_{\{\mu _1\mu _2\mathrm{}\mu _n\}}^q`$, spin polarization projection, and Dirichlet boundary conditions are as in Ref. Dolgov et al. (2002). In this present work, we relabel the time slices such that the source is at $`\tau _{\text{src}}=1`$, the sink is at $`\tau _{\text{snk}}=13`$ and boundaries are at $`\tau =9`$ and $`\tau =22`$.

The one-loop perturbative renormalization constants for the cases we consider in this work are taken from Ref. Dolgov et al. (2002). Note, that in the case of the operators $`𝒪_{\{\mu \nu \}}`$, Eqs. (A.1,30), the renormalization constant will depend on the index combination chosen. The diagonal and non-diagonal index combinations belong to distinct irreducible representations of the lattice point group and have different renormalization factors: $`Z_{\text{diag}}^{n=2}(\mu ^2=4\text{GeV}^2)=0.9768`$ and $`Z_{\text{non-diag}}^{n=2}(\mu =4\text{GeV}^2)=0.9884`$.

### III.1 Comparison of Minimally Determined and Overdetermined Measurements

To explore how the method works and to compare it with the conventional method of measuring a single operator for each observable, we consider the case of virtuality $`t_{\text{vlow}}=0.5925`$GeV<sup>2</sup>, which may be achieved by choosing any of the momentum combinations displayed in Tab. 1.

We start by considering only the combination of external momenta, $`\stackrel{}{P}^{}=(0,0,0)`$, $`\stackrel{}{P}=(p_l,0,0)`$, listed in the first column in Tab. 1. The system of equations, Eq. (15) then contains the operators

$$P^{}|𝒪_{\text{diag,1}}^{\text{u-d}}|P,P^{}|𝒪_{\text{diag,2}}^{\text{u-d}}|P,P^{}|𝒪_{\text{diag,3}}^{\text{u-d}}|P,P^{}|𝒪_{\{10\}}^{\text{u-d}}|P,P^{}|𝒪_{\{20\}}^{\text{u-d}}|P,P^{}|𝒪_{\{21\}}^{\text{u-d}}|P,$$
(20)
with the operator index conventions given in Appendix A, Eqs. (A.1,30).

Plateau plots of the ratios, $`R(\tau ,P^{},P)`$, are shown in Fig. 1 for three matrix elements, $`P^{}|𝒪_{\text{diag,1}}^{\text{u-d}}|P,P^{}|𝒪_{\text{diag,2}}^{\text{u-d}}|P,`$ and $`P^{}|𝒪_{\{20\}}^{\text{u-d}}|P`$. Measurements are averaged over the interval $`\left[\tau _{\text{min}}=5,\tau _{\text{max}}=9\right]`$ and this interval is denoted by vertical lines in the plots. The solid horizontal lines show the extracted plateau values and the dotted lines denote one standard deviation.

The resulting lattice matrix elements, together with their jackknife errors are as follows:

$`\overline{R}_{\text{diag,1}}^{\text{u-d}}(P^{},P)`$ $`=`$ $`0.2327\pm 0.0224,`$
$`\overline{R}_{\text{diag,2}}^{\text{u-d}}(P^{},P)`$ $`=`$ $`0.2970\pm 0.0235,`$
$`\overline{R}_{\text{diag,3}}^{\text{u-d}}(P^{},P)`$ $`=`$ $`0.00832\pm 0.01651,`$
$`\overline{R}_{\{10\}}^{\text{u-d}}(P^{},P)`$ $`=`$ $`0.09768\pm 0.01648,`$
$`\overline{R}_{\{20\}}^{\text{u-d}}(P^{},P)`$ $`=`$ $`0.1802\pm 0.0205,`$
$`\overline{R}_{\{21\}}^{\text{u-d}}(P^{},P)`$ $`=`$ $`0.05926\pm 0.01162.`$ (21)
We now consider a minimal set of three operators to determine, but not overdetermine, the three generalized form factors, $`A_{20}^{\text{u-d}},B_{20}^{\text{u-d}}`$, and $`C_{20}^{\text{u-d}}`$. The subset of diagonal operator equations alone are insufficient Schroers and QCDSF collaboration (2002), since the diagonal index combinations are linearly dependent and provide only a system of equations of rank two.

The minimal choice to determine, but not overdetermine, the form factors is to choose two diagonal operators and one off-diagonal operator, analogous to familiar procedures to calculate moments of parton distributionsGöckeler et al. (1996a); Dolgov et al. (2002) or electromagnetic form factors Göckeler et al. (2001); Göckeler et al. (2003a). The best choice for this purpose is the set of matrix elements $`P^{}|𝒪_{\text{diag,1}}^{\text{u-d}}|P`$, $`P^{}|𝒪_{\text{diag,2}}^{\text{u-d}}|P`$, and $`P^{}|𝒪_{\{20\}}^{\text{u-d}}|P`$, which have the smallest relative errors and provide a system of three linearly independent equations.

The form factors and associated errors determined this way are shown in Fig. 2 by the triangles plotted to at the left at $`N=0`$.

Still considering the same single choice of external momenta, the full set of six matrix elements in Eq. (20) provides an overdetermined set of equations that yield the form factors and errors labeled $`N=1`$ in Fig 2. We previously argued that the $`C`$ form factor is least accurately determined because its coefficients have the highest powers of $`\mathrm{\Delta }`$, and the minimally determined fit, $`N=0`$ bears out this expectation. Note that the overdetermined fit, $`N=1`$ leaves the errors in $`A`$ and $`B`$ roughly the same, but substantially reduces that of $`C`$. The results of the two fits agree within errors, but since the overdetermined fit exploits the maximal information in the lattice measurements, Eq. (III.1) and has smaller overall statistical errors, we consider it superior.

We now progress to the next level by including all the external momentum combinations enumerated in Tab. 1. The results are displayed in Fig. 2, where the abscissa denotes the total number of momentum combinations included in the fit. Entry $`N`$ corresponds to the inclusion of momentum sets $`1,2,\mathrm{}N`$ of Tab. 1. As expected, the errors for each form factor decrease significantly as new statistical information is included by adding momenta in new directions. The improvement is generally less than $`\sqrt{N}`$ because measurements of different momenta on the same lattices are correlated.

The overall result in Fig. 2 is a strong validation of our approach. Comparing the minimally determined, $`N=0`$ result at the far left with the overdetermined result with all seven momentum combinations, $`N=7`$ at the far right, we observe that our method reduces the error in the least accurately determined form factor $`C`$ by a factor of 5, the error in $`B`$ by a factor of 3, and the error in $`A`$ by a factor of 2.

### III.2 Generalized Form Factors and Quark Angular Momentum

Having validated our analysis technique, we now apply it to the full set of external momenta and virtualities listed in Appendix A, Tab. 2. The results of the full, overdetermined analysis for the flavor non-singlet generalized form factors $`A_{20}^{\text{u-d}}(t)`$, $`B_{20}^{\text{u-d}}(t)`$, and $`C_2^{\text{u-d}}(t)`$ are shown in Fig. 3. We observe that with the present analysis, all three generalized form factors are determined quite accurately throughout the full range of virtuality up to $`t=3.1`$GeV<sup>2</sup>. As expected, some combinations of external momenta induce more statistical noise than others, but the overall structure of the form factors is well determined. Although there is no fundamental argument for the functional dependence on $`t`$, it is useful to fit the form factors to the dipole form that provides a good phenomenological fit to the nucleon electromagnetic form factors, and the results of least-squares fit by dipole form factors are shown for $`A(t)`$ and $`B(t)`$. Although the individual up and down form factors $`C_2^{\text{u/d}}(t)`$ are non-vanishing, their difference, $`C_2^{\text{u-d}}(t)`$, is statistically consistent with zero. We note this flavor independence is a feature of the Chiral Quark Soliton Model in Ref. Kivel et al. (2001).

Although results are not yet available at other quark masses to enable extrapolation to the physical quark mass, these results clearly show the behavior of the generalized form factors in the $`900`$MeV pion world. For this non-singlet case, there are no corrections from disconnected diagrams.

It is useful to point out that in all cases for which $`t0`$, there were no zero (or nearly zero) singular values in our calculations. Hence, in these cases, the singular value decomposition is completely equivalent to minimizing Eq. (18) by conventional least-squares analysis. In the case of $`t=0`$, because of the explicit factors of $`\mathrm{\Delta }`$ appearing in Eq. (I), only the coefficients $`A_{n0}^q(0)`$ can be determined. As expected, the singular value decomposition handled this automatically, with the zero singular values resulting in minimization in the appropriate subspace of $`A_{n0}^q(0)`$. Note that this procedure introduces no bias or model assumptions. Rather, it simply specifies the correct physical subspace.

Since the total quark angular momentum is given by the zero virtuality limit, $`J_q=\frac{1}{2}[A_{20}^{\text{u+d}}(0)+B_{20}^{\text{u+d}}(0)]`$, it is particularly interesting to study this limit. Two issues need to be addressed. The first is extrapolation to $`t=0`$. This is not a problem for $`A`$, since it is just the moment of the spin-averaged parton distribution $`A_{20}(0)=x`$. However, like the electromagnetic form factor $`F_2`$, $`B`$ cannot be measured at $`t=0`$ since the kinematical factors always contain $`\mathrm{\Delta }`$, and therefore it must be extrapolated. From the behavior of the flavor non-singlet combination $`B_{20}^{\text{u-d}}(t)`$ in Fig. 3, one might expect significant extrapolation errors. However, as shown in Fig. 4, the cancellation between the up and down quark contributions is nearly complete and the flavor singlet combination $`B_{20}^{\text{u+d}}(t)`$ is consistent with zero. Hence, essentially the full contribution comes from $`A_{20}^{\text{u+d}}(0)`$.

The second issue is the fact that disconnected diagrams, which are beyond the scope of the present work, must also be included in these flavor singlet matrix elements.

It is interesting to note that light cone arguments indicate that the contribution of each Fock space sector to the gravitomagnetic moment $`B(0)=_{q,g}B_{20}^{q,g}(0)=0`$ is separately zero Brodsky et al. (2001). Although this may not necessarily imply that connected and disconnected diagrams for $`B(0)`$ separately cancel, it is suggestive of our result above that the connected contribution to $`B(0)`$ was consistent with zero. Note that the separate vanishing of quark and gluon contributions to $`B(0)`$ was also conjectured in Ref. Teryaev (1998).

With the caveat that we must omit disconnected diagrams for the present, we proceed to quote our lattice results from connected diagrams. Using the results of Ref. Dolgov et al. (2002) at $`\kappa =0.1560`$, the quark spin contribution to nucleon angular momentum is given by the lowest moment of the longitudinal spin distribution

$`{\displaystyle \frac{1}{2}}\mathrm{\Delta }\mathrm{\Sigma }`$ $`=`$ $`{\displaystyle \frac{1}{2}}[1_{\mathrm{\Delta }u}+1_{\mathrm{\Delta }d}]`$ (22)
$``$ $`{\displaystyle \frac{1}{2}}\mathrm{\hspace{0.17em}\hspace{0.17em}0.682}(18).`$
Since $`B_{20}^{\text{u+d}}(0)0`$, the total quark spin and orbital contribution to the nucleon angular momentum is

$`J_q`$ $`=`$ $`{\displaystyle \frac{1}{2}}[A_{20}^{\text{u+d}}(0)+B_{20}^{\text{u+d}}(0)]`$ (23)
$``$ $`{\displaystyle \frac{1}{2}}[x_u+x_d+0]`$
$``$ $`{\displaystyle \frac{1}{2}}\mathrm{\hspace{0.17em}\hspace{0.17em}0.675}(7).`$
Thus, since these two results are equal within errors, if we consider only connected diagrams in the world of $`900`$MeV pions, we conclude that the quark spin produces 68 percent of the nucleon angular momentum and the quark orbital angular momentum contribution is negligible. Although this is superficially consistent with expectations from relativistic quark models, we note that serious physical interpretation and quantitative comparison with quark models requires consideration of the disconnected diagrams. The striking fact for our present purposes, however, is that it is possible to measure the connected part to the order of one percent, under the assumption that $`B_{20}^{\text{u+d}}(0)=0\pm 0`$.

Note that quenched connected diagram contributions of $`\mathrm{\Delta }\mathrm{\Sigma }`$ and $`J_q`$ have been reported in Refs. Mathur et al. (2000); Göckeler et al. (2003b). Since those calculations were performed for pion masses ranging above $`500`$MeV, we would expect quenching effects to be negligible and hence that their results would be statistically consistent with ours. The chiral extrapolated results of Göckeler et. al. Göckeler et al. (2003b) that 66 $`\pm `$ 14 percent of the nucleon spin arises from the quark spin and that 6 $`\pm `$ 14 percent arises from quark orbital angular momentum are completely compatible with our results because of the small quark mass dependence. The results for $`\mathrm{\Delta }\mathrm{\Sigma }`$ by Mathur et. al. in Ref. Mathur et al. (2000) are also consistent with ours. Reference Dolgov et al. (2002) showed that the linear extrapolation of $`\mathrm{\Delta }\mathrm{\Sigma }`$ is essentially constant, only changing from 68 $`\pm `$ 2 percent of the nucleon spin at $`m_\pi =900`$ MeV to 69 percent in the chiral limit. This is consistent with the connected diagram chiral extrapolation by Mathur et. al. of 62 $`\pm `$ 8 percent. The results for $`J_q`$, however, are inconsistent. From the results of their Fig. 3 and their 1.045 renormalization factor, $`J_q`$ contributes 84 $`\pm `$ 8 percent of the spin, which disagrees with our result of 68 $`\pm `$ 1 percent. We believe the discrepancy arises because Mathur et. al. use the dipole prescription to extrapolate the sum $`A_{20}^{\text{u+d}}+B_{20}^{\text{u+d}}`$, whereas, as shown in our Fig. 4, we calculate $`A_{20}^{\text{u+d}}(0)`$ directly without extrapolation and only extrapolate the nearly vanishing $`B_{20}^{\text{u+d}}`$.

## IV Conclusions

In summary, we have presented a new method for calculating generalized form factors in lattice QCD that extracts the full information content from a given lattice configuration by measuring an overdetermined set of lattice observables. We demonstrated its effectiveness in an exploratory calculation of $`n=2`$ generalized form factors up to $`3`$GeV<sup>2</sup> and showed that it reduces errors to as small as one-fifth of those obtained by a conventional minimally-determined analysis. The final error bars for $`A_{20}^{\text{u-d}}(t)`$, $`A_{20}^{\text{u+d}}(t)`$, and $`B_{20}^{\text{u-d}}(t)`$ are typically of the order of five to ten percent, providing useful information about the Fourier transform of the transverse structure of the nucleon. Because of fortunate cancellation in the flavor singlet case, the connected diagram contributions to the total quark angular momentum are measured to the order of one percent. Although this exploratory calculation was performed for the heaviest of the three SESAM quark masses used in earlier calculation of moments of parton distributions, we expect the technique to be sufficiently robust to treat all the masses.

We note that the first calculation of generalized form factors were reported by Schroers in Schroers and QCDSF collaboration (2002) using the standard minimally-determined analysis, and more extensive calculations by the QCDSF collaboration using our method are being published simultaneously with this present work Göckeler et al. (2003b).

This work provides the foundation for a number of promising investigations. The $`n=3`$ relations introduced in this work enable us to calculate the $`t`$-dependence for three moments and thus explore the variation of the transverse structure of the proton with $`x`$. The same methodology can, of course, be used to calculate spin-dependent generalized form factors. In the longer term, extensive calculations with emerging multi-Teraflops computers dedicated to lattice QCD will enable the chiral and continuum calculations required to have definitive impact on the new generation of experiments currently being undertaken Airapetian et al. (2001); Stepanyan et al. (2001); Favart (2002).

###### Acknowledgements.

The authors wish to thank Matthias Burkardt, Meinulf Göckeler, Axel Kirchner, Andreas Schäfer, Oleg V. Teryaev and Christian Weiss for stimulating discussions. P.H. and W.S. are grateful for Feodor-Lynen Fellowships from the Alexander von Humboldt Foundation and thank the Center for Theoretical Physics at MIT for its hospitality. This work is supported in part by the U.S. Department of Energy (D.O.E.) under cooperative research agreement #DF-FC02-94ER40818.

## Appendix A Matrix elements in terms of generalized form factors

In order to find the relation between Minkowski-expressions, e.g. $`\left(𝒪^\text{M}\right)^{\{\mu _1\mu _2\mathrm{}\mu _n\}}`$, and the corresponding Euclidean quantities, $`\left(𝒪^\text{E}\right)_{\{\mu _1\mu _2\mathrm{}\mu _n\}}`$, we observe that the operators under consideration are constructed from gamma matrices and covariant derivatives $`D`$. The conventions we use, consistent with Ref. Montvay and Münster (1994), are

$`\left(\gamma ^\text{E}\right)_\mu `$ $`=`$ $`{\displaystyle \underset{\nu =1}{\overset{4}{}}}h_{\mu \nu }\left(\gamma ^\text{M}\right)^\nu ,h_{\mu \nu }=\text{diag}(\text{i},\text{i},\text{i},1),`$
$`\left(D^\text{E}\right)_\mu `$ $`=`$ $`{\displaystyle \underset{\nu =1}{\overset{4}{}}}d_{\mu \nu }\left(D^\text{M}\right)^\nu ,d_{\mu \nu }=\text{diag}(1,1,1,\text{i}),`$
and we will set $`\left(a^\text{M}\right)^{\nu =0}=`$ $`\left(a^\text{M}\right)^{\nu =4}`$. Since all operators are symmetrized, the general transformation rule for the operators with one gamma matrix and $`n1`$ covariant derivatives investigated in this work is given by

$$\left(𝒪^\text{E}\right)_{\{\mu _1\mu _2\mathrm{}\mu _n\}}=h_{\mu _1\nu _1}d_{\mu _2\nu _2}\mathrm{}d_{\mu _n\nu _n}\left(𝒪^\text{M}\right)^{\{\nu _1\nu _2\mathrm{}\nu _n\}},$$
(24)
where a summation over the indices $`\nu _i`$ is implicit. In the following, we will drop the label E or M, and it will be clear from the context if Euclidean or Minkowski expressions are used.

We write matrix elements in terms of Dirac spinors as follows

$$P^{}\left|𝒪\right|P=\overline{U}(P^{})𝒦_𝒪(P^{},P)U(P).$$

The ratio $`R_𝒪^{\text{cont}}(P^{},P)`$ can be written explicitly in Minkowski space as

$`R_𝒪^{\text{cont}}(P^{},P)|_{\text{ground state}}`$ $`=`$ $`\left(E(P^{})E(P)\right)^{1/2}`$ (25)
$`\times \left({\displaystyle \frac{1}{2}}\text{Tr}\left[\mathrm{\Gamma }_{\text{unpol}}\left(\mathit{}^{}+m\right)\right]{\displaystyle \frac{1}{2}}\text{Tr}\left[\mathrm{\Gamma }_{\text{unpol}}\left(\mathit{}+m\right)\right]\right)^{1/2}`$
$`\times {\displaystyle \frac{1}{4}}\text{Tr}\left[\mathrm{\Gamma }_{\text{pol}}\left(\mathit{}^{}+m\right)𝒦_𝒪(P^{},P)\left(\mathit{}+m\right)\right],`$
neglecting all excited states which would introduce a time dependence. Note that $`R_𝒪^{\text{cont}}(P^{},P)`$ is directly proportional to $`P^{}\left|𝒪\right|P`$. The ratio in Eq. (25) is expressed in a general form that applies to all quark bilinear operators and projectors $`\mathrm{\Gamma }_{\text{pol/unpol}}`$ relevant to this work.

The Minkowski space expressions for the QCDSF and the LHPC projectors used in our analysis are as follows

$`\mathrm{\Gamma }_{\text{unpol}}^{\text{LHPC}}`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(1+\gamma _0\right),\mathrm{\Gamma }_{\text{pol}}^{\text{LHPC}}={\displaystyle \frac{1}{4}}\left(1+\gamma _0\right)\left(1\gamma _5\gamma _3\right),`$ (26)
$`\mathrm{\Gamma }_{\text{unpol}}^{\text{QCDSF}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(1+\gamma _0\right),\mathrm{\Gamma }_{\text{pol}}^{\text{QCDSF}}={\displaystyle \frac{1}{2}}\left(1+\gamma _0\right)\gamma _5\gamma _2.`$ (27)
We also record the complete expressions for the kernels $`𝒦_𝒪(P^{},P)`$ in Eq. (25). Following the general parameterization of Ji (1997), we have for the lowest three moments $`n=1,2,3`$ in Minkowski space

$`P^{}\left|𝒪^{\{\mu _1\mu _2\mathrm{}\}}\right|P`$ $`=`$ $`\overline{U}(P^{})𝒦_𝒪^{\{\mu _1\mu _2\mathrm{}\}}(P^{},P)U(P)`$
$`P^{}\left|\overline{\psi }_q\gamma ^{\mu _1}\psi _q\right|P`$ $`=`$ $`\overline{U}(P^{})\left\{\gamma ^{\mu _1}A_{q,1}(\mathrm{\Delta }^2)+\text{i}{\displaystyle \frac{\sigma ^{\mu _1\alpha }\mathrm{\Delta }_\alpha }{2m}}B_{q,1}(\mathrm{\Delta }^2)\right\}U(P)`$
$`P^{}\left|\overline{\psi }_q\gamma ^{\{\mu _1}\stackrel{}{D}^{\mu _2\}}\psi _q\right|P`$ $`=`$ $`\overline{U}(P^{})\{(\text{i})\overline{P}^{\{\mu _1}\gamma ^{\mu _2\}}A_{q,2}(\mathrm{\Delta }^2)`$
$`+(\text{i})\text{i}\overline{P}^{\{\mu _1}{\displaystyle \frac{\sigma ^{\mu _2\}\alpha }\mathrm{\Delta }_\alpha }{2m}}B_{q,2}(\mathrm{\Delta }^2)`$
$`+(\text{i}){\displaystyle \frac{\mathrm{\Delta }^{\{\mu _1}\mathrm{\Delta }^{\mu _2\}}}{m}}C_{q,2}(\mathrm{\Delta }^2)\}U(P)`$
$`P^{}\left|\overline{\psi }_q\gamma ^{\{\mu _1}\stackrel{}{D}^{\mu _2}\stackrel{}{D}^{\mu _3\}}\psi _q\right|P`$ $`=`$ $`\overline{U}(P^{})\{(\text{i})^2\overline{P}^{\{\mu _1}\overline{P}^{\mu _2}\gamma ^{\mu _3\}}A_{q,30}(\mathrm{\Delta }^2)`$ (28)
$`+(\text{i})^2\mathrm{\Delta }^{\{\mu _1}\mathrm{\Delta }^{\mu _2}\gamma ^{\mu _3\}}A_{q,32}(\mathrm{\Delta }^2)`$
$`+\text{i}(\text{i})^2\overline{P}^{\{\mu _1}\overline{P}^{\mu _2}{\displaystyle \frac{\sigma ^{\mu _3\}\alpha }\mathrm{\Delta }_\alpha }{2m}}B_{q,30}(\mathrm{\Delta }^2)`$
$`+\text{i}(\text{i})^2\mathrm{\Delta }^{\{\mu _1}\mathrm{\Delta }^{\mu _2}{\displaystyle \frac{\sigma ^{\mu _3\}\alpha }\mathrm{\Delta }_\alpha }{2m}}B_{q,32}(\mathrm{\Delta }^2)\}U(P).`$
In practice, we evaluate the ratio $`R_𝒪^{\text{cont,}\{\mu _1\mu _2\mathrm{}\}}(P^{},P)`$ by inserting the kernels $`𝒦_𝒪^{\{\mu _1\mu _2\mathrm{}\}}`$ from Eq. (28) and the projectors in Eq. (26,27) into Eq. (25). For given values of the nucleon momenta $`P`$ and $`P^{}`$, the ratio is computed numerically, determining the kinematical coefficients of the GFFs $`A(t)`$, $`B(t)`$ and $`C(t)`$. Eventually, we transform the free indices $`\{\mu _1\mu _2\mathrm{}\}`$ to Euclidean space with the use of Eq. (24). The system of linear equations is then constructed by equating the Euclidean lattice result and the transformed Minkowskian parameterization.

### A.1 Lattice operators

On the lattice side, we have H$`\left(4\right)`$-invariance instead of Lorentz-invariance, which complicates the classification and the renormalization of the operators. Following Göckeler et al. (1996b), we choose appropriate linear combinations of the operators $`𝒪_{\{\mu _1\mu _2\mathrm{}\}}`$, corresponding to representations of $`H(4)`$.

For the case of one derivative, $`n=2`$, there are three diagonal operators

$`𝒪_1^{\text{diag,n=2}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[𝒪_{11}+𝒪_{22}𝒪_{33}𝒪_{44}\right],`$
$`𝒪_2^{\text{diag,n=2}}`$ $`=`$ $`{\displaystyle \frac{1}{2^{1/2}}}\left[𝒪_{33}𝒪_{44}\right],`$
$`𝒪_3^{\text{diag,n=2}}`$ $`=`$ $`{\displaystyle \frac{1}{2^{1/2}}}\left[𝒪_{11}𝒪_{22}\right],`$ (29)
which are symmetric and traceless by construction. Additionally there are six non-diagonal combinations

$$𝒪_{\mu _1,\mu _2}^{\text{non-diag,n=2}}=\frac{1}{2}\left[𝒪_{\mu _1\mu _2}+𝒪_{\mu _2\mu _1}\right]=𝒪_{\left\{\mu _1\mu _2\right\}},\mu _{1,2}=1\mathrm{}4;\mu _1<\mu _2,$$
(30)
giving in the most general case a set of nine independent operators. There is no operator mixing for $`n=2`$.

In the case of two derivatives, $`n=3`$, we take

$`𝒪_1^{\text{n=3}}`$ $`=`$ $`\left({\displaystyle \frac{3}{2}}\right)^{1/2}\left[𝒪_{\left\{122\right\}}𝒪_{\left\{133\right\}}\right],𝒪_2^{\text{n=3}}={\displaystyle \frac{1}{2^{1/2}}}\left[𝒪_{\left\{122\right\}}+𝒪_{\left\{133\right\}}2𝒪_{\left\{144\right\}}\right],`$
$`𝒪_3^{\text{n=3}}`$ $`=`$ $`\left({\displaystyle \frac{3}{2}}\right)^{1/2}\left[𝒪_{\left\{211\right\}}𝒪_{\left\{233\right\}}\right],𝒪_4^{\text{n=3}}={\displaystyle \frac{1}{2^{1/2}}}\left[𝒪_{\left\{211\right\}}+𝒪_{\left\{233\right\}}2𝒪_{\left\{244\right\}}\right],`$
$`𝒪_5^{\text{n=3}}`$ $`=`$ $`\left({\displaystyle \frac{3}{2}}\right)^{1/2}\left[𝒪_{\left\{311\right\}}𝒪_{\left\{322\right\}}\right],𝒪_6^{\text{n=3}}={\displaystyle \frac{1}{2^{1/2}}}\left[𝒪_{\left\{311\right\}}+O_{\left\{322\right\}}2𝒪_{\left\{344\right\}}\right],`$
$`𝒪_7^{\text{n=3}}`$ $`=`$ $`\left({\displaystyle \frac{3}{2}}\right)^{1/2}\left[𝒪_{\left\{411\right\}}𝒪_{\left\{422\right\}}\right],𝒪_8^{\text{n=3}}={\displaystyle \frac{1}{2^{1/2}}}\left[𝒪_{\left\{411\right\}}+𝒪_{\left\{422\right\}}2𝒪_{\left\{433\right\}}\right],`$ (31)
while the four non-diagonal combinations are denoted by

$$𝒪_{\mu _1,\mu _2,\mu _3}^{\text{non-diag,n=3}}=𝒪_{\left\{\mu _1\mu _2\mu _3\right\}},\mu _{1,2,3}=1\mathrm{}4;\mu _1<\mu _2<\mu _3,$$

giving all in all 12 traceless and symmetric linear combinations. It turns out that some of the operators in Eqns. (31) mix with different non-symmetric representations (cf. Ref. Göckeler et al. (1996b)) under renormalization. Fortunately, the mixing coefficient turns out to be negligible, at least for the specific case considered in Göckeler et al. (1996a). It has to be checked if this holds true for all potential mixing candidates in the $`n=3`$ case.

The choice of index combinations for the construction of the H(4) operators is not unique. However, in the continuum limit all options will lead to the same GFFs. The interesting question of determining an optimal set of operators that minimizes lattice artifacts is beyond the scope of this presentation.

### A.2 Lattice momenta

The possible values of the nucleon momenta $`P`$ and $`P^{}`$ and therefore the momentum transfer squared $`t`$ are severely restricted by the lattice momenta. The general three-momentum on a periodic lattice is given by

$$\stackrel{}{P}=\frac{2\pi }{aL_s}\stackrel{}{p}=\frac{2\pi }{aL_s}(p_1,p_2,p_3),p_i=L_s,L_s+1,\mathrm{},L_s1,L_s,$$

where $`a`$ is the lattice spacing in GeV<sup>-1</sup> (we set $`a=1`$ in all intermediate steps), and $`L_s`$ is the spatial lattice extension, which is in our case $`16`$. It is clear that the lattice spacing times the spatial extension determines the smallest possible non-zero momentum, while the lattice spacing alone fixes the largest available momentum. In practice, however, we are even more restricted to values $`p_iL_s`$, because large momenta lead to considerable noise. For the nucleon four momenta $`P`$ and $`P^{}`$, we use the continuum dispersion relation, $`P_0=P_4=\sqrt{m^2+\stackrel{}{P}^2}`$, where $`m`$ is the nucleon mass, determined in our lattice calculation. In order to construct the maximally determined set of linear equations, Eq. (II.2), we need to know all $`\stackrel{}{p},\stackrel{}{p}^{}`$ leading to the same momentum transfer squared $`t`$. So far we have propagators available for the two sink momenta $`\stackrel{}{p}^{}=(0,0,0)`$ and $`\stackrel{}{p}^{}=(1,0,0)`$, and the value for $`\stackrel{}{q}\stackrel{}{p}^{}\stackrel{}{p}`$ is restricted by $`\stackrel{}{q}^2<6`$. We plan to extend the allowed values in future investigations. The total number of different $`t`$ is then 16, including the forward case.

Note, that $`t`$ depends in a nontrivial manner on $`\stackrel{}{p}^{}`$ and $`\stackrel{}{p}`$ individually through the energies $`P_4,P_4^{}`$, and not only on the relative momentum $`\stackrel{}{q}`$. For 8 of the 16 possible values of $`t`$, we encountered highly visible fluctuations in the nucleon propagators with negative values for the two point functions at large $`\tau \tau _{\text{src}}`$. We excluded the corresponding momentum combinations from our analysis. A list of all remaining $`\stackrel{}{q},\stackrel{}{p}^{},\stackrel{}{p}`$ and $`t`$, is shown in table 2.