Nucleon Generalized Parton Distributions from Full Lattice QCD
Ph. Hägler, W. Schroers, J. Bratt, R.G. Edwards, M. Engelhardt, G.T. Fleming, B. Musch, J.W. Negele, K. Orginos, A.V. Pochinsky, D.B. Renner, D.G. Richards
aa r X i v : . [ h e p - l a t ] M a y JLAB-THY-07-651, MIT-CTP-3841, TUM-T39-07-09
Nucleon Generalized Parton Distributions from Full Lattice QCD
Ph. H¨agler
Institut f¨ur Theoretische Physik T39, Physik-Department der TU M¨unchen,James-Franck-Straße, D-85747 Garching, Germany
W. Schroers ∗ John von Neumann-Institut f¨ur Computing NIC/DESY, D-15738 Zeuthen, Germany
J. Bratt, J.W. Negele, and A.V. Pochinsky
Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139
R.G. Edwards and D.G. Richards
Thomas Jefferson National Accelerator Facility, Newport News, VA 23606
M. Engelhardt
Physics Department, New Mexico State University, Las Cruces, NM 88003-8001
G.T. Fleming
Sloane Physics Laboratory, Yale University, New Haven, CT 06520
B. Musch
Institut f¨ur Theoretische Physik T39, Physik-Department der TU M¨unchen,James-Franck-Straße, D-85747 Garching, Germany
K. Orginos
Department of Physics, College of William and Mary, Williamsburg VA 23187-8795and Thomas Jefferson National Accelerator Facility, Newport News, VA 23606
D.B. Renner † Department of Physics, University of Arizona, 1118 E 4th Street, Tucson, AZ 85721 (LHPC Collaboration) (Dated: October 26, 2018)We present a comprehensive study of the lowest moments of nucleon generalized parton distri-butions in N f = 2 + 1 lattice QCD using domain wall valence quarks and improved staggered seaquarks. Our investigation includes helicity dependent and independent generalized parton distribu-tions for pion masses as low as 350 MeV and volumes as large as (3.5 fm) , for a lattice spacing of0 .
124 fm. We use perturbative renormalization at one-loop level with an improvement based on thenon-perturbative renormalization factor for the axial vector current, and only connected diagramsare included in the isosinglet channel.
PACS numbers: 12.38.Gc,13.60.FzKeywords: Generalized parton distribution, lattice QCD, hadron structure
I. INTRODUCTION
Generalized parton distributions (GPDs) [1, 2, 3, 4] play a vital role in our understanding of the structure of thenucleon in terms of the fundamental building blocks of QCD, the quarks and gluons. Before the advent of GPDs,fundamental questions as to the origin of the spin of the nucleon, the decomposition of the nucleon total momentum,and the distribution and density of the nucleon constituents in position and momentum space seemed to be largelyunrelated. In some cases, it was even unclear how to formulate these questions in a theoretically sound way andhow to measure the underlying observables experimentally. With the introduction of GPDs, it is possible not only to ∗ Present address: Department of Physics, Center for Theoretical Sciences, National Taiwan University, Taipei 10617, Taiwan † Present address: DESY Zeuthen, Theory Group, Platanenallee 6, D-15738 Zeuthen, Germany x+ x x- x L , P ',' L P ),,( txH x FIG. 1: GPDs as part of a scattering amplitude. give precise definitions to quantities, such as the quark and gluon angular momentum contributions to the nucleonspin [5] and the probability densities of quarks and gluons in impact parameter space [6], but also to unify andextend the successful concepts of parton distribution functions (PDFs) and form factors. Nucleon generalized partondistributions are experimentally accessible in deeply virtual Compton scattering of virtual photons off a nucleon anda range of other related processes [7, 8, 9, 10] . Since these reactions involve in general convolutions of GPDs over thelongitudinal momentum fraction x , which makes it difficult if not impossible to map them over the whole parameterspace, the most stringent quantitative information on GPDs currently comes from quark PDFs and nucleon formfactors [11].Complementary to experimental efforts, lattice QCD offers a unique opportunity to calculate x -moments of GPDsfrom first principles. The first investigations of GPDs including studies of the quark angular momentum contributionsto the nucleon spin have been presented by the QCDSF collaboration in quenched QCD [12] and by LHPC/SESAMin N f = 2 lattice QCD [13]. Lattice results on nucleon GPDs published since then have provided important insightsinto the transverse structure of unpolarized nucleons [14], the lowest moments of polarized [15] and tensor GPDs [16],and transverse spin densities of quarks in the nucleon [17, 18]. With the exception of several initial studies [19, 20], allpreviously published lattice results on GPDs have been obtained from calculations in a two-flavor ”heavy pion world”with pion masses in the range of 550 to over 1000 MeV. In this work, we improve on previous studies by presenting acomprehensive analysis of the lowest three moments of unpolarized and polarized GPDs in N f = 2 + 1 lattice QCDwith pion masses as low as 350 MeV and volumes as large as (3 . .The paper is organized as follows. We begin with an introduction to the calculation of moments of GPDs in latticeQCD in section II. Section III describes the hybrid approach of using domain wall valence quarks with 2+1 flavorsof improved staggered sea quarks. In section IV we present our numerical results for the generalized form factors,including a discussion and interpretation of the quark orbital angular momentum and the transverse nucleon structure.Chiral extrapolations of selected lattice results to the physical pion mass are presented in section V. Conclusions aregiven in the final section. II. LATTICE CALCULATION OF MOMENTS OF GENERALIZED PARTON DISTRIBUTIONS
Generalized parton distributions determine off-forward matrix elements of gauge invariant light cone operators O Γ ( x ) = Z dλ π e iλx q ( − λn P e − ig R λ/ − λ/ dα n · A ( αn ) q ( λn , (1)where x is the momentum fraction, n is a light cone vector and Γ = /n or Γ = /nγ . The twist-2 tensor GPDs [21]related to Γ = n µ σ µj , j = 1 , H , E , e H and e E are defined by the parametrizations h P ′ , Λ ′ | O n ( x ) | P, Λ i = hh /n ii H ( x, ξ, t ) + n µ ∆ ν m hh iσ µν ii E ( x, ξ, t ) , (2)and h P ′ , Λ ′ | O nγ ( x ) | P, Λ i = hh /nγ ii e H ( x, ξ, t ) + n · ∆2 m hh γ ii e E ( x, ξ, t ) , (3)where we use the short-hand notation hh Γ ii = U ( P ′ , Λ ′ )Γ U ( P, Λ) for products of Dirac spinors U , and where ∆ = P ′ − P , t = ∆ and ξ = − n · ∆ /
2. In Eqs. (2) and (3) we suppress the dependence of the GPDs on the resolutionscale µ . An illustration of the GPDs parametrizing the lower part of the handbag diagram is given in Fig. (1). Themomentum fractions x and ξ both have support in the interval [ − , +1]. Depending on x , there are three kinematicregions, which offer different interpretations for the GPDs. For x ∈ [ ξ,
1] and x ∈ [ − , − ξ ], the GPDs describe theemission and reabsorption of quarks and anti-quarks, respectively. In the case that x lies in the interval [ − ξ, ξ ], theydescribe the emission of a qq -pair.In our lattice calculations, we do not work directly with the bi-local operators in Eq. (1) but instead considermoments, defined by the integral R − dx x n − f ( x ), of the operators in Eqs. (2,3), leading to towers of symmetrized,traceless local operators O { µ ...µ n } [ γ ] = q (0) γ { µ [ γ ] i ↔ D µ · · · i ↔ D µ n } q (0) , (4)where [ γ ] denotes the possible inclusion of the corresponding matrix, the curly brackets represent a symmetrizationover the indices µ i and subtraction of traces, and ↔ D = 1 / → D − ← D ). We relate nucleon matrix elements of the tower oflocal operators in Eq. (4) to x -moments of the twist-2 GPDs. To this end, we parametrize off-forward matrix elements h P ′ , Λ ′ |O { µ ...µ n } | P, Λ i in terms of the generalized form factors A ni ( t ), B ni ( t ), C n ( t ), e A ni ( t ), and e B ni ( t ). Apart frompotential difficulties related to lattice operator mixing for higher moments n , lattice measurements of the operatorsin Eq. (4) become increasingly noisy as the number of derivatives increases, and we therefore restrict our calculationsto n ≤
3. The parametrization of nucleon matrix elements of Eq. (4) in terms of generalized form factors (GFFs) for n = 1, 2 and 3 reads [4, 22] h P ′ |O µ | P i = hh γ µ ii A ( t ) + i m hh σ µ α ii ∆ α B ( t ) , h P ′ |O { µ µ } | P i = ¯ P { µ hh γ µ } ii A ( t ) + i m ¯ P { µ hh σ µ } α ii ∆ α B ( t ) + 1 m ∆ { µ ∆ µ } C ( t ) , h P ′ |O { µ µ µ } | P i = ¯ P { µ ¯ P µ hh γ µ } ii A ( t ) + i m ¯ P { µ ¯ P µ hh σ µ } α ii ∆ α B ( t )+ ∆ { µ ∆ µ hh γ µ } ii A ( t ) + i m ∆ { µ ∆ µ hh σ µ } α ii ∆ α B ( t ) , (5)for the unpolarized case, and h P ′ |O µ γ | P i = hh γ µ γ ii e A ( t ) + 12 m ∆ µ hh γ ii e B ( t ) , h P ′ |O { µ µ } γ | P i = ¯ P { µ hh γ µ } γ ii e A ( t ) + 12 m ∆ { µ ¯ P µ } hh γ ii e B ( t ) , h P ′ |O { µ µ µ } γ | P i = ¯ P { µ ¯ P µ hh γ µ } γ ii e A ( t ) + 12 m ∆ { µ ¯ P µ ¯ P µ } hh γ ii e B ( t )+ ∆ { µ ∆ µ hh γ µ } γ ii e A ( t ) + 12 m ∆ { µ ∆ µ ∆ µ } hh γ ii e B ( t ) , (6)for the polarized case. Here and in the following we set ¯ P = ( P ′ + P ) / H n ( ξ, t ) ≡ Z − dx x n − H ( x, ξ, t ) , E n ( ξ, t ) ≡ Z − dx x n − E ( x, ξ, t ) , e H n ( ξ, t ) ≡ Z − dx x n − e H ( x, ξ, t ) , e E n ( ξ, t ) ≡ Z − dx x n − e E ( x, ξ, t ) , (7)are given by polynomials in the longitudinal momentum transfer ξ and the GFFs. For the lowest three moments, thecorresponding relations read H n =1 ( ξ, t ) = A ( t ) , H n =2 ( ξ, t ) = A ( t ) + (2 ξ ) C ( t ) , H n =3 ( ξ, t ) = A ( t ) + (2 ξ ) A ( t ) ,E n =1 ( ξ, t ) = B ( t ) , E n =2 ( ξ, t ) = B ( t ) − (2 ξ ) C ( t ) , E n =3 ( ξ, t ) = B ( t ) + (2 ξ ) B ( t ) , e H n =1 ( ξ, t ) = e A ( t ) , e H n =2 ( ξ, t ) = e A ( t ) , e H n =3 ( ξ, t ) = e A ( t ) + (2 ξ ) e A ( t ) , e E n =1 ( ξ, t ) = e B ( t ) , e E n =2 ( ξ, t ) = e B ( t ) , e E n =3 ( ξ, t ) = e B ( t ) + (2 ξ ) e B ( t ) . (8)The aim of our calculation is to extract the GFFs as functions of the momentum transfer squared, t , from nucleontwo- and three-point-functions as described below. Once the GFFs have been obtained, the complete ξ -dependence =t P ),( Dt (cid:1) snk ,' t P ),( Dt (cid:0) snk ,' t P =t P a) b) FIG. 2: Connected (a) and disconnected (b) diagrams in unquenched lattice QCD with an operator insertion at τ and finitemomentum transfer ∆. of the moments of the GPDs is directly given by Eqs. (8). Let us note that the Mellin-moments in Eq. (7) are takenwith respect to the entire interval from x = − ξ = 0 correspond to sums and differences of contributions from quarks q and anti-quarks q .For example, H nq ( ξ = 0 , t ) = Z dx x n − (cid:0) H q ( x, , t ) + ( − n H q ( x, , t ) (cid:1) ,E nq (0 , t ) = Z dx x n − (cid:0) E q ( x, , t ) + ( − n E q ( x, , t ) (cid:1) , e H nq (0 , t ) = Z dx x n − (cid:0) e H q ( x, , t ) + ( − ( n − e H q ( x, , t ) (cid:1) , e E nq (0 , t ) = Z dx x n − (cid:0) e E q ( x, , t ) + ( − ( n − e E q ( x, , t ) (cid:1) . (9)Such a simple decomposition is not possible for non-zero longitudinal momentum transfer ξ = 0. We denote the forwardlimit values of the moments of H and e H in Eq. (9) by h x n − i q = H nq (0 ,
0) = A qn (0) and h x n − i ∆ q = e H nq (0 ,
0) = e A qn (0),where h x n − i q and h x n − i ∆ q correspond to the moments of unpolarized and polarized quark parton distributions.Below, we give a brief summary of the methods and techniques used to extract moments of generalized partondistributions in lattice QCD. For details, we refer the reader to [13, 23]. As usual, the matrix elements, Eqs. (5,6),are calculated from the ratio of nucleon three-point and two-point functions: C ( τ, P ) = X j,k (Γ unpol ) jk h N k ( τ, P ) N j ( τ src , P ) i ,C O ( τ, P ′ , P ) = X j,k (Γ pol ) jk h N k ( τ snk , P ′ ) O ( τ, ∆) N j ( τ src , P ) i , (10)where Γ unpol = (1 + γ ) / pol = (1 + γ )(1 + iγ γ ) /
2. The nucleon source, ¯ N ( τ, P ), and sink, N ( τ, P ), createand annihilate states with the quantum numbers of the nucleon. To maximize the overlap with the ground state, weused the smeared sources given in [23]. The three-point-function C O ( τ, P ′ , P ) with the operator insertion at τ isillustrated in Fig. (2) in terms of quark propagators, showing examples of connected and disconnected contributionsin an unquenched lattice calculation.Using the transfer matrix formalism, we can rewrite Eqs. (10) to obtain C ( τ, P ) = e − E ( P )( τ − τ src ) (cid:0) Z ( P ) Z ( P ) (cid:1) / E ( P ) + mE ( P ) + higher states , (11) C O ( τ, P ′ , P ) = e − E ( P )( τ − τ src ) − E ( P ′ )( τ snk − τ ) (cid:0) Z ( P ) Z ( P ′ ) (cid:1) / E ( P ′ ) E ( P ) Tr (cid:8) Γ pol ( i/P ′ − m ) (cid:0) aA ( t ) + bB ( t ) + · · · (cid:1) ( i/P − m ) (cid:9) + higher states , (12)where the factors a, b, . . . represent the prefactors (including Dirac-matrices) of the corresponding GFFs A ( t ) , B ( t ) , . . . in the parametrizations in Eqs. (5,6), transformed to Euclidean space. Higher states with energies E > E in Eqs. (11)and (12) are suppressed when τ snk − τ ≫ / ( E − E ) and τ − τ src ≫ / ( E − E ).In order to cancel the exponentials and Z -factors in Eq. (12) for zero and non-zero momentum transfer ∆, weconstruct the ratio of two- and three-point-functions R O ( τ, P ′ , P ) = C O ( τ, P ′ , P ) C ( τ snk , P ′ ) (cid:20) C ( τ snk − τ + τ src , P ) C ( τ, P ′ ) C ( τ snk , P ′ ) C ( τ snk − τ + τ src , P ′ ) C ( τ, P ) C ( τ snk , P ) (cid:21) / . (13)For an operator-insertion sufficiently far away from the source and the sink in the Euclidean time direction, higherstates are negligible, and the ratio R O ( τ, P ′ , P ) exhibits a plateau in τ . We finally average over the plateau regionfrom τ min to τ max to obtain an averaged ratio R O ( P ′ , P ). On a finite periodic lattice with spatial extent L s , three-momenta are given by ~P = 2 π/ ( aL s ) ~n with integer components n i = − L s / , . . . , L s /
2, and for the nucleon energy weuse the continuum dispersion relation P = p m + ~P . Therefore, the discrete lattice momenta result in a finite setof values for the momentum transfer squared t which can be realized in our calculation.In order to obtain symmetric and traceless operators O { µ µ ... } , we have to choose specific linear combinationsof the indices. For the diagonal operators typical examples are O n =2 i =1 = ( O + O − O − O ) / O n =3 i =2 =( O + O − O ) / √
2, where O n =21 belongs to the 3-dimensional irreducible H (4)-representation τ (3)1 for n = 2and O n =3 i =2 is a member of the 8-dimensional representation τ (8)1 for n = 3 [24]. The sets of operators O ni we are usingare the same as in [13]. Altogether, there are 9 linearly independent index combinations for n = 2 and 12 for n = 3. Inorder to be able to compare our results with experiment, the operators have to be renormalized and transformed to theMS-scheme at a renormalization scale µ . In general, operators mix under renormalization, and the renormalizationmatrix Z O is non-diagonal. We will denote the renormalized operators in the MS-scheme by O n, MS i = Z O ij O nj . Somedetails concerning the renormalization procedure and numerical results for the renormalization constants will bediscussed at the end of the next section.Based on the renormalized operators, we compute the averaged ratio R O ( P ′ , P ) and equate it with the continuumparametrization in terms of the GFFs given in Eq. (12). This is done simultaneously for all momentum combinations P and P ′ corresponding to the same momentum transfer squared t and all contributing symmetric and tracelessoperators O n, MS i , giving a finite set of linear equations R O ,k ( P ′ , P ) = c A ( t ) + c B ( t ) + . . . ,R O ,l ( P ′ , P ) = c A ( t ) + c B ( t ) + . . . ,R O ,m ( P ′ , P ) = c A ( t ) + c B ( t ) + . . . , · · · , (14)where ( P ′ j − P j ) = t for all j = 1 , , . . . . The coefficients c ij in Eqs. (14) only depend on the nucleon mass m andthe momenta P, P ′ and are calculated from the traces in Eq. (12). Finally, the set of equations (14), which in generalis overdetermined, is solved numerically to extract the GFFs. The statistical errors for the GFFs are obtained froma jackknife analysis. III. LATTICE SIMULATION USING DOMAIN WALL VALENCE QUARKS WITH STAGGERED SEAQUARKS
Since calculations at physical quark masses are prohibitively expensive with current algorithms and machines,we have used dynamical quark configurations at the lightest masses available, and where feasible, have used chiralperturbation theory to extrapolate to the physical mass. Staggered sea quarks with the Asqtad improved action werechosen because the computational economy of staggered quarks enabled the MILC collaboration to generate largesamples of configurations at low masses on large spatial volumes [25, 26], which they freely made available to thelattice community.Chiral symmetry is crucial for avoiding some operator mixing, convenient for operator renormalization, and valuablefor chiral extrapolation. Furthermore, the four tastes associated with staggered fermions immensely complicatecalculating operator matrix elements in nucleon states. Hence, we chose a hybrid action utilizing chirally symmetricvalence domain wall fermions (DWF) on an improved staggered fermion sea. Although this hybrid scheme breaksunitarity at finite lattice spacing, given the arguments that the valence and sea actions separately approach the physicalcontinuum limit[27], we expect that the hybrid action also approaches the physical continuum limit. Furthermore,partially quenched mixed action chiral perturbation theory calculations are now becoming available for quantitativecontrol of the continuum limit. We also note that hybrid actions have been successfully used in other contexts where,for example, the NRQCD action for valence quarks was combined with improved staggered sea quarks[28] and wassuccessful in predicting mass splitting in heavy quark systems.In our calculation, we used MILC configurations [29] both from the NERSC archive and provided directly by thecollaboration. We then applied HYP-smearing [30] and bisected the lattice in the time direction. We have chosengauge fields separated by 6 trajectories. Furthermore, we alternate between the first temporal half (time slices 0 to31) and the second temporal half (time slices 32 to 63) on successive gauge configurations. In these samples we didnot find residual autocorrelations. The scale is set by the lattice spacing a = 0 . L . They preserve the Ward-Takahashiidentity [35] even at finite lattice spacing in the limit L → ∞ . At finite values of L a residual explicit breakingof chiral symmetry is still present which can be parameterized by an additional mass term in the Ward-Takahashiidentity [36, 37]. In our calculations, we have kept this additional mass, ( am ) res , at least an order of magnitudesmaller than the quark mass, ( am ) DWF q [19]. To the extent that ( am ) res is negligible, perturbative renormalizationof O [ γ ] is independent of γ in the chiral limit and the non-perturbative renormalization of quark bilinear currentsyields the same renormalization coefficients for the axial and the vector currents in the chiral limit [38].We now consider the parameters entering the DWF action. The domain wall action realizes chiral symmetry byproducing right-handed states on one domain wall that decay exponentially away from the wall and exponentiallydecaying left handed states on the other wall. To the extent that no low eigenmodes associated with dislocations (orrough fields) destroy the exponential decay and that the fifth dimension L is large enough, chiral symmetry will benearly exact and ( am ) res will be small. HYP smearing was essential to reduce the effect of low eigenmodes, but itis still necessary to use spectral flow to determine a value of the domain wall mass, M , for which the density of loweigenmodes is as small as possible. This was done on an ensemble of test configurations, with the result that we use M = 1 .
7. As discussed below, L was then tuned to keep ( am ) res below 10% of the quark mass and to have negligibleeffect on our lattice observables. Finally, the quark mass was tuned to produce a pion mass equal to the Goldstonepion mass for the corresponding MILC configurations. A. Tuning the fifth dimension
The extent of the fifth dimension, L , has been adjusted such that the residual mass, ( am ) res is at least an orderof magnitude smaller than the quark mass itself. This tuning is most relevant at the lightest quark mass since inthat case the computational cost is largest and thus our tuning should be optimal. In addition, the breaking of chiralsymmetry is also expected to be the largest and the resulting L provides a minimum value needed for our calculationsat the higher masses. The residual explicit chiral symmetry breaking characterized by ( am ) res is obtained from [36]∆ µ A aµ = 2 m q J aq ( x ) + 2 J a q ( x ) , (15)where J a q ( x ) ≈ m res J a ( x ) , (16)which holds up to O ( a ).We have run simulations using two samples of 25 configurations with volume Ω = 20 ×
32: three degeneratedynamical Asqtad quarks with bare masses ( am ) Asqtad,sea q = 0 .
050 (denoted as “heavy” and corresponding to m π ∼ am ) Asqtad,sea q = 0 . , .
050 (termed “light” and corresponding to m π ∼
350 MeV). The corresponding bare DWF masses have been adjusted to ( am ) DWF q = 0 . . L = 16 just fulfills our requirement, while in the heavy quark case L = 16 more than satisfies it. This confirms ourexpectation that the value of L chosen at the lightest quark mass sets the lower limit for the other masses as well.One quantitative check that L = 16 is adequate is provided by the dependence of masses on L as shown inFigs. (4), (5), (6) and (7). The leading effect of m res is to shift the quark mass, so that when L is sufficiently largethat this is the only effect, m π ∝ ( m q + m res ). Figs. (4) and (5) show the difference in the ratio m π / ( m q + m res )at a general value of L and at L = 16, and indicate that the difference is essentially consistent with zero for L >
16. We expect the shifts in the nucleon mass induced by these small shifts in the pion mass to be negligible,and indeed, Figs. (6) and (7) show that the differences between the nucleon mass at a general value of L and at L = 16 are consistent with zero for L >
16. Hence, we choose L = 16 to be a good compromise between accuracyand performance. ( a m ) r e s heavylight FIG. 3: Residual quark mass as a function of L for the two samples (heavy and light) of 25 configurations each.FIG. 4: Dependence of the pion mass on the extentof the fifth dimension L for heavy quarks. FIG. 5: Dependence of the pion mass on the extentof the fifth dimension L for light quarks. B. Tuning the quark mass
We define the light quark masses in our hybrid theory by matching the pion mass in two calculations: (i) using twoplus one flavors of dynamical Asqtad sea fermions and Asqtad valence fermions[29] and (ii) using the pion mass inour hybrid calculation with Asqtad dynamical sea fermions and valence domain-wall fermions with L = 16. Becauseof the four tastes and correspondingly sixteen light pseudoscalar mesons in the staggered theory, it is necessaryto choose between matching the lightest pseudoscalar mass, corresponding to the Goldstone pion of the theory, or FIG. 6: Dependence of the nucleon mass on theextent of the fifth dimension L for heavy quarks. FIG. 7: Dependence of the nucleon mass on theextent of the fifth dimension L for light quarks. dataset Ω am ) Asqtad q ( am ) DWF q ( am ) Asqtad π ( am ) DWF π ( am ) Asqtad N ( am ) DWF N m DWF π [MeV]1 20 ×
32 425 0 . / .
050 0 . . . . . . . . / .
050 0 . . . . . . . . / .
050 0 . . . . . . . . / .
050 0 . . . . . . . . / .
050 0 . . . . . . . ×
32 270 0 . / .
050 0 . . . . . some appropriately defined average. In this work, we have chosen to match the Goldstone pion, and the results oftuning the domain wall quark mass such that the domain wall pion mass agrees within one percent with the AsqtadGoldstone pion mass are shown in Table I. The substantial difference between the bare quark masses for Asqtad andDWF valence quarks reflects the significant difference in renormalization for the two actions. An observable physicaldifference is the fact that once the DWF quark masses have been adjusted to fit the Asqtad Goldstone pion masses,the DWF nucleon masses are approximately 6% lower than the corresponding Asqtad nucleon masses. We attributethis to the range of pseudoscalar masses in the staggered theory and note that had we used a heavier quark massso that the DWF pion fit some appropriately weighted average of the staggered pion masses, then the DWF nucleonwould have been heavier. C. Operator renormalization
The quark bilinear operators in Eq. [4] are renormalized using a combination of one-loop perturbation theory andnon-perturbative renormalization of the axial vector current.Our lattice calculations using lattice regularization with cutoff 1 /a are related to physical observables at scale µ in the MS renormalization scheme in 1-loop perturbation theory by O MS i ( µ ) = X j (cid:18) δ ij + g π N c − N c (cid:16) γ MS ij log( µ a ) − ( B LAT Tij − B MS ij ) (cid:17)(cid:19) · O LAT Tj ( a )= Z O ij · O LAT Tj ( a ) , (17)where the anomalous dimensions γ ij and the finite constants B ij have been calculated for domain wall fermions withHYP smearing in Refs. [39, 40]. Because the renormalization factors for operators with and without γ are identicalat quark mass zero, we use mass independent renormalization with all renormalization constants defined at quarkmass zero. The renormalization factors Z O ij for domain wall mass M = 1 . g / (12 π ) = 1 / .
64 from Refs. [39, 40]. By virtue of the operator H (4) Z O , pert ¯ q [ γ ] γ { µ D ν } q τ (3)1 . q [ γ ] γ { µ D ν } q τ (6)1 . q [ γ ] γ { µ D ν D ρ } q τ (4)2 . q [ γ ] γ { µ D ν D ρ } q τ (8)1 . µ = 1 /a . suppression of loop integrals by HYP smearing, the ratio of the one-loop perturbative renormalization factor for ageneral bilinear operator to the renormalization factor for the axial current is within a few percent of unity, suggestingadequate convergence for this ratio at one-loop level. Since one element in the calculation common to all operatorsarising from the wave function renormalization in the fifth dimension is not small, it is desirable to determine thisone common factor non-perturbatively. This is accomplished using the fact that the renormalization factor, Z A , forthe four dimensional axial current operator A µ = ¯ qγ µ γ q may be calculated using the five dimensional conservedaxial current for domain wall fermions A µ by the relation[36] hA µ ( t ) A µ (0) i = Z A h A µ ( t ) A µ (0) i . Hence the completerenormalization factor is written as the exact axial current renormalization factor times the ratio of the perturbativerenormalization factor for the desired operator divided by the perturbative renormalization factor for the axial current. (cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:8)(cid:9)(cid:10)(cid:11) (cid:12)(cid:13)(cid:14)(cid:15) (cid:16)(cid:17)(cid:18)(cid:19)(cid:20)(cid:21)(cid:22)(cid:23)(cid:24)(cid:25)(cid:26)(cid:27) (cid:28)(cid:29)(cid:30)(cid:31) !" FIG. 8: Unpolarized (vector) generalized n = 2 form factors for the flavor combinations u − d (left) and u + d (right).Disconnected contributions are not included. That is, Z O = Z O , pert Z pert A · Z nonpert A . (18)In the continuum, because of Lorentz invariance, the totally symmetric operator ¯ q [ γ ] γ { µ D ν D ρ } q cannot mix withthe mixed symmetry operator ¯ q [ γ ] γ [ µ D { ν ] D ρ } q , where the square brackets denote antisymmetrization. In contrast,on the lattice, both operators appear in the same representation, τ (8)1 , so that they can and do mix. However, themixing coefficient[39, 40], Z O ij = 2 . × − , is very small, so that we have ignored the contribution of the mixedsymmetry operator in this present work.All results below have been transformed to a scale of µ = 4 GeV .0 m œß , u (cid:252) dA (cid:253) B (cid:254) m (cid:255)(cid:1) , u (cid:4) dA (cid:3) B (cid:0) m (cid:2)(cid:5) , u (cid:6) dA (cid:7) B (cid:8) m (cid:9)(cid:10) , u (cid:11) dA (cid:12) B (cid:13) m (cid:14)(cid:15) , u (cid:16) dA (cid:17) B (cid:18) m (cid:19)(cid:20) , u (cid:21) dA (cid:22) B (cid:23) m (cid:24)(cid:25) , u (cid:26) dA (cid:27) B (cid:28) m (cid:29)(cid:30) , u (cid:31) dA B ! m " , u $ dA % B & m ’( , u ) dA * B + m ,- , u . dA / B m , u dA B t GeV t : GeV ; < t = GeV > ? t @ GeV A A B , B C A D , B E A F , B G A H , B I A J , B K A L , B M A N , B O A P , B Q A R , B S A T , B U A V , B W A X , B Y FIG. 9: Polarized (axial vector) generalized n = 2 form factors for the flavor combinations u − d (left) and u + d (right).Disconnected contributions are not included. IV. NUMERICAL RESULTS FOR THE GENERALIZED FORM FACTORS
Since two point functions taken at the sink τ = τ snk , C ( τ snk , P ′ ) and C ( τ snk , P ) in the ratio Eq. (13) decayexponentially for the full Euclidean distance τ snk − τ src , they are particularly subject to statistical noise. In theworst case, they may become negative, which we observe for three values of the momentum transfer for the dataset m = 0 . , . The corresponding datapoints are excluded from our analysis. Our numerical results for the completeset of unpolarized and polarized n = 1 , , A A A m Z[ , u \ d A A A m ]^ , u _ dA A m ‘a , u b d A A m cd , u e dA A A m fg , u h d A A A m ij , u k dA A A m lm , u n d A A A m op , u q dA A A m rs , u t d A A A m uv , u w dA A A m xy , u z d A A A m {| , u } d A , A , A A , A , A A , A , A A , A , A A , A , A A , A , A A , A , A A , A , A A , A , A A , A , A A , A , A A , A , A ~ t (cid:127) GeV (cid:128) (cid:129) t (cid:130) GeV (cid:131) (cid:132) t (cid:133) GeV (cid:134) (cid:135) t (cid:136) GeV (cid:137) FIG. 10: Flattening of the slope of the A n GFFs with increasing n for flavor combinations u − d (left) and u + d (right). Thesolid curves and error bands correspond to dipole fits described in the text. Disconnected contributions are not included. In Figs. (8) and (9) we show results for the vector generalized form factors A , B , C and axial vector GFFs e A , e B as functions of the momentum transfer squared t . We observe that the absolute values in the isovector andisosinglet channels are in qualitative agreement with the predictions from large N c counting rules, see e.g. [41], forthe unpolarized GFFs | A u + d | ∼ N c ≫ | A u − d | ∼ N c , | B u − d | ∼ N c ≫ | B u + d | ∼ N c , | C u + d | ∼ N c ≫ | C u − d | ∼ N c . (19)In the polarized case, the inequalities from the counting rules are not satisfied nearly as strongly. Whereas the countingrules predict: | e A u − d | ∼ N c ≫ | e A u + d | ∼ N c , | e B u − d | ∼ N c ≫ | e B u + d | ∼ N c , (20)2 Π @ GeV D < r ¦ > (cid:144) @ f m D n = = = FIG. 11: Two dimensional rms radii of the vectorGFFs versus m π for the flavor combination u − d . The results for m π = 354 MeV, L = 20 aredisplayed in gray. Π @ GeV D < r ¦ > (cid:144) @ f m D n = = = FIG. 12: Two dimensional rms radii of the axialvector GFFs versus m π for the flavor combination u − d . The results for m π = 354 MeV, L = 20 aredisplayed in gray and slightly shifted to the rightfor clarity. our results for e A u − d are only slightly larger than e A u + d , and although the errors are large, e B u − d appears to becomparable to e B u + d rather than dominating it. Finally, our results disagree with the predicted hierarchy betweendifferent types of GFFs: | B u − d | ∼ N c ≫ | A u + d | ∼ N c , (21)since the lattice results (at non-zero t ) clearly give A u + d > B u − d . It would be valuable to understand why thesecounting rules are only partially satisfied.For future reference, it is important to note that the GFF C , which gives rise to the ξ -dependence of the n = 2moment of the GPDs H ( x, ξ, t ) and E ( x, ξ, t ), is compatible with zero for u − d , over the full range of momentumtransfer squared t ≈ − . . . . − . . Similarly, the isosinglet GFF B u + d , which is one of the terms in thecontribution of the total angular momentum to the nucleon spin, is compatible with zero within errors. We will studyboth these GFFs in detail in section V.We now consider the behavior of the slopes of the GFFs A n and their relation to the transverse size of the nucleon.Since R − dx x n − H ( x, ξ = 0 , t ) = A n ( t ), it is evident that the GFFs for increasing n correspond to increasing averagemomentum fractions h x i . As the average momentum fraction gets larger, or equivalently as n → ∞ , we expect the t -slope of the GFFs A n to flatten. This may be understood in terms of the light cone Fock representation [42, 43] bythe fact that the final state nucleon wavefunction for a struck quark with momentum fraction x and initial transversemomentum k in ⊥ depends on the transverse momentum k ⊥ = k in ⊥ − (1 − x )∆ ⊥ . Hence, a large transverse momentumtransfer t = − ∆ ⊥ can be better absorbed without causing breakup of the bound state by quarks with large momentumfraction x . Additional insight is obtained by considering the impact parameter dependent GPD, H ( x, b ⊥ ), whichhas a probability interpretation and is the Fourier transform [6] with respect to transverse momentum transfer of H ( x, ξ = 0 , t = − ∆ ⊥ ): H ( x, b ⊥ ) = Z d ∆ ⊥ (2 π ) e − ib ⊥ · ∆ ⊥ H ( x, ξ = 0 , t = − ∆ ⊥ ) , (22)where ∆ ⊥ is the transverse momentum transfer. The impact parameter b ⊥ corresponds to the distance of the activequark from the center of momentum of the nucleon. As x →
1, a single quark will carry all the longitudinal momentumof the nucleon and therefore represent its center of momentum, so that the impact parameter distribution in this limitis strongly peaked around the origin, H ( x → , b ⊥ ) ∝ δ ( b ⊥ ). The corresponding flattening of the GFFs in themomentum transfer t is clearly visible in Fig. (10), where we compare the slopes of the GFFs A ( n =1 , , which havebeen normalized to unity at t = 0.Dipole fits to the GFFs in Fig. (10), denoted by the solid lines and statistical error bands, enable us to determinethe slopes of the form factors and thus express the two- and three-dimensional rms radii h r i / ⊥ and h r i / in termsof the dipole masses m D h r ⊥ i = 23 h r i = 8 m D . (23)3 m (cid:138)(cid:139) , u (cid:140) d m (cid:141)(cid:142) , u (cid:143) dm (cid:144)(cid:145) , u (cid:146) d m (cid:147)(cid:148) , u (cid:149) dm (cid:150)(cid:151) , u (cid:152) d m (cid:153)(cid:154) , u (cid:155) dm (cid:156)(cid:157) , u (cid:158) d m (cid:159)(cid:160) , u ¡ dm ¢£ , u ⁄ d m ¥ƒ , u § dm ¤' , u “ d m «‹ , u › d A fi A A fl A A (cid:176) A A – A A † A A ‡ A A · A A (cid:181) A A ¶ A A • A A ‚ A A „ A ” t » GeV … ‰ t (cid:190) GeV ¿ (cid:192) t ` GeV ´ ˆ t ˜ GeV ¯ FIG. 13: Ratio of generalized form factors A ( t ) /A ( t ) for the flavor combinations u − d (left) and u + d (right) comparedwith the parametrization in Ref. [11]. Disconnected contributions are not included. Since the range of values for the momentum transfer t is much smaller for the large volume ( L = 28 ) dataset, wehave restricted the dipole fits for all datasets to the overlapping region of t = 0 . . . − . . Our results for the2d rms radii versus the pion mass squared are presented in Figs. (11) and (12). These results confirm the dramaticdependence of the transverse rms radius on the moment and thus the average momentum fraction as first observed[14]for pion masses 750 MeV and higher, and show that this dependence increases as the pion mass decreases. Indeed,considering the ratio of the n = 3 moment to the n = 1 moment, which both correspond to the same sum or differenceof quarks and antiquarks, we observe that for vector GFFs this ratio decreases from approximately 0 .
58 to 0 .
22 as thepion mass decreases from 750 MeV to 350 MeV, and for axial vector GFFs, it decreases from roughly 0 .
71 to 0 . A ( t ) /A ( t )and e A ( t ) / e A ( t ) to the parametrization by Diehl et al.[11] as function of the momentum transfer squared t . As thepion mass decreases, the slope of our results approaches that of the phenomenological parametrization. Our resultsclearly indicate that a factorized ansatz for the GPDs in x and t , which would lead to constant ratios in Figs. (13)and (14) breaks down already for small values of the momentum transfer squared | t | ≪ .4 m ˘˙ , u ¨ d m (cid:201)˚ , u ¸ dm (cid:204)˝ , u ˛ d m ˇ— , u (cid:209) dm (cid:210)(cid:211) , u (cid:212) d m (cid:213)(cid:214) , u (cid:215) dm (cid:216)(cid:217) , u (cid:218) d m (cid:219)(cid:220) , u (cid:221) dm (cid:222)(cid:223) , u (cid:224) d m Æ(cid:226) , u ª dm (cid:228)(cid:229) , u (cid:230) d m (cid:231)Ł , u Ø d A Œ º A (cid:236) A (cid:237) (cid:238) A (cid:239) A (cid:240) æ A (cid:242) A (cid:243) (cid:244) A ı A (cid:246) (cid:247) A ł A ø œ A ß A (cid:252) (cid:253) A (cid:254) A (cid:255) (cid:1) A (cid:0) A (cid:2) (cid:3) A (cid:4) A (cid:5) (cid:6) A (cid:7) A (cid:8) (cid:9) A (cid:10) A (cid:11) (cid:12) A (cid:13) (cid:14) t (cid:15) GeV (cid:16) (cid:17) t (cid:18) GeV (cid:19) (cid:20) t (cid:21) GeV (cid:22) (cid:23) t (cid:24) GeV (cid:25) FIG. 14: Ratio of polarized generalized form factors e A ( t ) / e A ( t ) for the flavor combinations u − d (left) and u + d (right)compared with the parametrization in Ref. [11]. Disconnected contributions are not included. The GFFs A q ( t = 0) = h x i q and B q ( t = 0) enable us to compute the total quark angular momentum contributionto the nucleon spin [5], J q = 1 / A q (0) + B q (0)). Figures (15) and (16) show results for the quark spin e A q ( t =0) / q / L q = J q − ∆Σ q / S = 1 / q based on self-consistently improved one-loopChPT [44, 45, 46, 47, 48] for ∆Σ u + d and ChPT including the ∆ resonance [49, 50] for g A = ∆Σ u − d and are shownas shaded bands. The values for B q ( t = 0) have been obtained from a linear extrapolation of B u + d ( t ) and a dipoleextrapolation for B u,d ( t ) in t . The resulting uncertainty in B q ( t = 0), which contributes to the uncertainty in L q ,depends on the details of the corresponding fit, such as the functional form and range of t , and is therefore partiallysystematic. To allow the reader to assess the absolute statistical errors, we represent the errors for L q coming fromthe extrapolation in t by error bands around the m π -axis in Figs. (15) and (16). Experimental results for the quarkspin fractions ∆Σ u + d and ∆( u, d ) = ∆Σ u,d are represented by open stars for the prediction given in the HERMESpublication from 1999 [51] and filled stars for the 2007 HERMES results [52]. The significant difference between the5new HERMES results, which are consistent with recent COMPASS results[53], and the values given in [51] is probablyto a large extent due to the simple Regge-parametrization which has been used in [51] to compute the contributionto ∆Σ coming from the low x -region. It is gratifying that the new values are much closer to our lattice results.These results reveal two remarkable features of the quark contributions to the nucleon spin. The first is that themagnitude of the orbital angular momentum contributions of the up and down quarks, L u and L d , are separatelyquite substantial, starting at 0.15 at m π = 750 MeV and increasing to nearly 0.20 at 350 MeV, and yet they cancelnearly completely at all pion masses. The second is the close cancellation between the orbital and spin contributionsof the d quarks, L d and ∆Σ d / Π @ GeV D c on t r i bu t i on s t onu c l eon s p i n L u + d DS u + d (cid:144) Π @ GeV D c on t r i bu t i on s t onu c l eon s p i n FIG. 15: Total quark spin and orbital angular momentumcontributions to the spin of the nucleon. The filled andopen stars represent values given in HERMES 2007 [52] and1999 [51] respectively and open symbols represent earlierLHPC/SESAM calculations. The error bands are explainedin the text. Disconnected contributions are not included. Π @ GeV D - c on t r i bu t i on s t onu c l eon s p i n L u L d DS u (cid:144) DS d (cid:144) Π @ GeV D - c on t r i bu t i on s t onu c l eon s p i n FIG. 16: Quark spin and orbital angular momentum contri-butions to the spin of the nucleon for up and down quarks.Squares and triangles denote ∆Σ u and ∆Σ d respectively, anddiamonds and circles denote L u and L d respectively. Thefilled and open stars represent values given in HERMES 2007[52] and 1999 [51] respectively and open symbols representearlier LHPC/SESAM calculations. The error bands are ex-plained in the text. Disconnected contributions are not in-cluded. V. CHIRAL EXTRAPOLATIONS
Our ultimate goal is to use the combination of full QCD lattice calculations in the chiral regime and chiral pertur-bation theory to extrapolate to the physical pion mass, to extrapolate to infinite volume, to extrapolate in momentumtransfer, and to correct for lattice artifacts, with all the relevant low energy constants being determined solely fromlattice data. Significant progress has been made in many aspects of chiral perturbation theory (ChPT) relevant tothe nucleon observables addressed in this work [47, 50, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66]. Althoughimportant developments have been made in correcting for our hybrid action [61, 67, 68, 69, 70] and finite volume[50, 60], results for the relevant GFFs are not yet available. In this work we will focus on ChPT treatment of the pionmass and momentum dependence.The basic problem is that currently, there is not yet unambiguous evidence supporting a particular counting schemeand its convergence criteria, leading to a range of alternative re-summations, and there is similar ambiguity concerningthe choice of degrees of freedom, such as when and if it is essential to include the ∆ resonance. When complete resultsfor the observables of interest are available, it will be interesting to compare four approaches: heavy baryon ChPT(HBChPT)[62, 63, 64], covariant ChPT in the baryon sector (BChPT)[65], self-consistently improved one-loop ChPT[44, 45, 46, 47, 48], and ChPT with finite-range regulators[57, 59, 66]. Although self-consistent improvement byutilizing values of parameters like f π and g A calculated on the lattice at the relevant pion mass and finite-rangeregulators appear to improve the behavior of ChPT at larger values of the pion mass, based on the results availablein the literature, we will focus on the two formulations HBChPT and BChPT.Heavy baryon ChPT, which we will subsequently always refer to as HBChPT, assumes that m π and the magnitudeof the spatial three-momentum, p , are much smaller than the nucleon mass in the chiral limit, m N ∼
890 MeV,and the chiral scale Λ χ = 4 πf π = 1 .
17 GeV, and simultaneously expands in powers of the four quantities ǫ = { p Λ χ , m π Λ χ , pm N , m π m N } . In contrast, covariant baryon ChPT, which, slightly changing the notation of Ref. [65], we willsubsequently always refer to as CBChPT, does not treat m N and Λ χ as comparable scales, but rather keeps all powers6 f π [GeV] m N [GeV] g A h x i π, u + d c [GeV − ] c [GeV − ] c [GeV − ] g πN ∆ ∆ = m − m N [GeV]0 .
092 0 .
89 1 . / (2 / ) g A ( m N ) n generated by the couplings included in the ChPT Lagrangian. Thus, it is a resummation that includes termsthat would contribute in higher order to HBChPT and may be thought of as recoil corrections. The HBChPT results ofrefs. [62, 63, 64] and the CBChPT results of ref. [65] have the desirable property that they use the same regularizationscheme, so that truncation of the higher order terms in CBChPT yields the corresponding result in HBChPT, afeature that we will utilize below. One of our primary objectives will be to assess the regimes of applicability of thesetwo alternative formulations. For notational convenience, we will refer to the generic momentum in both theoriesas p , and order both theories in powers of p n . We would like to note, however, that the counting scheme of theHBChPT-approaches in [62, 63, 64] differs from the one used in the CBChPT-approach of [65].The HBChPT[62, 63, 64] and CBChPT[65] results for GFFs, including the dependence on the momentum transfersquared, t , enable us to investigate for the first time possible non-trivial correlations in the m π - and t -dependence ofGFFs. It is interesting to note that to order O ( p ) in HBChPT, unpolarized and polarized GFFs are independentof each other and depend on separate chiral limit values and counter-terms. In contrast, in CBChPT, the pion-massdependence of the isovector momentum fraction of quarks h x i u − d is simultaneously controlled by both the chirallimit values h x i u − d and h x i u − ∆ d , as is also the case for h x i ∆ u − ∆ d . To O ( p ), however, CBChPT does not includeinsertions from pion operators, and it turns out that the t -dependence for the isosinglet (and isovector) case to thisorder is therefore essentially linear and decouples from the pion mass dependence. Once CBChPT calculations havebeen pushed to higher orders, it will be interesting to study the combined non-analytic ( t, m π )-dependence of thefull set of polarized and unpolarized GFFs based on this approach. For the time being, we will investigate possiblenon-trivial correlations in m π and t in the framework of covariant CBChPT by including partial O ( p )-corrections asdiscussed below.The low energy constants used for the chiral extrapolations are summarized in Table III. Ultimately, all these valueswill be determined simultaneously by a global fit to a full set of lattice calculations of all the relevant observables usingthe same lattice action, but at present they are chosen as follows. We will use f π = 0 .
092 GeV, m N = 0 .
89 GeV and g A = 1 . f π , m N , and g A . In addition,we need the low energy constants c , c and c for the chiral extrapolation of C u + d in the framework of HBChPTin section V J. Here, c = 3 . − and c = − . − are taken from Refs. [65, 71] and c = − .
90 GeV − has been obtained from a CBChPT fit to our nucleon mass lattice data, which provides us with a parametrizationfor the pion mass dependence of the nucleon mass in our simulation needed for some of the chiral extrapolationsbelow. Depending on the order of ChPT, diagrams with insertions of pion operators contribute for the isosingletGFFs, which introduces the momentum fraction of quarks in the pion in the chiral limit, h x i π, u + d , as an additional lowenergy constant. Ultimately, we will calculate this quantity from chiral fits to hybrid lattice results for the pion, butfor now we use h x i π, u + d = 0 . . N c relation g πN ∆ = 3 / (2 / ) g A for the pion-nucleon-∆ coupling g πN ∆ . This latter result will soon be superseded by extrapolation of lattice calculations of the N − ∆ transition formfactor[78, 79, 80].Our chiral extrapolations are organized as follows, and summarized in Table IV. We begin in section V A with acomparison of CBChPT and HBChPT extrapolations of the isovector GFF A u − d ( t ) and show that whereas CBChPTyields a satisfactory fit over the range of pion masses used in the lattice calculations, HBChPT only produces fits to thelowest few points. Hence, in sections V B through V F, we study the pion mass and t -dependence of the isovector GFFs B u − d and C u − d and of the isoscalar GFFs A u + d , B u + d and C u + d . This is followed in section V G by a discussion of ourresults for the angular momentum of quarks, based on the CBChPT extrapolations for A ( t = 0) and B ( t = 0). Inthe counting scheme of [62, 63, 64] insertions of pion operators occur at O ( p ) in HBChPT, leading to a non-analyticcombined dependence on m π and t for the GFFs B u + d ( t ) and C u + d ( t ) , which we study in sections V H, V I, and V J.Finally, in section V K, we study the pion mass dependence of the total quark angular momentum J q in HBChPT,both including [55] and excluding explicit ∆ degrees of freedom and compare with the corresponding CBChPT results.The chiral extrapolations of B u + d ( t ) and C u + d ( t ) in sections V F, V I, and V J are the first parametrizations of theircombined ( m π , t )-dependence, providing valuable insights into non-trivial correlations of m π and t . Although we havea considerable amount of lattice data available for reasonably small values of | t | ≤ . and m π ≤ section GFF HBChPT CBChPT expected dependence on m π , t A A u − d O ( p ) O ( p ) non-analytic in m π , ≈ linear in t B B u − d O ( p ) + corr. of O ( p ) non-analytic in m π , ≈ linear in t C C u − d O ( p ) + corr. of O ( p ) non-analytic in m π , ≈ linear in t D A u + d O ( p ) + corr. of O ( p ) non-analytic in m π and t I, E B u + d O ( p ) O ( p ) + O ( p )-CTs non-analytic in m π and t F C u + d O ( p ) + corr. of O ( p , ) non-analytic in m π and t G J u + d = 1 / A + B ) u + d O ( p ) + corr. of O ( p )H E u + d = ( A + t/ (4 m N ) B ) u + d O ( p ) linear in m π and t H M u + d = ( A + B ) u + d O ( p ) non-analytic in m π and t J C u + d O ( p ) non-analytic in m π and t K J u + d = 1 / A + B ) u + d O ( p )K J u + d = 1 / A + B ) u + d O ( p ) with ∆TABLE IV: Overview of different approaches to the ( m π , t )-dependence of GFFs in ChPT studied in sections V A-V K. use an extended set of results for | t | < . , m π < A. CBChPT extrapolation of A u − d ( t ) The O ( p ) CBChPT result[65] for the isovector GFF A u − d ( t ) is A u − d ( t, m π ) = A ,u − d (cid:18) f u − dA ( m π ) + g A π f π h A ( t, m π ) (cid:19) + e A ,u − d j u − dA ( m π ) + A m π ,u − d m π + A t,u − d t , (24)where f u − dA ( m π ), h A ( t, m π ) and j u − dA ( m π ) contain the non-analytic dependence on the pion mass and momentumtransfer squared and A ,u − d ≡ A u − d ( t = 0 , m π = 0). Because of the small prefactor, the term ∝ h A ( t, m π ) isof O (10 − ) for m π ≤
700 MeV, | t | < and therefore numerically negligible. Thus, there are essentially nocorrelations of t and m π present, and the dependence on t is only due to the counter term ( A t,u − d t ). We use the value e A ,u − d = 0 .
17 obtained from a chiral fit to our lattice results for e A u − d ( t = 0) = h x i ∆ u − ∆ d [47]. Since the low energyconstant A ,u − d is a common parameter in the CBChPT-formulae for the GFFs A u − d , B u − d and C u − d , we performeda simultaneous fit based on Eq. (24) (for A u − d ), Eq. (25) (for B u − d ) and Eq. (26) (for C u − d ) with a total of 9 (1common and 8 separate) fit parameters to over 120 lattice datapoints. The details of the CBChPT-extrapolationsand the results for the GFFs B u − d and C u − d will be discussed below in sections V B and V C, respectively. Wefind A ,u − d = 0 . h x i u − d = A u − d ( t = 0 , m π, phys ) = 0 . h x i MRST2001 u − d = 0 . h x i CTEQ6 u − d = 0 . A u − d ( t ) onthe momentum transfer squared is presented in Figs. (19) and (20), where we again obtain a good description of thelattice data.To study the difference between HBChPT and CBChPT, we took the heavy baryon limit of CBChPT whilekeeping the same values of the fit parameters, and obtained the dotted line in Fig. (19). This curve only overlaps theCBChPT curve for m π < m π, phys and drops off sharply for m π > m π, phys , indicating the quantitative importanceof the truncated terms when using the coefficients from the CBChPT fit. In addition, it is important to ask theseparate question of how well the lattice data can be fit with the HBChPT expression when the coefficients aredetermined directly by a best fit to the data. The dashed curve in Fig. (17) shows the result of fitting our lattice datafor | t | < . and m π < . O ( p ). In the absence of the requisite full ChPT analysis, we have studied uncertainties in the chiralextrapolations by repeating the fit for different maximal values of the included pion masses. Figure (21) shows acomparison of the chiral extrapolations of A u − d , based on fits to the lattice data in the regions m π < Π @ GeV D A - d H t = . G e V L Π @ GeV D A - d H t = . G e V L FIG. 17: Lattice results for A u − d at t = 0 GeV versus m π . The error band is the result of a globalsimultaneous chiral fit using Eqs. (24), (25) and(26). The phenomenological result from CTEQ6 isindicated by the star. The heavy-baryon-limit ofthe CBChPT fit is shown by the dotted line, and aHBChPT fit to the lattice data for | t | < . and m π < . Π @ GeV D A - d H t »- . G e V L Π @ GeV D A - d H t »- . G e V L FIG. 18: Lattice results for A u − d at t ≈ − . versus m π together with the result of aglobal simultaneous chiral fit using Eqs. (24), (25)and (26). - t @ GeV D A - d H t L m Π > - t @ GeV D A - d H t L FIG. 19: Lattice results for A u − d ( t ) at m π ≈ . - t @ GeV D A - d H t L m Π > - t @ GeV D A - d H t L FIG. 20: Lattice results for A u − d ( t ) at m π ≈ . point is not driven by the large pion mass region, where O ( p ) corrections would be largest. Furthermore, all fourchiral fits are, within statistical errors, in agreement with the lattice data points at large pion masses. This indicatesthat the present statistical error envelope is comparable to any systematic effects due to higher order corrections.Another prescription to estimate O ( p ) corrections that has been advocated in the literature Ref. [65] is simplyadding a single m π term and, assuming both “naturalness” of the coefficient and the lack of other functional forms,seeing what error band arises from varying the coefficient from -1 to +1. Thus, to explore this possibility, followingRef. [65], we have added to the result in Eq. (24) the term δ (3) ,u − dA m π / (Λ χ m N ), where Λ χ ≈ . δ (3) ,u − dA in the range − , . . . , +1. The results of fits to the lattice datafor A u − d including this additional term are shown in Fig. (22), where the error band corresponds to δ (3) ,u − dA = 0, thedashed line corresponds to δ (3) ,u − dA = +1, and the dotted line corresponds to δ (3) ,u − dA = −
1. From the figure, we notethat this m π term alone with coefficients +1 and -1 is clearly inconsistent with the behavior of the data. Theoretically,this is not unreasonable, since the foundation of the IR regularization scheme is a resummation of classes of terms,and here a single cubic term has been arbitrarily singled out. Analogous fits with similar qualitatively inconsistentbehavior were also obtained for A u + d , treated in a later section. Hence, although this prescription may provide usefulestimates in other contexts, we do not believe it is useful in this work, and hence we do not include it in subsequentfits.9 Π @ GeV D A - d H t = G e V L Π @ GeV D A - d H t = G e V L FIG. 21: Lattice results for A u − d at t = 0 GeV versus m π together with chiral fits based onEq. (24). The four different error bands representchiral fits to lattice results including pion massesin the regions m π < Π @ GeV D A - d H t = L Π @ GeV D A - d H t = L FIG. 22: Chiral fits including the terms in Eq. (24)plus an additional m π term as described in the text.The narrow band is the original band omitting thisterm, the dashed line corresponds to δ (3) ,u − dA = 1,and the dotted line to δ (3) ,u − dA = − Π @ GeV D B - d H t »- . G e V L Π @ GeV D B - d H t »- . G e V L FIG. 23: Lattice results for B u − d at t ≈ − . versus m π together with the result of a globalsimultaneous chiral fit using Eqs. (24), (25) and(26). - t @ GeV D B - d H t L m Π > - t @ GeV D B - d H t L FIG. 24: Lattice results for B u − d ( t ) at m π ≈ − t ) together with the result of a globalsimultaneous chiral fit using Eqs. (24), (25) and(26). B. CBChPT extrapolation of B u − d ( t ) The O ( p ) CBChPT calculation [65] for the isovector B GFF gives B u − d ( t, m π ) = m N ( m π ) m N B ,u − d + A ,u − d h u − dB ( t, m π ) + m N ( m π ) m N (cid:26) δ t,u − dB t + δ m π ,u − dB m π (cid:27) , (25)where m N ( m π ) is the pion mass dependent nucleon mass, B ,u − d ≡ B u − d ( t = 0 , m π = 0), and where we haveincluded estimates of O ( p )-corrections in form of ( δ t,u − dB t ) and ( δ m π ,u − dB m π ). The low energy constants B ,u − d , δ t,u − dB and δ m π ,u − dB are treated as free parameters and may be obtained from a fit to the lattice data. The non-analyticdependence on m π and t is given by h u − dB ( t, m π ), but it turns out that this function is approximately independentof t , h u − dB ( t, m π ) ≈ h u − dB ( m π ) for m π ≤
700 MeV, | t | < . The t -dependence is therefore in practice lineardue to the O ( p )-correction term. For the m π -dependent nucleon mass we use O ( p ) CBChPT [65, 82] fitted to ourlattice results for m N . The chiral extrapolation of B u − d is based on a global simultaneous fit as discussed in theprevious section with A ,u − d as common fit parameter in Eqs. (24), (25) and (26). We obtain B ,u − d = 0 . B u − d ( t = 0 , m π, phys ) = 0 . C. CBChPT extrapolation of C u − d ( t ) The pion mass dependence of the isovector GFF C in CBChPT to O ( p ) is very similar to that of the isovector B above and given by [65] C u − d ( t, m π ) = m N ( m π ) m N C ,u − d + A ,u − d h u − dC ( t, m π ) + m N ( m π ) m N (cid:26) δ t,u − dC t + δ m π ,u − dC m π (cid:27) , (26)where C ,u − d ≡ C u − d ( t = 0 , m π = 0). As in the case of B u − d , ( δ t,u − dC t ) and ( δ m π ,u − dC m π ) represent O ( p )-correction Π @ GeV D - - - C - d H t »- . G e V L Π @ GeV D - - - C - d H t »- . G e V L FIG. 25: Lattice results for C u − d at t ≈ − . versus m π together with the result of a globalsimultaneous chiral fit using Eqs. (24), (25) and(26). - t @ GeV D - - C - d H t L m Π > - t @ GeV D - - C - d H t L FIG. 26: Lattice results for C u − d ( t ) at m π ≈ − t ) together with the result a globalsimultaneous chiral fit using Eqs. (24), (25) and(26). terms, and it turns out that h u − dC ( t, m π ) is practically independent of t . From a global simultaneous chiral fit based onEqs. (24), (25) and (26) with common fit parameter A ,u − d as discussed in section V A, we obtain C ,u − d = − . C u − d ( t = 0 , m π, phys ) = − . C u − d ( t ) is roughly one order of magnitude smaller than A u − d ( t ) and B u − d ( t ) and fully compatible with zerowithin errors, C u − d ( t ) ≈
0. This implies a rather mild dependence of the n = 2 moment of the GPDs H u − d ( x, ξ, t )and E u − d ( x, ξ, t ) on the longitudinal momentum transfer ξ (at least for small ξ ), so that H n =2 u − d ( ξ, t ) ≈ A u − d ( t ) and E n =2 u − d ( ξ, t ) ≈ B u − dn ( t ). D. CBChPT extrapolation of A u + d ( t ) The (total) isosinglet momentum fraction of quarks, A u + d ( t = 0) = h x i u + d is not only an important hadron structureobservable on its own but is in addition an essential ingredient for the computation of the total angular momentumcontribution of quarks to the nucleon spin, J u + d = 1 / A u + d (0) + B u + d (0)). The combined ( t, m π )-dependence inCBChPT is given by [65]: A u + d ( t, m π ) = A ,u + d (cid:18) f u + dA ( m π ) − g A π f π h A ( t, m π ) (cid:19) + A m π ,u + d m π + A t,u + d t + ∆ A u + d ( t, m π ) + O ( p ) , (27)where A ,u + d ≡ A u + d ( t = 0 , m π = 0), f u + dA ( m π ) and h A ( t, m π ) contain the non-analytic dependence on the pionmass and momentum transfer squared, and the constants A m π ,u + d and A t,u + d may be obtained from a fit to thelattice data. In this counting scheme, contributions from operator insertions in the pion line proportional to themomentum fraction of quarks in the pion in the chiral limit, h x i π, u + d , are of order O ( p ). However, in order to seeif such contributions could be relevant for the pion masses and values of the momentum transfer squared accessiblein our calculation, we include the estimate of the O ( p )-contribution ∆ A u + d provided in [65] in the fit to the latticedata points.Similar to the isovector case discussed in the previous sections, the low energy constant A ,u + d is a commonparameter in the chiral extrapolation formulae for the isosinglet GFFs A u + d , B u + d and C u + d . Using h x i π, u + d = 0 . B u + d and C u + d we refer to the sections V E and V F below. The chiral fit gives A ,u + d = 0 . h x i u + d = A u + d ( t = 0 , m π, phys ) = 0 . h x i MRST2001 u + d = 0 . h x i CTEQ6 u + d = 0 . h x i π, u + d by ±
10% only leads to a small change in A ,u + d ( t = 0) of O (1%), which is significantly smaller than the statistical error of ≈ m π , andtherefore the good agreement with the phenomenological value, is due to the O ( p )-contribution ∆ A u + d . It has to beseen if this somewhat unusual curvature persists once the full O ( p ) contribution is available and fitted to the latticeresults. The inclusion of contributions from disconnected diagrams could also require a different extrapolation in m π .The dependence of A u + d ( t ) on t at fixed values of m π is presented in Figs. (29) and (30).As in the case of A u − d , we also consider the heavy baryon limit of the CBChPT fit, giving the result A u + d ( t =0 , m π ) = A ,u + d + A m π ,u + d m π represented by the dotted line in Fig. (27), which agrees with the CBChPT result onlyover a very limited range at low pion masses. Notably, while the lattice results for A u + d are rising for larger pionmasses, the heavy-baryon-limit curve has the opposite slope with negative A m π ,u + d . However, a direct HBChPT fitwith free coefficients (see also section V H) shown by the dashed curve leads to a positive A m π ,u + d and a reasonabledescription of the lattice datapoints. Π @ GeV D A + d H t = . G e V L Π @ GeV D A + d H t = . G e V L FIG. 27: Lattice results for A u + d at t = 0 GeV versus m π . The error band is the result of a globalsimultaneous chiral fit using Eqs. (27), (28) and(29). The phenomenological value from CTEQ6 isdenoted by a star. The heavy-baryon-limit of theCBChPT fit is shown by the dotted line, and aHBChPT fit to the lattice data for | t | < . and m π < . Π @ GeV D A + d H t »- . G e V L Π @ GeV D A + d H t »- . G e V L FIG. 28: Lattice results for A u + d ( t ) at t ≈ − . versus m π together with the result of aglobal simultaneous chiral fit using Eqs. (27), (28)and (29). - t @ GeV D A + d H t L m Π > - t @ GeV D A + d H t L FIG. 29: Lattice results for A u + d ( t ) at m π ≈ . - t @ GeV D A + d H t L m Π > - t @ GeV D A + d H t L FIG. 30: Lattice results for A u + d ( t ) at m π ≈ . E. CBChPT extrapolation of B u + d ( t ) The dependence on the pion mass and the momentum transfer squared of the isosinglet B GFF in O ( p ) CBChPTis given by [65] B u + d ( t, m π ) = m N ( m π ) m N B ,u + d + A ,u + d h u + dB ( t, m π ) + ∆ B u + d ( t, m π ) + m N ( m π ) m N (cid:26) δ t,u + dB t + δ m π ,u + dB m π (cid:27) + O ( p ) , (28)where B ,u + d ≡ B u + d ( t = 0 , m π = 0), and the terms ∆ B , δ t,u + dB t and δ m π ,u + dB m π are of O ( p ) and represent onlya part the full O ( p ) contribution. The a priori unknown constants B ,u + d , δ t,u + dB and δ m π ,u + dB may be obtained froma fit to the lattice data. A fit to our lattice results based on Eq. (28) turns out to be unstable and produces largevalues for the counter term parameter δ m π ,u + dB ≈
15. This can be seen as indication that other higher order correctionterms of O ( p ) not yet included in Eq. (28) are numerically important and needed to stabilize the extrapolation. Wenote that the counting scheme of [65] suggests that ∆ B is not a dominant O ( p )-contribution concerning the pionmass dependence, at least for t = 0. This can be seen to some extent from the heavy-baryon-limit of Eq. (28), whichdoes not reproduce the full coefficient, ∝ ( A + B ) ,u + d , of the m π log( m π )-term in HBChPT (see e.g. [55, 62, 63]),but rather gives a term ∝ A ,u + d m π log( m π ) without the B ,u + d m π log( m π ) contribution. Π @ GeV D - - B + d H t »- . G e V L Π @ GeV D - - B + d H t »- . G e V L FIG. 31: Lattice results for B u + d at t ≈ − . versus m π together with the result of a globalsimultaneous chiral fit based on Eqs. (27), (29) anda variant of Eq. (28), as described in the text. - t @ GeV D - - - B + d H t L m Π > - t @ GeV D - - - B + d H t L FIG. 32: Lattice results for B u + d ( t ) at m π ≈ − t ) together with the result of a globalsimultaneous chiral fit based on Eqs. (27), (29) anda variant of Eq. (28), as described in the text. Since the instability of the fit can be traced back to the term ∆ B u + d ( t, m π ), we performed the final fit droppingthis contribution but keeping the counter terms ∝ t and ∝ m π . Based on this approach, we find that a globalsimultaneous fit to the GFFs A u + d , B u + d and C u + d , using Eqs. (27) and (29) as described in the previous section,leads to a stable chiral extrapolation of all three GFFs. In particular the counter term parameters δ t,u + dB and δ m π ,u + dB turn out to be very small and fully compatible with zero within errors. We obtain B ,u + d = − . B u + d, ( t = 0 , m π, phys ) = − . O ( p ) contributions to B u + d will lead to a qualitatively differentdependence on t and m π in the region where lattice results are available, so that the results above should to be takenwith due caution.In section V H, we will study the combined ( t, m π )-dependence of B u + d in HBChPT at O ( p ) [64]. One-loop graphswith insertions of pion operators are fully included, leading to a non-analytic dependence on the pion mass and themomentum transfer that is quite different from that shown in Figs. (31) and (32). F. CBChPT extrapolation of C u + d ( t ) The ( t, m π )-dependence of the isosinglet GFF C in CBChPT to O ( p ) is given by [65] C u + d ( t, m π ) = m N ( m π ) m N C ,u + d + A ,u + d h u + dC ( t, m π ) + ∆ C u + d ( t, m π ) + O ( p ) , (29)where C ,u + d ≡ C u + d ( t = 0 , m π = 0), and the term ∆ C u + d ∝ h x i π, u + d is a part of the full O ( p )-corrections [65].In this counting scheme, counter terms of the form δ t,u + dC t and δ m π ,u + dC m π first appear at O ( p ). In order to get3a first idea about the possible t - and m π -dependence of C u + d in CBChPT, we have included the formally higherorder counter terms δ t,u + dC t and δ m π ,u + dC m π in the fit to our lattice data, resulting in a stable chiral extrapolation.With h x i π, u + d = 0 . A u + d and Π @ GeV D - - - C + d H t »- . G e V L Π @ GeV D - - - C + d H t »- . G e V L FIG. 33: Lattice results for C u + d at t ≈ − . versus m π together with the result of a globalsimultaneous chiral fit using Eqs. (27), (28) and(29). - t @ GeV D - - - - - C + d H t L m Π > - t @ GeV D - - - - - C + d H t L FIG. 34: Lattice results for C u + d ( t ) at m π ≈ − t ) together with the result of a globalsimultaneous chiral fit using Eqs. (27), (28) and(29). B u + d discussed in sections V D and V E, respectively, we obtain from a global simultaneous fit C ,u + d = − . C u + d ( t = 0 , m π, phys ) = − . δ t,u + dC and δ m π ,u + dC turn out to be small. Changing h x i π, u + d by ±
10% results in a variation of C ,u + d ( t = 0)by ± n = 2 moment of the isosinglet GPDs H and E on the longitudinal momentum transfer ξ .A three-dimensional plot showing the combined ( t, m π )-dependence of C u + d is presented in Fig. (35), where theerror bars of the lattice data points are illustrated by the stretched cuboids. The lattice data are superimposed withthe result from the chiral fit discussed above, which is shown as a surface. The statistical error bars of the fit areshown for clarity only as bands for t = 0 and m π = 0, respectively. In section V J below, we will compare the resultsbased on CBChPT with a fit to C u + d in the framework of HBChPT. G. Quark angular momentum J in CBChPT The forward limit values of the isovector and isosinglet GFFs A ( t = 0) and B ( t = 0) we have studied insections V A, V B, V D and V E allow us to compute the angular momentum contributions of up- and down-quarksto the spin of the nucleon, J q = 1 / A q (0) + B q (0)) = 1 / h x i q + B q (0)). From the separate chiral extrapolationsof the isosinglet A and B in CBChPT, we find for the total u + d quark angular momentum at the physicalpion mass J u + d ( m π, phys ) = 0 . S = 1 /
2. Together with thechiral extrapolations for the isovector u − d combination, we obtain the surprising result that the total quark angularmomentum is carried by the up-quarks, J u ( m π, phys ) = 0 . J d ( m π, phys ) = − . d / L d appears to be systematicfor all m π , and it will be interesting to understand whether this is accidental or has a physical origin. Taking intoaccount preliminary results for the quark spin e A q / t = 0) = ∆Σ q /
2, as obtained from a ChPT extrapolationincluding the ∆ resonance [49, 50] of g A = ∆Σ u − d and a self-consistently improved one-loop ChPT extrapolation of∆Σ u + d [48], we find that the quark orbital angular momentum L q = J q − ∆Σ q / L u = − . L d = 0 . m π, phys , wherewe find L u + d = 0 . h i ∆ q , h x i q and h x i ∆ q have been included in the extrapolations, and that we have omitted disconnected diagrams in the lattice calculations.We will compare these CBChPT results with corresponding HBChPT results below in section V J.4 Π @ GeV D - t @ GeV D - - - - Π @ GeV D FIG. 35: Combined ( t, m π )-dependence of C u + d from a global simultaneous chiral fit using Eqs. (27), (28) and (29) comparedto lattice data. H. HBChPT extrapolation of E u + d and M u + d In heavy baryon chiral perturbation theory [62, 64] to O ( p ), the combined ( t, m π )-dependence of the GFF-combination E u + d ( t ) = A u + d ( t ) + t/ (4 m N ) B u + d ( t ) is quite different from that of M u + d ( t ) = A u + d ( t ) + B u + d ( t ),which in the forward limit is equal to two times the total quark angular momentum 2 J q = M u + d ( t = 0). While atthis order M u + d shows a non-analytic dependence on t and m π as discussed below, E u + d is constant up to analytictree-level contributions, E u + d ( t, m π ) = E ,u + d + E m π ,u + d m π + E t,u + d t. (30)A fit to our lattice results based on Eq. (30) is shown in Figs. (36), (37), (38) and (39). The linear dependence Π @ GeV D H A + d + t (cid:144) H M N L B + d LH t = . G e V L Π @ GeV D H A + d + t (cid:144) H M N L B + d LH t = . G e V L FIG. 36: Lattice results for E u + d at t = 0 versus m π together with the result of a global chiral fitusing Eq. (30). Π @ GeV D H A + d + t (cid:144) H M N L B + d LH t »- . G e V L Π @ GeV D H A + d + t (cid:144) H M N L B + d LH t »- . G e V L FIG. 37: Lattice results for E u + d at − t ≈ . versus m π together with the result of a globalchiral fit using Eq. (30). on t and m π works well even beyond the fitted range, i.e. for − t ≥ .
48 GeV . This is not surprising since E u + d is clearly dominated by the GFF A u + d , which does not show a strong curvature in t as seen in Figs. (8) and (10).However, it is obvious that the HBChPT result in Eq. (30) lacks structures which would allow for an upwards bending5of E u + d ( t = 0) = A u + d ( t = 0) at small pion masses towards the phenomenological value, in contrast to the covariantapproach studied in section V A. The fit gives E ,u + d = 0 . - t @ GeV D H A + d + t (cid:144) H M N L B + d LH t L m Π > - t @ GeV D H A + d + t (cid:144) H M N L B + d LH t L FIG. 38: Lattice results for E u + d at m π ≈ − t together with the result of a globalchiral fit using Eq. (30). - t @ GeV D H A + d + t (cid:144) H M N L B + d LH t L m Π > - t @ GeV D H A + d + t (cid:144) H M N L B + d LH t L FIG. 39: Lattice results for E u + d at m π ≈ − t together with the result of a globalchiral fit using Eq. (30). chiral extrapolations based on HBChPT of the total angular momentum and the anomalous gravitomagnetic moment B u + d ( t = 0) below. At the physical pion mass, we find E u + d ( t = 0 , m π, phys ) = h x i u + d = 0 . h x i CTEQ6 u + d = 0 . h x i MRST2001 u + d = 0 . M u + d ( t ) for non-zero t is given by [62, 64] M u + d ( t, m π ) = M ,u + d (cid:26) − g A m π (4 πf π ) ln (cid:18) m π Λ χ (cid:19) (cid:27) + M (2 ,π )2 ( t, m π ) + M m π ,u + d m π + M t,u + d t, (31)with new counter terms M m π ,u + d and M t,u + d . The non-analytic dependence on t and m π in M (2 ,π )2 ( t, m π ) resultsfrom pion-operator insertions and is directly proportional to the (isosinglet) momentum fraction of quarks in the pionin the chiral limit, h x i π, u + d . We use h x i π, u + d = 0 . M ,u + d = 0 . Π @ GeV D H A + d + B + d LH t »- . G e V L Π @ GeV D H A + d + B + d LH t »- . G e V L FIG. 40: Lattice results for M u + d at | t | ≈ . versus m π together with the result of a globalchiral fit using Eq. (31). - t @ GeV D H A + d + B + d LH t L m Π > - t @ GeV D H A + d + B + d LH t L FIG. 41: Lattice results for M u + d at m π ≈ − t ) together with the result of aglobal chiral fit using Eq. (31). and M u + d ( t = 0 , m π, phys ) = 0 . I. HBChPT extrapolation of B u + d Total momentum and angular momentum conservation implies that the total anomalous gravitomagnetic momentof quarks and gluons in the nucleon has to vanish, P q,g B ( t = 0) = 0. An interesting question is whether theindividual quark and gluon contributions to B are separately zero or very small, as previously speculated [83, 84].6The first lattice calculations [12, 13] showed that B u + d is compatible with zero for ( − t ) ≥ . at pion massesof m π ≥
600 MeV. Based on our new results and the ChPT fits performed above, we are now in a position to study B u + d more carefully as a function of t and m π . The GFF B u + d can be written as a linear combination of Eq. (30) andEq. (31). A separate fit to the data with fixed E ,u + d = 0 . B u + d ( t = 0 , m π, phys ) = 0 . M u + d and E u + d above that in combination give ( M − E ) u + d ( t = 0 , m π, phys ) = 0 . B u + d ( t = 0) is again rather small, we note that the sign is different from that foundin section V E based on the CBChPT fit. A 10% variation of the input parameter h x i π, u + d leads to change of 0.023 in B u + d ( t, m π, phys ) at t = 0, which is well below the statistical error of 0.049, and a change of 0.008 at a momentumtransfer of | t | ≈ .
24 GeV , which is well below the statistical error of 0.031.The dependence of B u + d on t and on the pion mass is shown in Figs. (42) and (43). The dependence on the momen-tum transfer squared turns out to be somewhat different from the CBChPT result in Fig. (32) where contributionsfrom pion operator insertions ∝ h x i π, u + d have not been included, but the two results are statistically consistent. The Π @ GeV D - - B + d H t »- . G e V L Π @ GeV D - - B + d H t »- . G e V L FIG. 42: Lattice results for B u + d at | t | ≈ . versus m π together with the result of a globalchiral fit using a linear combination of Eqs. (30) and(31). - t @ GeV D - - B + d H t L m Π > - t @ GeV D - - B + d H t L FIG. 43: Lattice results for B u + d at m π ≈
350 MeVversus ( − t ) together with the result of a globalchiral fit using a linear combination of Eqs. (30)and (31). combined dependence of B u + d on t and m π from the HBChPT fit compared to lattice data points is presented ina 3d-plot in Fig. (44). It is interesting to note that the fit based on Eq. (31) leads to a clearly visible non-analyticdependence of B u + d on the pion mass and the momentum transfer. In particular, we find a non-zero, negative B u + d for | t | > . , m π / . from the chiral extrapolation, which may be confirmed in future lattice calculationsor experimental measurements. J. HBChPT extrapolation of C u + d At order O ( p ), the pion mass dependence of the GFF C u + d ( t ) is given by [62, 63, 64] C u + d ( t, m π ) = 11 − t/ (4 m N ) (cid:26) C ,u + d + E (1 ,π )2 ( t, m π ) + E (2 ,π )2 ( t, m π ) + C m π ,u + d m π + C t,u + d t (cid:27) , (32)where C ,u + d ≡ C u + d ( t = 0 , m π = 0). The terms E (1 ,π )2 ( t, m π ) and E (2 ,π )2 ( t, m π ) contain non-analytic terms in t and m π that come from insertions of pion operators proportional to h x i π, u + d . Additionally, E (2 ,π )2 ( t, m π ) depends on the lowenergy constants c , c and c . We fix these parameters to the values given in Table III. The result of a fit to our latticedata as a function of the pion mass squared for fixed t ≈ − .
24 GeV is shown in Fig. 45. The t -dependence at a pionmass of ≈
350 MeV is presented in Fig. (46). We find C ,u + d = − . C u + d ( t = 0 , m π, phys ) = − . C u + d for ( − t ) → h x i π, u + d by 10% results in a ≈
5% change of C ,u + d ( t = 0), which is belowthe 11% statistical error.Fig. (47) shows our combined lattice results for C u + d versus ( − t ) and m π in a single three-dimensional plot, togetherwith the result from the HBChPT fit discussed above, which is shown as a surface. The statistical error bars are7 Π @ GeV D - t @ GeV D - Π @ GeV D FIG. 44: Combined ( t, m π )-dependence of the quark anomalous gravitomagnetic moment B u + d from a global chiral fit using alinear combination of Eqs. (30) and (31) compared to lattice data. Π @ GeV D - - - - - C + d H t »- . G e V L Π @ GeV D - - - - - C + d H t »- . G e V L FIG. 45: Lattice results for C u + d at | t | ≈ . versus m π together with the result of a globalchiral fit using Eq. (32). - t @ GeV D - - - - - C + d H t L m Π > - t @ GeV D - - - - - C + d H t L FIG. 46: Lattice results for C u + d at m π ≈
350 MeVversus − t together with the result of a global chiralfit using Eq. (32). shown for clarity as bands for t = 0 GeV and m π = 0 MeV, respectively. It is interesting to note that the overallshape of the extrapolation surface is similar to the CBChPT result in section V F. The only behavior that differsby more than the statistical errors is the slightly stronger bending towards negative values of C u + d at the origin inFig. (47). K. HBChPT extrapolation of quark angular momentum J From our results for M above, we find a total quark angular momentum J u + d = 0 . J u + d by first extrapolating B ( t, m π ) to t = 0 and combining it with A ( t = 0 , m π ) to obtain J ( m π ), and then extrapolating the values of J ( m π ) to the physical pion mass using HBChPT8 Π @ GeV D - t @ GeV D - - Π @ GeV D FIG. 47: Combined ( t, m π )-dependence of C u + d from a global chiral fit using Eq. (32) compared to lattice data. that explicitly includes the ∆ resonance [55]. Evaluating Eq. (31) at t = 0 yields J u + d ( m π ) = 12 (cid:26) ( A + B ) ,u + d + 3 (cid:18) h x i π, u + d − ( A + B ) ,u + d (cid:19) g A m π (4 πf π ) ln (cid:18) m π Λ χ (cid:19)(cid:27) + J m π ,u + d m π , (33)which agrees with [55]. Note that in the notation of [55] we have b qN = ( A + B ) ,u + d = M ,u + d . Including explicitly Π @ GeV D J u + d Π @ GeV D J u + d FIG. 48: Chiral extrapolation of J u + d usingHBChPT including the ∆ resonance, Eq. (34).The fit and error band on the axis are explainedin the text. Π @ GeV D J u + d Π @ GeV D J u + d FIG. 49: Comparison of a global fit to A u + d + B u + d using Eq. (31) and the lattice results for J u + d . Thefit and error band on the axis are explained in thetext. the ∆ resonance in the calculation, the ChPT result then reads [55] J u + d ( m π ; ∆) = J u + d ( m π ) − (cid:18)
92 ( A + B ) ,u + d + 3 h x i π, u + d − b q ∆ (cid:19) × g πN ∆ πf π ) ( ( m π − ) ln (cid:18) m π Λ χ (cid:19) + 2∆ p ∆ − m π ln ∆ − p ∆ − m π ∆ + p ∆ − m π !) , (34)9where ∆ = m ∆ − m N denotes the ∆-nucleon mass difference. In order to reduce the number of free parametersin the fit, we use ∆ = 0 . N c relation g πN ∆ = 3 / (2 / ) g A from Table III. We first extrapolate B u + d ( t, m π ) linearly in t to t = 0. A fit based on Eq. (34) to the lattice results for ( A + B ) u + d ( t = 0 , m π ) with m π ≤
700 MeV then gives ( A + B ) ,u + d = b qN = 0 . b q ∆ = 0 . J u + d ( m π, phys ; ∆) = 0 . J u + d ( m π, phys ) = 0 . J u + d ( m π ; ∆) as a function of the pion mass is shown in Fig. (48), where the error due to the linearextrapolation of B ( t, m π ) to t = 0 is indicated by the error band at the m π -axis as explained at the end of sectionIV and the error bars on the lattice data points for J only include the errors arising from A ( t = 0). In Fig. (49),we compare the result of the chiral extrapolation of M u + d ( t ) for non-zero t from the previous section with the latticeresults for J u + d corresponding to the extrapolated B u + d ( t = 0). The two different ways of fitting and extrapolating J u + d in m π are compatible within errors, J u + d from t =0 = 0 . J u + dt → ( m π, phys ; ∆) = 0 . m π and therefore to a smaller central value for J u + d at the physical point. Together with preliminary results for the quark spin e A u + d / t = 0) = ∆Σ u + d / L u + d = J u + d − ∆Σ u + d / L u + d = 0 . J u + dt → ( m π, phys ; ∆) = 0 . L u + d = 0 . J u + d from t =0 = 0 . VI. SUMMARY AND CONCLUSIONS
This work presents the first comprehensive full lattice QCD study of the lowest three moments of unpolarized andpolarized GPDs in the chiral regime with pion masses as low as 350 MeV. We find good overall agreement with existingexperimental results. We note, however, that we have omitted disconnected diagrams, which in principle contributeto isoscalar observables.As in our previous study of the axial vector coupling constant [85], the consistency of these moments in latticevolumes of (2 . and (3 . at m π = 350 MeV indicates that finite volume effects are negligible within statisticalerrors.One significant result of this work is the clear indication that the transverse size of the nucleon, as characterizedby the transverse 2d rms radius, is a strongly decreasing function of the momentum fraction x carried by the quarks.At the lightest quark mass, the isovector transverse rms radius drops by almost 60% between the zeroth moment,which roughly corresponds to an average momentum fraction[86] h x i ≈ . h x i ≈ .
4. This decrease in the chiral regime is even stronger than ouroriginal observation of the decrease of the transverse size with momentum fraction in the “heavy pion world”[14].In a first direct comparison with phenomenological parametrizations of the GPDs H ( x, ξ = 0 , t ) and ˜ H ( x, ξ = 0 , t )constrained by structure function and form factor data[11], we find qualitative consistency for the ratios of GFFs inboth the isovector and isosinglet cases.Our results provide insight into the contributions of the spin and orbital angular momentum of u and d quarks tothe spin of the proton. Although the individual orbital angular momentum contributions of the u and d quarks aresizeable, L d ≈ − L u ≈ L u + d ≈
0. In addition, the spinand orbital contributions of the d quark also cancel within errors, so J d ≈
0. The total quark angular momentumcontribution is J u + d ≈ −
50% at our lowest pion masses. A u − d B u − d C u − d A u + d B u + d C u + d J u + d J u J d covariant BChPT 0 . . − . . − . − . . . − . . . − . . . . . t = 0 and m π, phys from chiral extrapolations, in the MS scheme at a scale of µ = 4 GeV . More quantitatively, we performed a variety of chiral fits to the unpolarized n = 2 moments using covariant BChPT[65] and HBChPT results [55, 62, 63, 64]. A summary of our results for various observables at the physical pion massand vanishing momentum transfer t = 0 is given in Table V, where the quoted errors are statistical only. We note thata consistent inclusion of all O ( p ) terms in our fits, which may involve additional constants, would have the potentialto increase the statistical errors on the physical quantities at the chiral limit. However, for the reasons described insection V A, the prescription of adding a single m π term with extremal values of “natural” coefficients does not leadto reasonable χ fits, so we have not attempted to include quantitative estimates of these errors in the present work.0With the exception of a fit to J q based on HBChPT including the ∆ resonance, we have consistently extrapolatedthe GFFs simultaneously in the pion mass and the momentum transfer squared t . The simultaneous covariant BChPTfits to the GFFs A , B and C , which include approximately 120 lattice data points and typically 9 unknownlow energy constants, produce reasonable parametrizations of the ( t, m π )-dependences of the generalized form factorsin the ranges m π ≤ | t | ≤ . . In particular, the covariant extrapolations for the isovectorand isosinglet momentum fractions h x i yield values at the physical point remarkably close to phenomenology. Thisrepresents a significant advance in our understanding of the pion mass dependence of these important observables.The first exploration of the combined non-analytic dependence of the isosinglet GFFs B u + d and C u + d on t and m π was made using covariant BChPT and HBChPT, and visualizations of the resulting ( t, m π )-dependences of theseGFFs in 3d-plots reveal interesting non-linear correlations in the pion mass and the momentum transfer squared.In spite of the overall success of the chiral extrapolations, it is clear that the ChPT analysis has not yet beencarried to sufficiently high order to be applicable to the full range of pion masses included in this work. The facts thatHBChPT fits to A u − d and A u + d cannot describe the behavior of the lattice data over as large a range or as accuratelyas covariant BChPT, and that the fitted counter terms are so different indicate that the higher order terms in 1 /m N included in CBChPT are important for these observables. Similarly, the fact that it was important to include someparticular terms of order O ( p ) and O ( p ) in some CBChPT fits indicates the need to fully determine these orders inChPT so that they can be consistently included in fits to lattice data. Finally, the significant effect of including the ∆in J u + d , its known importance in the axial charge, and large N c arguments[87] indicate the desirability of consistentinclusion of the ∆. Thus, future progress requires both the extension of lattice calculations to lower pion mass andthe inclusion of higher order terms and ∆ degrees of freedom in ChPT. Acknowledgments
The authors wish to thank M. Diehl, M. Dorati, T. Gail, T. Hemmert, J.-W. Chen and A. Manohar for stim-ulating discussions. This work was supported in part by U.S. DOE Contract No. DE-AC05-06OR23177 underwhich JSA operates Jefferson Laboratory, by the DOE Office of Nuclear Physics under grants DE-FG02-94ER40818,DE-FG02-04ER41302, DE-FG02-96ER40965, the Jeffress Memorial Trust grant J-813, the DFG (ForschergruppeGitter-Hadronen-Ph¨anomenologie) and in part by the EU Integrated Infrastructure Initiative Hadron Physics (I3HP)under contract number RII3-CT-2004-506078. Ph. H. and B. M. acknowledge support by the Emmy-Noether programof the DFG. Ph. H. and W. S. would like to thank the A.v. Humboldt-foundation for support by the Feodor-Lynenprogram. It is a pleasure to acknowledge the use of computer resources provided by the DOE through the USQCDproject at Jefferson Lab and ORNL and through its support of the MIT Blue Gene/L. We are indebted to membersof the MILC and SESAM Collaborations for providing the dynamical quark configurations which made our full QCDcalculations possible.1
APPENDIX A: TABLES
Lattice parameters of the datasets 1 , . . . , µ = 4 GeV . -t[GeV ˜ A u − d B u − d A u − d B u − d C u − d A u − d B u − d ˜ A u − d B u − d A u − d B u − d C u − d A u − d B u − d -t[GeV ˜ A u − d B u − d A u − d B u − d C u − d A u − d B u − d ˜ A u − d B u − d A u − d B u − d C u − d A u − d B u − d ˜ A u − d B u − d A u − d B u − d C u − d A u − d B u − d -t[GeV ˜ A u − d B u − d A u − d B u − d C u − d A u − d B u − d ˜ A u + d B u + d A u + d B u + d C u + d A u + d B u + d ˜ A u + d B u + d A u + d B u + d C u + d A u + d B u + d -t[GeV ˜ A u + d B u + d A u + d B u + d C u + d A u + d B u + d ˜ A u + d B u + d A u + d B u + d C u + d A u + d B u + d ˜ A u + d B u + d A u + d B u + d C u + d A u + d B u + d -t[GeV ˜ A u + d B u + d A u + d B u + d C u + d A u + d B u + d ˜ e A u − d e B u − d e A u − d e B u − d e A u − d e B u − d ˜ e A u − d e B u − d e A u − d e B u − d e A u − d e B u − d -t[GeV ˜ e A u − d e B u − d e A u − d e B u − d e A u − d e B u − d ˜ e A u − d e B u − d e A u − d e B u − d e A u − d e B u − d , 101 (1994), hep-ph/9812448.[2] X.-D. Ji, Phys. Rev. D55 , 7114 (1997), hep-ph/9609381.[3] A. V. Radyushkin, Phys. Rev.
D56 , 5524 (1997), hep-ph/9704207.[4] M. Diehl, Phys. Rept. , 41 (2003), hep-ph/0307382.[5] X.-D. Ji, Phys. Rev. Lett. , 610 (1997), hep-ph/9603249.[6] M. Burkardt, Phys. Rev. D62 , 071503 (2000), hep-ph/0005108.[7] S. Stepanyan et al. (CLAS), Phys. Rev. Lett. , 182002 (2001), hep-ex/0107043.[8] A. Airapetian et al. (HERMES), Phys. Rev. Lett. , 182001 (2001), hep-ex/0106068.[9] S. Chekanov et al. (ZEUS), Phys. Lett. B573 , 46 (2003), hep-ex/0305028.[10] A. Aktas et al. (H1), Eur. Phys. J.
C44 , 1 (2005), hep-ex/0505061.[11] M. Diehl, T. Feldmann, R. Jakob, and P. Kroll, Eur. Phys. J.
C39 , 1 (2005), hep-ph/0408173.[12] M. G¨ockeler et al. (QCDSF), Phys. Rev. Lett. , 042002 (2004), hep-ph/0304249.[13] P. H¨agler et al. (LHPC), Phys. Rev. D68 , 034505 (2003), hep-lat/0304018.[14] P. H¨agler et al. (LHPC), Phys. Rev. Lett. , 112001 (2004), hep-lat/0312014.[15] W. Schroers et al. (LHPC), Nucl. Phys. Proc. Suppl. , 907 (2004), hep-lat/0309065.[16] M. G¨ockeler et al. (QCDSF), Phys. Lett. B627 , 113 (2005), hep-lat/0507001.[17] M. Diehl et al. (QCDSF) (2005), hep-ph/0511032.[18] M. G¨ockeler et al., Phys. Rev. Lett. , 222001 (2007), hep-lat/0612032.[19] D. B. Renner et al. (LHP), Nucl. Phys. Proc. Suppl. , 255 (2005), hep-lat/0409130. -t[GeV ˜ e A u − d e B u − d e A u − d e B u − d e A u − d e B u − d ˜ e A u − d e B u − d e A u − d e B u − d e A u − d e B u − d LAT2005 , 056 (2005), hep-lat/0509185.[21] M. Diehl, Eur. Phys. J.
C19 , 485 (2001), hep-ph/0101335.[22] P. H¨agler, Phys. Lett.
B594 , 164 (2004), hep-ph/0404138.[23] D. Dolgov et al. (LHPC), Phys. Rev.
D66 , 034506 (2002), hep-lat/0201021.[24] M. G¨ockeler et al., Phys. Rev.
D54 , 5705 (1996), hep-lat/9602029.[25] K. Orginos and D. Toussaint (MILC), Phys. Rev.
D59 , 014501 (1999), hep-lat/9805009.[26] K. Orginos, D. Toussaint, and R. L. Sugar (MILC), Phys. Rev.
D60 , 054503 (1999), hep-lat/9903032.[27] S. R. Sharpe, PoS
LAT2006 , 022 (2006), hep-lat/0610094.[28] C. T. H. Davies et al. (HPQCD), Phys. Rev. Lett. , 022001 (2004), hep-lat/0304004.[29] C. W. Bernard et al., Phys. Rev. D64 , 054506 (2001), hep-lat/0104002.[30] A. Hasenfratz and F. Knechtli, Phys. Rev.
D64 , 034504 (2001), hep-lat/0103029.[31] C. Aubin et al., Phys. Rev.
D70 , 094505 (2004), hep-lat/0402030.[32] D. B. Kaplan, Phys. Lett.
B288 , 342 (1992), hep-lat/9206013.[33] Y. Shamir, Nucl. Phys.
B406 , 90 (1993), hep-lat/9303005.[34] R. Narayanan and H. Neuberger, Phys. Lett.
B302 , 62 (1993), hep-lat/9212019.[35] V. Furman and Y. Shamir, Nucl. Phys.
B439 , 54 (1995), hep-lat/9405004.[36] T. Blum et al., Phys. Rev.
D69 , 074502 (2004), hep-lat/0007038.[37] T. Blum, Nucl. Phys. Proc. Suppl. , 167 (1999), hep-lat/9810017.[38] T. Blum et al., Phys. Rev. D66 , 014504 (2002), hep-lat/0102005.[39] B. Bistrovic, Ph.D. thesis, MIT (2005).[40] B. Bistrovic, R. Erwin, J.W. Negele and M. Ramsey-Musolf, to be published.[41] K. Goeke, M. V. Polyakov, and M. Vanderhaeghen, Prog. Part. Nucl. Phys. , 401 (2001), hep-ph/0106012. -t[GeV ˜ e A u + d e B u + d e A u + d e B u + d e A u + d e B u + d ˜ e A u + d e B u + d e A u + d e B u + d e A u + d e B u + d B596 , 99 (2001), hep-ph/0009254.[43] M. Diehl, T. Feldmann, R. Jakob, and P. Kroll, Nucl. Phys.
B596 , 33 (2001), hep-ph/0009255.[44] S. R. Beane, P. F. Bedaque, K. Orginos, and M. J. Savage (NPLQCD), Phys. Rev.
D73 , 054503 (2006), hep-lat/0506013.[45] S. R. Beane et al., Phys. Rev.
D74 , 114503 (2006), hep-lat/0607036.[46] S. R. Beane, P. F. Bedaque, K. Orginos, and M. J. Savage, Phys. Rev.
D75 , 094501 (2007), hep-lat/0606023.[47] R. G. Edwards et al., PoS
LAT2006 , 121 (2006), hep-lat/0610007.[48] D.B. Renner et al., LHPC, to be published.[49] T. R. Hemmert, M. Procura, and W. Weise, Phys. Rev.
D68 , 075009 (2003), hep-lat/0303002.[50] S. R. Beane and M. J. Savage, Phys. Rev.
D70 , 074029 (2004), hep-ph/0404131.[51] K. Ackerstaff et al. (HERMES), Phys. Lett.
B464 , 123 (1999), hep-ex/9906035.[52] A. Airapetian et al. (HERMES), Phys. Rev.
D75 , 012007 (2007).[53] V. Y. Alexakhin et al. (COMPASS), Phys. Lett.
B647 , 8 (2007), hep-ex/0609038.[54] J.-W. Chen and X.-D. Ji, Phys. Lett.
B523 , 107 (2001), hep-ph/0105197.[55] J.-W. Chen and X.-D. Ji, Phys. Rev. Lett. , 052003 (2002), hep-ph/0111048.[56] D. Arndt and M. J. Savage, Nucl. Phys. A697 , 429 (2002), nucl-th/0105045.[57] W. Detmold, W. Melnitchouk, J. W. Negele, D. B. Renner, and A. W. Thomas, Phys. Rev. Lett. , 172001 (2001),hep-lat/0103006.[58] A. V. Belitsky and X. Ji, Phys. Lett. B538 , 289 (2002), hep-ph/0203276.[59] R. D. Young, D. B. Leinweber, and A. W. Thomas, Prog. Part. Nucl. Phys. , 399 (2003), hep-lat/0212031.[60] W. Detmold and C. J. D. Lin, Phys. Rev. D71 , 054510 (2005), hep-lat/0501007.[61] J.-W. Chen, D. O’Connell, R. S. Van de Water, and A. Walker-Loud, Phys. Rev.
D73 , 074510 (2006), hep-lat/0510024.[62] M. Diehl, A. Manashov, and A. Schafer, Eur. Phys. J.
A31 , 335 (2007), hep-ph/0611101. -t[GeV ˜ e A u + d e B u + d e A u + d e B u + d e A u + d e B u + d ˜ e A u + d e B u + d e A u + d e B u + d e A u + d e B u + d D74 , 094013 (2006), hep-ph/0602200.[64] M. Diehl, A. Manashov, and A. Sch¨afer, Eur. Phys. J.
A29 , 315 (2006), hep-ph/0608113.[65] M. Dorati, T. A. Gail, and T. R. Hemmert, Nucl. Phys.
A798 , 96 (2008), nucl-th/0703073.[66] P. Wang, D. B. Leinweber, A. W. Thomas, and R. D. Young, Phys. Rev.
D75 , 073012 (2007), hep-ph/0701082.[67] B. C. Tiburzi, Phys. Rev.
D72 , 094501 (2005), hep-lat/0508019.[68] O. Bar, C. Bernard, G. Rupak, and N. Shoresh, Phys. Rev.
D72 , 054502 (2005), hep-lat/0503009.[69] C. Aubin, J. Laiho, and R. S. Van de Water, Phys. Rev.
D75 , 034502 (2007), hep-lat/0609009.[70] J.-W. Chen, D. O’Connell, and A. Walker-Loud, Phys. Rev.
D75 , 054501 (2007), hep-lat/0611003.[71] M. Procura, B. U. Musch, T. Wollenweber, T. R. Hemmert, and W. Weise, Phys. Rev.
D73 , 114510 (2006), hep-lat/0603001.[72] C. Best et al., Phys. Rev.
D56 , 2743 (1997), hep-lat/9703014.[73] D. Br¨ommel et al., PoS
LAT2005 , 360 (2006), hep-lat/0509133.[74] S. Capitani et al., Phys. Lett.
B639 , 520 (2006), hep-lat/0511013.[75] P. J. Sutton, A. D. Martin, R. G. Roberts, and W. J. Stirling, Phys. Rev.
D45 , 2349 (1992).[76] M. Gl¨uck, E. Reya, and I. Schienbein, Eur. Phys. J.
C10 , 313 (1999), hep-ph/9903288.[77] K. Wijesooriya, P. E. Reimer, and R. J. Holt, Phys. Rev.
C72 , 065203 (2005), nucl-ex/0509012.[78] C. Alexandrou et al., Phys. Rev.