Nucleon generalized form factors from two-flavor lattice QCD

We determine the generalized form factors, which correspond to the second Mellin moment (i.e., the first $x$-moment) of the generalized parton distributions of the nucleon at leading twist. The results are obtained within lattice QCD calculations with $N_f=2$ non-perturbatively improved Wilson fermions, employing a range of quark masses down to an almost physical value with a pion mass of about 150 MeV. We also present results for the isovector quark angular momentum and for the first $x$-moment of the transverse quark spin density. We compare two different fit strategies and find that directly fitting the ground state matrix elements to the functional form expected from Lorentz invariance and parameterized in terms of form factors yields comparable, and usually more stable results than the traditional two step approach.


I. INTRODUCTION
The understanding of hadron structure has greatly evolved over the last decades. The collected knowledge is parameterized by a large number of functions. Generalized parton distributions (GPDs) are one set of such functions. They parameterize, e.g., the transverse coordinate distribution of partons in a fast moving hadron and contain information on how these distributions depend on the parton or hadron spin direction. Pinning down all these multi-variable functions experimentally is unrealistic at present. Therefore, lattice QCD has to substitute some of the missing experimental data. With this article we contribute to the effort of various lattice groups to provide some of these needed results [1][2][3][4][5][6][7][8][9][10][11][12][13].
From the experimental point of view, GPDs play a similarly important role for the description of exclusive hadronic reactions as parton distribution functions (PDFs) do for inclusive reactions. The most extensively studied channel is deeply-virtual Compton scattering (DVCS), i.e., Compton scattering with a highly virtual incoming photon and a correspondingly large, spacelike momentum transfer Q 2 = −q 2 . One advantage of DVCS is that the GPD matrix element interferes with the well-known Bethe-Heitler cross section for which the final state photon is emitted from the scattered lepton. Thus the measured cross sections provide not only information on the absolute value of the DVCS correlators but also on their signs. In all generality, including spin effects, the experimental analysis becomes somewhat involved, as is, e.g., illustrated by the publications [14,15] * gunnar.bali@ur.de † sara.collins@ur.de ‡ meinulf.goeckeler@ur.de § rudolf.roedl@ur.de ¶ andreas.schäfer@ur.de * * andre.sternbeck@uni-jena.de of the Hermes experiment. For a recent careful theoretical analysis and references to experimental work see Ref. [16].
The theoretical understanding of GPDs and their moments, the generalized form factors (GFFs), has already a long history and is presented in the seminal work of Refs. [17][18][19][20][21]. More recent reviews can be found in Refs. [22,23]. The interest in some of the nucleon GPDs (there exist in total eight) is increased by the fact that they provide information on the elusive orbital angular momentum of partons in the nucleon. However, the physical interpretation in this case is not straightforward, because there exist inequivalent definitions of orbital angular momentum [19,24]. For recent discussions of this topic see, e.g., Refs. [25][26][27] and the articles cited therein. In this article we will not review the many fascinating aspects of GPDs but concentrate on our lattice calculation of the nucleon GFFs using well established techniques for the calculation of Mellin moments of GPDs; see, e.g., Ref. [28].
We remark that recently new methods have been proposed to "directly" calculate the Bjorken-x dependence of parton distribution functions (PDFs), distribution amplitudes (DAs), transverse momentum dependent PDFs (TMDPDFs) and GPDs, see, e.g., Refs. [29][30][31][32][33][34]. In these approaches Euclidean correlation functions are computed and then matched within collinear factorization to lightcone distribution functions, employing continuum perturbative QCD. For the example of DAs [33], some of us are involved in calculations with these new techniques, using the "momentum smearing" technique [35] to enable large hadron momenta to be realized, and found results that are consistent with, but less accurate than those obtained from the lowest non-trivial moment. This may change as smaller lattice spacings and larger computers become available. Here we will only determine the first x-moment, i.e., the second Mellin moment, to constrain the nucleon GPDs. This paper is organized as follows. In Section II we shortly review definitions and the operator product expansion for Mellin moments of GPDs. The lattice QCD techniques used to extract GFFs are introduced in Section III followed by a discussion of the numerical methods. In Sections V and VI we present our results. Some preliminary results have been reported in Refs. [5,7,10]. Finally, we investigate the transverse spin density of the nucleon in Section VII.

II. BASIC PROPERTIES OF GPDS
The starting point is the off-forward nucleon matrix element of a bilocal operator with quark flavor q The Wilson line U in Eq. (2) connects −λn/2 and +λn/2 on the light-cone (n 2 = 0). Depending on the Dirac structure, indicated by the symbol Γ in Eqs. (1) and (2), one can parametrize the matrix element M in terms of GPDs. For leading twist these read (see, e.g., Refs. [28,36]): with σ µν = i [γ µ , γ ν ]/2 and the nucleon spinors U (p , σ ) and U (p, σ). The GPDs, e.g., H q and E q , and the corresponding tensor structures γ µ and iσ µν ∆ ν /(2m N ) are written as vectors, where we apply a standard scalar product to simplify the notation and introduce the kinematic variables For the anti-symmetrization of indices we use the notation [. . .], e.g., The GPDs are functions of the three variables (x, ξ, t), such that H q = H q (x, ξ, t) etc. We define where t is the total momentum transfer squared which is related to the virtuality Q 2 = −t. The longitudinal momentum fraction x varies between −1 and 1 and the skewness ξ between 0 and 1. Negative values of x correspond to plus or minus (depending on the GPD) times the corresponding antiquark GPD at −x. In this work we restrict ourselves to the isovector case and therefore we only consider the above eight quark GPDs. An analogous set of gluonic GPDs exists, which we will not address here. For a more detailed discussion we refer the reader to Refs. [22,23,28,[36][37][38].
In physical terms (for |x| > ξ) GPDs parameterize the probability amplitude for a hadron to stay intact if a parton is removed at the light-cone point −λ/2 and replaced by a parton with different momentum at light-cone time λ/2. In practice, it is of crucial importance to find effective parameterizations of GPDs with a minimum number of parameters which are then fitted to experimental data see, e.g., Ref. [16]. Lattice input in principle allows one to pin down the values of these parameters, however, at present the accuracy of such studies is for many GPDs not yet sufficient to make a decisive impact.
As time is analytically continued to imaginary time to enable the numerical evaluation on the lattice, the lightcone loses its meaning. The operator product expansion (OPE) relates, however, Mellin moments of GPDs to local matrix elements that are amenable to lattice calculation. For H q and E q , for instance, these x-moments read (see, e.g., Refs. [28,36] where the real functions A q (t), B q (t) and C q (t) in the ξ-expansion on the RHS are the GFFs. The case n = 1 corresponds to the electromagnetic form factors F q 1 (t) = A q 10 (t) and F q 2 (t) = B q 10 (t). For n = 2 and t = 0 we obtain the average quark momentum fraction A q 20 = x q + , where, for this example, we indicated q ± = q ±q. Below we will drop this distinction since in the case of the vector and tensor GPDs the even moments automatically give the q + combination and the odd moments q − , while for axial GPDs it is the opposite.
In principle one can determine Mellin moments of GPDs for any n on the lattice, in practice one is restricted to the lowest few n. The reason for this restriction is twofold. On the one hand the signal to noise ratio becomes worse for an increasing number of covariant derivatives. On the other hand as n increases, mixing with lower dimensional operators will take place, resulting in divergences that are powers of the inverse lattice spacing a −1 . In this study we focus on the case n = 2, where such mixing does not occur. Similarly to elastic form factors, the respective GFFs are extracted from lattice calculations of two-and three-point correlation functions where the currents are the local twist-2 operators Here S µν and A µν denote symmetrization (also subtracting traces and dividing by n! for n indices) and antisymmetrization operators, respectively, and is the symmetric covariant derivative.
In the continuum we can decompose the matrix elements with the nucleon four-momentum (p µ ) = (E N ( p ), p ). In Section III we will show how we extract the matrix elements from the temporal dependence of the three-point correlation functions. The desired GFFs are contained in the Dirac structures: Some aspects of GFFs have been more intensively discussed in the literature than others, in particular: • As has already been mentioned above, in the forward limit (t = 0), A q 20 equals the average quark momentum fraction. Similar limits exist for A q 20 and A q T 20 and the polarized and transversity PDFs, respectively.
• Furthermore, in this limit A q 20 and B q 20 add up to twice the total angular momentum of the quark q plus that of the antiquarkq in the nucleon (the Ji sum rule [19]) such that represents the quark contribution to the nucleon spin. Combining J q with the quark spin contribution 1 2 ∆Σ q , one can also obtain the quark orbital angular momentum L q = J q − 1 2 ∆Σ q . We remark that this decomposition is not unique [24].
• The five GFFs A 20 , B 20 , A T 20 , B T 20 and A T 20 parameterize, after Fourier-transformation to impact parameter space, the first x-moment of the transverse spin density of a quark in a fast-moving nucleon [39].

III. EXTRACTING GENERALIZED FORM FACTORS
On the lattice, the GFFs are extracted from combinations of hadronic two-and three-point correlation functions in Euclidean space-time. The two-point function reads where the nucleon destruction and creation interpolators N and N are appropriate combinations of u and d (anti)quark fields C is the charge conjugation matrix. The lattice threepoint function is expressed as In this work we only consider isovector currents O, therefore, all quark lines are connected. To improve the overlap of our interpolators in Eqs. (13) with the physical ground state we employ the combination of APE and Wuppertal (Gauss) smearing techniques described in Refs. [40][41][42]. This procedure reduces the impact of excited states substantially. For the computation of Eq. (14), we use the sequential propagator method [43] which implies fixing the sink time t . We use the projector and contract it with the open spin indices of Eq. (14) to realize different spin projections and positive parity. For ρ = 1, 2, 3 we obtain the difference of the spinpolarization with respect to the quantization axis ρ, while ρ = 4 corresponds to the unpolarized case. The positive parity projection is only correct for zero momentum, however, excited state contributions (including states of In the case of the vector operator we work with multiplets that transform according to two distinct irreducible representations of the hypercubic group H(4) labeled as v 2,a and v 2,b . These are combinations of the operators in Eq. (16a) given by and respectively. The renormalized operators read where we use µ = 2 GeV as the renormalization scale. Note that the renormalization factors depend on the multiplet, i.e., they slightly differ for v 2,a and v 2,b . Similarly, the axial operators are renormalized with factors Z r2,a MS and Z (18) and (19). The tensor operators are renormalized with Z h1a MS and Z h 1b MS . The operator multiplets used in this case are listed in Appendix A.
A detailed description of the renormalization procedure, that consists of first non-perturbatively matching from the lattice to the RI -MOM scheme [45,46] and then translating perturbatively to the MS scheme, may be found in Ref. [44]. To make the article self-contained we summarize the basic steps in Appendix B, where we also address the error propagation from the renormalization constants to the GFFs. The relevant renormalization factors are summarized in Table I. They result from a reanalysis of the data presented in Ref. [44] and correspond to the physical input r 0 = 0.5 fm [47] and r 0 Λ MS = 0.789 [48]. Table II lists the relative errors on the renormalized GFFs, associated with the uncertainties in the renormalization constants; these amount to about 2%.
In the following we demonstrate the extraction procedure for the vector GFFs. The axial and tensor GFFs are treated analogously. We start by expanding Eq. (14) in terms of energy eigenstates where the ground state amplitude reads The exponentials contain the energy of the nucleon as a function of the considered spatial momentum, the Euclidean operator insertion time τ , and the sink time t . Up to lattice artifacts, the matrix elements of an operator O µν V,q (z) can be decomposed according to the Euclidean versions of Eqs. (9a) and (10a). In doing so, it is necessary to distinguish TABLE III. Parameters of the N f = 2 lattice ensembles used in this study. Latin numerals in the first column serve as ensemble identifiers. After the number of configurations N conf we list in parentheses the number of independent (randomly chosen) source positions that we average over within each gauge configuration. Wherever this is indicated by parentheses after the sink-source separation t /a, a smaller number of sources was used for this value. For more information about our set-up we refer to Ref. [41]. between the two multiplets v 2,a and v 2,b (cf. Eqs. (17) and (18)). The decomposition can be written as Applying the projection operator P ρ (cf. Eq. (15)) to C 3pt yields with and / p := iE N ( p )γ 4 + p · γ . The Z factors in Eq. (23) depend on the overlap of our nucleon interpolation operators with the nucleon ground state. They vary with momentum and smearing and can be extracted from the two-point correlation function C 2pt .
The right hand side of Eq. (24) contains the desired GFFs. The prefactors can be computed by inserting the respective Euclidean γ-matrices. Here we restrict ourselves to the final momentum p = 0. Taking all available combinations of operators [cf. Eqs. (17) and (18)], projections P ρ and momenta p for a fixed virtuality we obtain a linear system of equations with the GFF vector g V = (A 20 (t), B 20 (t), C 20 (t)) . The coefficient matrix M V consists of the prefactors calculated from Eq. (24) and F V is extracted from a fit of Eqs. (20) and (23) to lattice data for C 2pt and C 3pt . The number of columns of M V is equal to the number of unknown GFFs (in this case 3), but the number of rows depends on the available combinations. In almost all the cases this yields an overdetermined system of equations, meaning that the number of elements in F V , denoted with dim F V , is larger than the number of GFFs. Note that the individual rows of M V are either real or imaginary. 1 For a given ensemble this system of equations has to be solved separately for each virtuality to yield the GFFs as functions of t. In the general case we write Eq. (26) as FIG. 1. Overview of the nucleon energies for our ensembles. We compare the energies EN ( p ) and the errors extracted from a two-exponential fit shown as black error bars with the energies E c N expected from the continuum dispersion relation, which are depicted as colored boxes.
where Γ can take the values V , A, T and g q Γ is the vector of the respective GFFs (cf. Eqs. 10). Due to equivalent combinations of momenta and polarizations most rows in the matrix M Γ are equal or differ by a sign only. We average the corresponding correlation functions, which improves the signal-to-noise ratio considerably and reduces the number of equations.

A. Gauge ensembles
Our analysis is based on the large set of gauge configurations produced by the QCDSF and the RQCD (Regensburg QCD) collaborations using the standard Wilson gauge action with two mass-degenerate nonperturbatively improved clover fermions, see Table III. We have three different lattice spacings 0.081 fm, 0.071 fm and 0.060 fm. Despite the O(a) improved action, we expect discretization effects linear in the lattice spacing for our matrix elements since the currents are not improved. The pion masses range from about 490 MeV down to 150 MeV. In terms of Lm π we cover values from about 3.4 up to 6.7.

B. Fitting two-point correlation functions
We parameterize our two-point correlation functions with a two-exponential fit ansatz in order to create bootstrap ensembles for the fit parameters A( p ), E N ( p ), X( p ) and Y ( p ). To improve the signal, we average over all momentum combinations which lead to the same p 2 . Subsequently, we use Eq. (28b) to fix the overlap factors Z( p ) and Z( p ) which are needed to factor out F Γ from the three-point correlation functions (cf. Eq. (23)). The fit parameters X( p ) and Y ( p ) are introduced in order to parameterize the contributions from excited states. The parameter E N ( p ) represents the nucleon energy (we do not assume a functional form for the energy). However, our analysis assumes continuum symmetries. Therefore we restrict our lattice calculations to momenta whose fitted values for E N ( p ) are consistent with the continuum dispersion relation (cf. Fig. 1) The statistical errors are estimated by virtue of 500 bootstrap ensembles. We carefully study the fit-range dependence of the fit parameters. Therefore we consider the start time slices t s /a ∈ {2, 3} and vary the final time slice t f /a. We find that the impact of t s /a on the values for the GFFs is rather mild and therefore we fix t s /a = 2 in the following. In Fig. 2 we demonstrate how we choose the final time slice t f /a. We also try single exponential fits and find that they give similar results if one adjusts the fit ranges appropriately. However, the resulting errors on A( p ) are larger. Hence we use the two-exponential fit ansatz for our final analysis.

C. Three-point correlation functions
For the lattice calculations of three-point functions we use the sequential source method where we set the outgoing nucleon momentum p = 0 for all our ensembles. We parameterize the data using Eqs. (20) and (23) with E N ( p ) = m N . The initial energy E N ( p ) is determined from the continuum dispersion relation (29). The momentum restriction, which we discussed in the previous section, translates to a range 0 ≤ Q 2 < 0.6 GeV 2 for the three-point functions. With Z( p ) and Z( p ) having been determined from the two-point correlation functions, the only free parameter left is F q Γ . To achieve ground state dominance, one has to make sure that aN T t τ 0 (cf. Eq. (20)). We consider τ ∈ [τ s , τ e ] where τ s is well above zero and τ e well below t . The sink times vary with the ensemble (see the last column of Table III). In Section IV E we examine possible excited state contaminations.

D. Determination of the GFFs
As explained above, for every current Γ = V , A or T , quark flavor q and virtuality −t, we need to solve the linear system Eq. (27), i.e., F = M · g, to extract the relevant form factors g from the vector of inequivalent matrix elements F that correspond to non-vanishing rows of M . Here we drop all indices like the quark flavor q and Γ for convenience. In what follows m denotes the number of independent form factors while n ≥ m is the length of F. Consequently, M is a n × m matrix of maximal rank, i.e., rank(M ) = m.
The determination of the form factors is carried out in two ways. The first method consists of two steps: First we extract the ground state nucleon matrix elements F j from the lattice three-point function data c τ j , restricted to the range of insertion times τ ∈ [τ s , τ e ], through the numerical minimization of the χ 2 -function where δc τ j is the difference between the lattice data and the three-point function parametrization Eq. (23). The inverse covariance matrix cov −1 j depends on the insertion times τ and τ . One can easily generalize the fit to the situation of multiple source-sink distances t if this is required or include excited state contributions. The index j ∈ {1, . . . , n} runs over all possible polarizations ρ and initial momenta p (keeping the virtuality Q 2 fixed), which give non-vanishing contributions.
Once the fit parameters F j are determined, one can minimize (32) to determine the form factors g. The total number of parameters for this method is m + n and, in particular for large virtualities, this number can be quite large (up to 50). This is not the only problem but it can happen that the resulting value is quite large and it is not clear how one should deal with such a situation.
Ideally, should be zero but this is only possible if F is in the image of M (cf. Eq. (32)). Motivated by this observation, we carry out our fits employing a single step method, which combines the two subsequent steps into a single minimization problem, restricting the number of fit parameters to the relevant degrees of freedom. We start from the singular value decomposition with orthogonal matrices U ∈ R n×n , V ∈ R m×m and the matrix Σ ∈ R n×m , which has nonvanishing entries only on the diagonal. The pseudo-inverse Σ + is a m×n matrix that can easily be obtained, computing the inverses of the diagonal elements of Σ. Each vector F within the image of M can be uniquely expressed as a linear combination of the first m column vectors of U . Note that m = rank(M ). Substituting F → F( α ) in Eq. (31) (and thereby Eq. (30)), we obtain a modified χ 2 -function that depends on the parameters α i , where i ∈ {1, . . . , m}. Finally, we convert the extracted vector α to the desired GFF vector where in the last step Σ + is truncated to a m × m square matrix. In Fig. 3 we show for one example on the nearly physical quark mass ensemble VIII that this method works very well. In this case eight different lattice channels are well described in terms of three fit parameters. A comparison of the two fit methods shows that the results are consistent within errors for all GFFs and for all ensembles. The single step method, however, results in somewhat smaller statistical errors and a smoother Q 2 dependence, especially for the induced GFFs. In Fig. 4 we directly compare the two methods. For the final results we only use the single step method. In Fig. 5 we show all χ 2 dof values of all fits used in this paper to extract all considered GFFs: The correlated single step fits provide a very satisfactory description of the data.

E. Excited states
For some of our ensembles we have three-point function data for different source-sink separations. This allows us to analyze the influence of excited states on the GFFs. Our analysis is based on ensemble IV with five sourcesink separations in the range t /a ∈ [7, 17] and on ensemble VIII with three source-sink separations in the range t /a ∈ [9,15]. In physical units t = 15a corresponds to about 1 fm. Ensemble VIII has data for eight values of Q 2 and this ensemble corresponds to an almost physical pion mass. We show results only for this ensemble, but our findings are consistent for both ensembles.
For the tensor and axial GFFs we find that within statistical errors the Q 2 dependence is not affected by a variation of t . Only in the vector case, especially for A u−d 20 , excited state contaminations are visible (see Fig. 6). We have tried to parameterize these excitedstate contributions to the three-point function with various multi-exponential fit ansätze. This, however, introduces additional fit parameters, in particular the mass and the energy of the first excited state. The first excitation in the three-point function can be a multi-hadron state and hence its energy will in general not be well approximated by the single particle continuum dispersion relation. To parameterize excited state contributions clearly several source-sink separations are required. However, within present statistical errors little movement is visible for t 0.9 fm, even in the A 20 channel where we achieve the highest accuracy, see Fig. 6 for an example. We therefore have restricted our GFF fits to ranges of τ where the data are well described by a single exponential (cf. Fig. 3). In all the cases t is larger than 1 fm.

V. NUCLEON GFFS
Below we show results for the nucleon GFFs on a subset of the ensembles listed in Table III. We restrict ourselves to m π < 300 MeV and m π L > 3.4 and analyze the quark mass, volume and lattice spacing dependence. All results refer to the MS scheme at µ = 2 GeV.

A. Vector and axial GFFs
Results for the vector GFFs, and C u−d 2 , are shown in Fig. 7 (left) as a function of Q 2 = −t. We see that the discretization effects are negligible within errors (comparing ensembles I and XI, which give about the same pion mass and a similar value for Lm π ). Also the volume dependence (cf. V and VI) is small, although there is a slight trend towards larger values for B u−d is zero within errors. This agrees with the leading t-dependence expected from covariant baryon chiral perturbation theory (BChPT, see below). For large Q 2 we expect that A u−d 20 exhibits a dipole-like Q 2 -dependence, which we saw in our former study (cf. Fig. 2 of Ref. [5]).
Results for the axial GFFs are shown in the right panel of Fig. 7 We see that a change of volume, quark mass or lattice spacing has almost no effect on the data. Within errors these effects cannot be resolved. Both form factors grow approximately linearly for Q 2 → 0. For B u−d 20 the statistical errors become larger for Q 2 → 0 whereas the errors for A u−d 20 are nearly independent of Q 2 .

B. Tensor GFFs
Continuing with the tensor GFFs, we show results for Fig. 8 , are smaller in comparison and, besides a few outliers, are best described by a constant. However, a final conclusion cannot be drawn as the statistical errors for both GFFs are rather large. We also study the linear combination which corresponds to the combination of GPDs E T +2H T that is related to the Boer-Mulders function h ⊥ 1 [49]. We find that the statistical error of B q T 20 is significantly smaller compared to the individual errors of B q T 20 and A q T 20 (see Fig. 9). We will take advantage of this observation when looking at the transverse spin of the nucleon in Section VII. The results for B u−d T 20 are shown with the tensor GFFs in Fig. 8 for the same ensembles. The anticorrelations we find for B u−d T 20 are present for all ensembles.

VI. EXTRACTION OF J u−d
The GFFs A u−d 20 (t) and B u−d 20 (t) are of particular interest since for t → 0 they are related to the total angular momentum [19] In order to estimate J u−d at the physical pion mass we analyze our data for A u−d 20 (t) and B u−d 20 (t), employing the BChPT formulas of Ref. [50], which, however, we truncate at order m 3 π , The fit parameters T A 1 and T B 1 are added since our data extend up to virtualities −t ≈ (770 MeV) 2 m 2 π , however, these terms would naturally appear at the next order of BChPT. We determine the parameters by carrying out combined fits to our data sets for A u−d 20 (t, m π ) and B u−d 20 (t, m π ). The remaining parameters in Eqs. (38) and (39) are constrained to g A = 1.256, f π = 92.4 MeV and µ = 1.0 GeV.
Since it is not clear up to what values of −t and m π BChPT is applicable, we perform fits to all ensembles (set A) as well as fits using only ensembles with m π ≤ 300 MeV (set B). The result is shown in Fig. 10.     increases with m π → m phy π . For both sets we obtain values for χ 2 dof of about 0.75, hence we cannot use the χ 2 dof value to discriminate between the fit ranges. Instead, one may interpret the difference between fits A and B as a systematic uncertainty of the parameters. We also study the effect of the uncertainties of the renormalization constants using the strategy described in Appendix B 2. The final results are collected in Table IV. Within the errors, our values agree with the isovector results of Ref. [51].

VII. NUCLEON TOMOGRAPHY
We use our lattice results for the vector GFFs A 20 (t), B 20 (t) and the linear combination B T 20 (t) (cf. Eq. (36)) to investigate the transverse spin density of the nucleon. To this end, we transform these GFFs to the impact parameter space where we use the p-pole ansatz [52,53] for the interpolation of our lattice results. The impact parameter b ⊥ is defined in the transverse x-y plane. It measures the transverse distance from the "center of momentum" where x i is the momentum fraction of the ith parton [52,54]. We define To compute the transverse spin density, we also have to evaluate the derivative of G(b 2 ⊥ ) with respect to b 2 ⊥ : The Fourier transform (40)  K ν [52]: The transverse spin density ρ q (x, b ⊥ , s ⊥ , S ⊥ ) describes the probability to find a quark with longitudinal momentum fraction x, flavor q and transverse spin s ⊥ at a distance b ⊥ from the center of momentum of the nucleon with transverse spin S ⊥ . The explicit definition in terms of GPDs is given in Eq. (8) of Ref. [52]. Here we consider the two transverse spin combinations where the first line describes a transversely polarized quark in an unpolarized nucleon and the second an unpolarized quark in a transversely polarized nucleon. In terms of GFFs the first moment of ρ q (x, b ⊥ , s ⊥ , S ⊥ ) for these spin combinations reads For arbitrary spins S ⊥ and s ⊥ Eq. (47) will contain additional terms and we refer the reader to Refs. [52,53]. We fit the GFFs for ensemble VI to the p-pole ansatz Eq. (41). Due to the limited number of data points at our disposal, where we restricted ourselves to the kinematic range −t ≤ 0.6 GeV 2 , we find it impossible to simultaneously determine all three fit parameters, p, m p and G 0 . In particular the exponent p is strongly correlated with the pole mass m p . This is demonstrated in Fig. 11: An increase of p results in a larger value of m p , whereas χ 2 dof does not significantly change. Therefore, we cannot constrain p.
This arbitrariness means it is difficult to obtain reliable, parametrization independent results for the moment ρ q (b ⊥ , s ⊥ , S ⊥ ) as a function of b ⊥ . This distribution has been studied in the past (see, e.g., [53]), but we find that its shape strongly depends on the value of p. In Fig. 12 we show ρ q (b ⊥ , s ⊥ , S ⊥ ) for s ⊥ = (1, 0) and S ⊥ = (0, 0) for four distinct values if p ranging from 1.45 up to 3.0. We see that with increasing p the density becomes less localized in the impact parameter plane and the maximum of the density is shifted away from the center. This also holds for other spin combinations.
We found, however, meaningful integrated quantities which have a much milder p-dependence, namely the half b ⊥ -integrated moments with normalization factor The integrated moment ρ q + (s ⊥ , S ⊥ ) is the probability, weighted with the longitudinal momentum fraction x, to find a quark with flavor q in the upper part (b y ≥ 0) of the impact parameter space and ρ q − (s ⊥ , S ⊥ ) is the x-weighted probability to find a quark with flavor q in the lower part (b y ≤ 0). These integrated moments are a measure for the asymmetry of the transverse spin density. They depend much less on the value of p than ρ q (b ⊥ , s ⊥ , S ⊥ ) does. This is demonstrated in Fig. 13, where ρ d + and ρ d − are shown as functions of p for the transverse spin combination in Eq. (46a). Both integrated moments change only by 5% and 15%, respectively, while p doubles. We find this mild p-dependence for all considered transverse spin and flavor combinations and consider these integrated moments as the better candidates for reliable lattice estimates. Our results for ρ q ± for up and down quark for our two transverse spin combinations (Eq. (46)) are shown in Fig. 14. This figure corresponds to the value p = 2 but other values give very similar results.
We see the probability of a transversely polarized uor d-quark in an unpolarized nucleon is higher (∼ 70%) in the b y > 0 part of the impact parameter space than in the b y < 0 part (∼ 30%). For a transversely polarized nucleon however the probabilities of an unpolarized u-or d-quark differ: The unpolarized d-quark is more likely in the b y < 0 part (67%), while a u-quark is more likely in the b y > 0 part (60%) of the impact parameter space.

VIII. SUMMARY
We have calculated all quark GFFs, corresponding to operators with one derivative, of the nucleon GPDs at leading twist-2. Our lattice calculation includes the dominating connected contributions and neglects contributions from disconnected diagrams. The available gauge ensembles cover a wide range of quark masses and volumes. However, the three available lattice spacings only vary from 0.081 fm down to 0.060 fm. Within errors, all GFFs show a mild dependence on the quark mass, lattice spacing and volume.
We have compared two different fitting strategies for the GFFs and found that the direct fit method appears to be more reliable. With this method the number of fit parameters is reduced to the relevant degrees of freedom.
We recommend to use this method in future studies. The final results for the GFFs are shown in Figs. 7 and 8.
We have also studied the total angular momentum and the transverse spin density of quarks in the nucleon. Both quantities can be extracted from fits to our GFF data. For the total angular momentum we obtain a similar estimate in the isovector case as ETMC in Ref. [51]. Contributions from disconnected diagrams are not included in our lattice calculation. From Ref. [51] we know that these are small. Nevertheless, in the isoscalar case they should definitely be taken into account. For the second moment of the transverse spin density we have found that its distribution in impact parameter space strongly depends on the t-dependence of the GFF data. The shape of the distribution depends on the value of p that is used within a p-pole ansatz. High precision data at small and large values of −t would be required to eliminate this ambiguity. For integrated moments this situation improves. In Fig. 14 we provide lattice estimates for the x-weighted probabilities of a transverselypolarized (unpolarized) light quark in the upper or lower part of the impact parameter space, within an unpolarized (transversely-polarized) nucleon. Contributions from higher moments are not yet available but constitute an interesting object for future study.

Propagation of renormalization constant errors
Our estimates for the renormalization factors carry an uncertainty which has to be propagated into the GFFs. We do this in a very naive but conservative way by carrying out the whole analysis both using the central values of the renormalization factors and adding the error of these factors to their central values. The difference between these two sets of results is then due to the uncertainty of the renormalization. This procedure is applied to all ensembles and to all the available virtualities Q 2 . We find that the relative error is almost independent of Q 2 and the considered ensemble. Hence, for each GFF we decided to take the largest value of this uncertainty as an estimator of the error. These relative uncertainties are shown in Table II.