Lattice QCD evaluation of the Compton amplitude employing the Feynman-Hellmann theorem
K. U. Can, A. Hannaford-Gunn, R. Horsley, Y. Nakamura, H. Perlt, P. E. L. Rakow, G. Schierholz, K. Y. Somfleth, H. Stüben, R. D. Young, J. M. Zanotti
AADP-20-16/T1126, DESY 20-111, Liverpool LTH 1238
Lattice evaluation of the Compton amplitude employing the Feynman-Hellmanntheorem
K. U. Can, A. Hannaford-Gunn, R. Horsley, Y. Nakamura, H. Perlt, P. E. L. Rakow, G. Schierholz, K. Y. Somfleth, H. St¨uben, R. D. Young, and J. M. Zanotti (QCDSF/UKQCD/CSSM Collaborations) CSSM, Department of Physics, The University of Adelaide, Adelaide SA 5005, Australia School of Physics and Astronomy, University of Edinburgh, Edinburgh EH9 3JZ, UK RIKEN Center for Computational Science, Kobe, Hyogo 650-0047, Japan Institut f¨ur Theoretische Physik, Universit¨at Leipzig, 04103 Leipzig, Germany Theoretical Physics Division, Department of Mathematical Sciences,University of Liverpool, Liverpool L69 3BX, United Kingdom Deutsches Elektronen-Synchrotron DESY, 22603 Hamburg, Germany Regionales Rechenzentrum, Universit¨at Hamburg, 20146 Hamburg, Germany (Dated: July 6, 2020)The forward Compton amplitude describes the process of virtual photon scattering from a hadronand provides an essential ingredient for the understanding of hadron structure. As a physicalamplitude, the Compton tensor naturally includes all target mass corrections and higher twisteffects at a fixed virtuality, Q . By making use of the second-order Feynman-Hellmann theorem,the nucleon Compton tensor is calculated in lattice QCD at an unphysical quark mass across arange of photon momenta 3 ≲ Q ≲ . This allows for the Q dependence of the low momentsof the nucleon structure functions to be studied in a lattice calculation for the first time. Theresults demonstrate that a systematic investigation of power corrections and the approach to partonasymptotics is now within reach. Keywords: nucleon structure, parton distributions, Feynman Hellmann, compton amplitude, power correc-tions, scaling, lattice QCD
I. INTRODUCTION
Understanding the internal structure of hadrons fromfirst principles remains one of the foremost tasks in par-ticle and nuclear physics. It is an active field of researchwith important phenomenological implications in high-energy, nuclear and astroparticle physics. The staticproperties of hadrons, from the hybrid structure of quarkand meson degrees of freedom at low energies down tothe partonic structure at short distances, are encodedin structure functions. The tool for computing hadronstructure functions from first principles is lattice QCD.The connection between nucleon structure functionsand the quark structure of the nucleon is commonly ren-dered by the parton model. Although providing an in-tuitive language in which to interpret the deep-inelasticscattering data, the parton model represents an idealcase, valid if the partons are scattered elastically andincoherently by the incoming lepton. In the operatorproduct expansion (OPE) of the Compton amplitude theoperators are classified according to twist. The partonmodel accounts for twist-two contributions only, and doesnot accommodate power corrections arising from opera-tors of higher twist. It has been known for a long timethough that contributions from operators of higher twistare inseparably connected with the contributions of lead-ing twist, as a result of operator mixing and renormal-ization [1, 2].So far lattice calculations of nucleon structure func-tions have largely been limited to matrix elements of leading twist. That includes the calculation of a fewlower moments of parton distribution functions (PDFs),see Refs. [3, 4] for pioneering work and Ref. [5] for a sum-mary of the latest activity. More recently, there has beensignificant activity in the calculation of so-called quasi-PDFs [6–8] and its extensions, pseudo-PDFs [9] — seeRef. [10] for a review. Various other strategies have alsobeen proposed to overcome shortcomings of the study ofstructure functions on the lattice [11–16] and other re-lated inclusive processes [17–21].In the present work, we build upon a recent Letter[22] outlining a procedure to determine nucleon structurefunctions from a lattice calculation of the forward Comp-ton amplitude. By working with the physical amplitude,this approach overcomes issues of operator mixing andrenormalization, and the restriction to light-cone opera-tors [23, 24]. Here we establish the theoretical foundationof the approach and present results for the Compton am-plitude across a range of kinematics. The calculations areperformed at the SU(3) flavor symmetric point [25] at anunphysical pion mass. Results are reported on the low-est four moments of the unpolarized structure functionsof the nucleon for photon momenta Q ranging from ap-proximately 3–7 GeV . The variation of Q demonstratesthe potential to provide a quantitative test of the twistexpansion on the lattice for the first time.In terms of the practical computation, the determi-nation of the Compton amplitude takes advantage of theFeynman-Hellmann [26–30] approach to hadron structure— see also Refs. [31–35]. The use of Feynman-Hellmann a r X i v : . [ h e p - l a t ] J u l allows one to avoid the need to compute 3- or 4-pointfunctions and, thereby, obtain statistically cleaner sig-nals as well as maintain better control over excited-statecontamination. Here we also present a derivation of thesecond-order Feynman-Hellmann theorem necessary forthe present work — a related derivation has been pre-sented in Ref. [36].This paper is organized as follows: formal definitionsof the Compton amplitude and the structure functions,along with the connection between the OPE and the dis-persion relation are given in Section II. We explicitlyderive the second order Feynman-Hellmann theorem inSection III. Our lattice setup and the implementationdetails are given in Section IV. Results for the Comptonamplitude and the moments of the structure functionsare presented in Section V. We summarize our findingsin Section VI. II. FORWARD COMPTON AMPLITUDE ANDTHE STRUCTURE FUNCTIONSA. Notation
At leading order in the electromagnetic interaction,the general description for the inclusive scattering of acharged lepton from a hadronic target, e.g. eN → e ′ X ,is encoded in the hadron tensor. Conventionally, thehadron tensor is expressed as a matrix element of thecommutator of electromagnetic current operators [37–39], W µν ( p, q ) = π ∫ d z e iq ⋅ z ρ ss ′ ⟨ p, s ′ ∣[J µ ( z ) , J ν ( )]∣ p, s ⟩ , (1)for a hadron of momentum p and (virtual) photon mo-mentum q . For the present discussion, we will only con-sider spin-averaged observables by taking ρ ss ′ = δ ss ′ .The current operator takes the familiar form as thecharge-weighted sum of the quark vector currents, J µ =∑ f Q f J fµ , with Q f being the charge of quark flavor f .The flavor decomposition will be discussed in further de-tail in a later section.The spin-averaged nucleon tensor can be decomposedas W µν ( p, q ) = (− g µν + q µ q ν q ) F ( x, Q )+ ( p µ − p ⋅ qq q µ ) ( p ν − p ⋅ qq q ν ) F ( x, Q ) p ⋅ q , (2)which is defined such that Lorentz-invariant structurefunctions, F , , match onto their conventional partonicinterpretation in the deep inelastic scaling region. These In this section, we work in Minkowski space. structure functions are expressed as functions of theBjorken scaling variable ( x = Q /( p ⋅ q ) ) and Q = − q .While the inelastic structure functions are not directlyaccessible within a conventional Euclidean lattice for-mulation, we highlight that the spacelike component ofthe Compton tensor can be studied within a Euclideanframework—as also discussed in Refs. [11, 14, 17].The (spin-averaged) Compton tensor is defined simi-larly to Equation (1), T µν ( p, q ) = i ∫ d z e iq ⋅ z ρ ss ′ ⟨ p, s ′ ∣T {J µ ( z )J ν ( )}∣ p, s ⟩ , (3)where T is the time-ordering operator. This tensor canbe decomposed in precisely the same way as W µν inEquation (2), which defines the analogous scalar func-tions F , ( ω, Q ) . For our purposes, it is convenient toexpress these in terms of the inverse Bjorken variable ω = p ⋅ q / Q . These Compton structure functions are re-lated to the corresponding ordinary structure functionsvia the optical theorem, which states:Im F ( ω, Q ) = πF ( x, Q ) , (4)Im F ( ω, Q ) = πF ( x, Q ) . (5)Analyticity and crossing symmetry means we can writea dispersion relation for F [40] F ( ω, Q ) − F ( , Q ) = ω π ∫ ∞ dω ′ Im F ( ω ′ , Q ) ω ′ ( ω ′ − ω − i(cid:15) ) , (6) F ( ω, Q ) = ωπ ∫ ∞ dω ′ Im F ( ω ′ , Q ) ω ′ − ω − i(cid:15) . (7)To accommodate the subtraction necessary in F , we willmake use of the bar notation to denote the dispersivepart, F ( ω, Q ) = F ( ω, Q ) − F ( , Q ) . The dispersionintegrals can be directly connected to the hadron tensorby the optical theorem, giving: F ( ω, Q ) = ω ∫ dx x F ( x, Q ) − x ω − i(cid:15) , (8) F ( ω, Q ) = ω ∫ dx F ( x, Q ) − x ω − i(cid:15) . (9)The nature of the dispersion integral makes it clearthat whenever ∣ ω ∣ < i(cid:15) becomes irrelevant.Hence the current-current correlation remains spacelikeand there is no distinction between the Euclidean andMinkowski amplitudes. Physically, the condition ∣ ω ∣ < B. Operators displaced in time
In the Feynman-Hellmann approach employed in thiswork, the matrix elements calculated involve current–current correlations which are displaced in Euclideantime. It is therefore instructive to further clarify therelationship between the Compton tensor, as definedin Minkowski space, and the corresponding calculationwithin a Euclidean framework.We start by separating Equation (3) into two distincttime orderings by first defining the amplitude at fixedtemporal separation between the currents:˜ T M µν ( p, q, t ) = iρ ss ′ ∫ d z e i ( q + i(cid:15) ) t e − i q ⋅ z × ⟨ p, s ′ ∣J µ ( z , t )J ν ( )∣ p, s ⟩ . (10)From this definition it is straightforward to recover thefull Compton amplitude by integrating over t : T µν ( p, q ) = ∫ ∞ dt [ ˜ T M µν ( p, q, t ) + ˜ T M νµ ( p, − q, t )] . (11)To isolate the explicit t dependence in Equation (10),we insert a complete set of states and exploit transla-tional invariance in the usual way:˜ T M µν ( p, q, t ) = iρ ss ′ ⨋ X ∫ d z e i ( q + E p − E X + i(cid:15) ) t × e − i ( q + p − P X )⋅ z ⟨ p, s ′ ∣J µ ( )∣ X ⟩⟨ X ∣J ν ( )∣ p, s ⟩ . (12)The completeness integral, I = ⨋ X ∣ X ⟩⟨ X ∣ , describes a fullintegral over the entire state space, implicitly includingall possible momenta over all possible configurations ofparticles.Similarly to Equation (10), one can write down an ex-pression where the current insertions are separated inEuclidean time [11, 39]:˜ T E µν ( p, q, τ ) = ρ ss ′ ∫ d z e q τ e − i q ⋅ z × ⟨ p, s ′ ∣J µ ( z , τ )J ν ( )∣ p, s ⟩ . (13)Within the Feynman-Hellmann approach, by construc-tion, there is no external energy transfer, q = T E µν ( p, q, τ ) = ρ ss ′ ⨋ X ∫ d z e ( q + E p − E X ) τ × e − i ( q + p − P X )⋅ z ⟨ p, s ′ ∣J µ ( )∣ X ⟩⟨ X ∣J ν ( )∣ p, s ⟩ . (14)It is evident that Equations (12) and (14) differ non-trivially in their dependence on the temporal coordinate— a similar point has been made in Ref. [41]. However,upon integrating with respect to time, the Minkowskiand Euclidean expressions are easily equated — subjectto the caveat that we are below the elastic threshold.In particular, provided that E X ( p ± q ) > E p ± q for allnon-vanishing contributions to the completeness sum ⨋ X ,then the i(cid:15) prescription becomes irrelevant and we have: ∫ ∞ dτ ˜ T E µν ( p, ± q, τ ) = ∫ ∞ dt ˜ T M µν ( p, ± q, t ) . (15) By summing the two terms of Equation (11), this makesit clear that the Euclidean and Minkowski Comptonamplitudes are identical in the unphysical region, eventhough the current insertions are allowed to be separatedin time. We of course note that if the intermediate statescan go on shell, E X ( p ± q ) = E p ± q , then the Euclideanintegral in Equation (15) is not well defined [42]. The i(cid:15) factor in Equation (10) then becomes essential in orderto define the analytic continuation required to render theintegral finite. However, the fixed- t Euclidean matrix el-ements are perfectly well defined, and there is no restric-tion on the kinematic thresholds—as has been studiedin [11, 39, 42].Although we have just been through this carefulconsideration of the t dependence, we note that theFeynman-Hellmann technique relies on resolving thespectrum of a perturbed Hamiltonian—which itself doesnot make any reference to the nature of the temporalcorrelations. In the present work, there is hence no needto explicitly resolve the τ dependence of Equation (13).The connection to the Minkowski amplitude lies in Equa-tion (15), as the current insertions act at all times. C. Moments and the OPE
As will become clear in the next section, within theFeynman-Hellmann formalism, each external current mo-mentum vector q of interest requires a unique propagatorinversion. However, for each inversion, the variation ofhadron Fourier momenta p allows access to multiple dis-tinct ω = p . q / Q values. An ensemble of q and p valuestherefore provides a wealth of kinematic points to re-solve the Compton amplitude. It is therefore convenientto present a summary of the kinematic coverage in termsof the moments of the structure functions. We considerthe Taylor series expansion of Equation (8) at fixed Q : F ( ω, Q ) = ∞ ∑ n = ω n M ( ) n ( Q ) , (16) F ( ω, Q ) = ∞ ∑ n = ω n − M ( ) n ( Q ) , (17)with the moments being defined by M ( ) n ( Q ) = ∫ dx x n − F ( x, Q ) , (18) M ( ) n ( Q ) = ∫ dx x n − F ( x, Q ) . (19)From the perspective of the present lattice calculation,one can proceed by calculating the Compton tensor fora number of ω values and extract the moments of thestructure functions. By treating Q as an external scale,the method connects directly to the physical amplitudesof interest, and therefore circumvents the operator mix-ing issues discussed above. Although the described lat-tice calculations relate directly to the physical moments,we note that at asymptotically large Q the momentsbecome dominated by their leading-twist contributions,namely the moments of the familiar parton distributionfunctions, v n , M ( ) n ( Q ) = ∑ f C ( ) f, n ( Q µ , g ( µ )) v f n ( µ ) + O ( Q ) , (20) M ( ) n ( Q ) = ∑ f C ( ) f, n ( Q µ , g ( µ )) v f n ( µ ) + O ( Q ) , (21)where the sum runs over partonic flavors f . The shortdistance structure of the operator product in Equa-tion (3) is encoded in the Wilson coefficients, C , and thelong-distance hadronic features are encoded in the ma-trix elements of local operators, v , renormalized at somescale µ . For completeness, our notation is summarized inAppendix A. III. SECOND ORDER FEYNMAN-HELLMANNTHEOREM
The purpose of applying the Feynman-Hellmann theo-rem to lattice QCD is to relate matrix elements of interestto energy shifts in weak external fields. In the case of ageneralized Compton amplitude, described by a matrixelement of two (non-local) current insertions, the conven-tional approach would require the evaluation of lattice 4-point functions. The application of Feynman-Hellmannthen reduces the problem to a more straightforward anal-ysis of 2-point correlation functions using spectroscopictechniques.In order to compute the forward Compton amplitudevia the Feynman-Hellmann relation, we introduce the fol-lowing perturbation to the fermion action, S ( λ ) = S + λ ∫ d z ( e i q ⋅ z + e − i q ⋅ z )J µ ( z ) , (22)where λ is the strength of the coupling between thequarks and the external field, J µ ( x ) = Z V ¯ q ( x ) γ µ q ( x ) isthe electromagnetic current coupling to the quarks alongthe µ direction, q is the external momentum inserted bythe current and Z V is the renormalization constant forthe local electromagnetic current.The general strategy for deriving Feynman-Hellmannin a lattice QCD context is to consider the general spec-tral decomposition of a correlator in the presence of thebackground field. The differentiation of this correlationfunction with respect to the external field reveals a dis-tinct temporal signature for the energy shift. By explicitevaluation of the perturbed correlator, one is able to iden-tify this signature and hence resolve the desired relation-ship between the energy shift and matrix element. Ourprinciple theoretical result here is that for the perturbed action described in Equation (22), the second-order en-ergy shift of the nucleon is found to be: ∂ E N λ ( p ) ∂λ ∣ λ = = − T µµ ( p, q ) + T µµ ( p, − q ) E N ( p ) , (23)where T is the Compton amplitude defined in Equa-tion (3), q = ( q , ) is the external momentum encodedby Equation (22), and E N λ ( p ) is the nucleon energy atmomentum p in the presence of a background field ofstrength λ . In the following we sketch the main stepsof the derivation, and refer the interested reader to Ap-pendix B for further details.In the presence of the external field introduced inEquation (22), we define the two-point correlation func-tion projected to definite momentum as, G ( ) λ ( p ; t ) ≡ ∫ d xe − i p ⋅ x Γ ⟨ Ω λ ∣ χ ( x , t ) χ ( )∣ Ω λ ⟩ , (24)where here and in the following, a trace over Dirac indiceswith the spin-parity projection matrix Γ is understood,and ∣ Ω λ ⟩ is the vacuum in the presence of the externalfield. The asymptotic behavior of the correlator at largeEuclidean times takes the familiar form, G ( ) λ ( p ; t ) ≃ A λ ( p ) e − E Nλ ( p ) t , (25)where E N λ ( p ) is the energy of the ground state nucleonin the external field and A λ ( p ) the corresponding overlapfactor.For the purpose of current presentation, a nucleon in-terpolating operator is assumed for χ . However, thederivation applies to any ground-state hadron, providedthe ground state in the presence of the external fieldis perturbatively close to the free-field state. A simplecounter example could be a Σ baryon in the presenceof a strangeness-changing current, where at λ = e − E Σ t but at any finite λ this willeventually be dominated by e − E N t (kinematics permit-ting).It is for a similar physical reason that one must workwith nucleon states that have the least possible kineticenergy among all states connected to any number of cur-rent insertions. This same condition guarantees the con-nection between the Euclidean and Minkowski Comptonamplitudes described in the previous section. In the pres-ence of the background field, the Hamiltonian of the sys-tem will mix momentum states p ± n q . We hence choosethe Fourier projection of our correlation function, Equa-tion (24), such that p corresponds to the lowest kineticenergy, ∣ p ∣ < ∣ p + n q ∣ ∀ n ∈ Z . For notational purposes,at non-zero λ , the label p can be understood to describethe space of all momentum states connected by an integermultiple of q .Assuming that first-order perturbations of the energyvanish, as ensured by the chosen kinematics, the second-order derivative of Equation (24), evaluated at λ = ∂ G ( ) λ ( p ; t ) ∂λ ∣ λ = = ( ∂ A λ ( p ) ∂λ − tA ( p ) ∂ E N λ ( p ) ∂λ )× e − E N ( p ) t . (26)The derivatives of A λ ( p ) and E N λ ( p ) are assumed to beevaluated at λ =
0. The first term corresponds to theshift in the overlap factor and the second order energyshift is identified in the t -enhanced (or time-enhanced)term. It is this t enhancement that leads to a relation-ship between the energy shift and matrix element. Hence to complete the derivation, we differentiate the path inte-gral representation directly to identify the time-enhancedcontributions to the correlator.The path integral expression for the 2-point correlator,Equation (24), in the background field is given by G ( ) λ ( p ; t ) = ∫ d x e − i p ⋅ x Γ λ ⟨ χ ( x , t ) χ ( )⟩ λ , (27)where λ ⟨⋯⟩ λ denotes the full path integral over all fields,using the perturbed action S ( λ ) given in Equation (22)—an absence of λ subscript is taken to imply λ →
0. Bydifferentiating twice with respect to λ and evaluated at λ →
0, one finds (see Appendix B for details) ∂ G ( ) λ ( p ; y ) ∂λ RRRRRRRRRRR λ = = ∫ d x e − i p ⋅ x Γ ⎡⎢⎢⎢⎢⎣⟨ χ ( x , t ) χ ( ) ( ∂S ( λ ) ∂λ ) ⟩ + ⟨ χ ( x , t ) χ ( )⟩⟨( ∂S ( λ ) ∂λ ) ⟩⎤⎥⎥⎥⎥⎦ . (28)To arrive at this form, it is assumed that the vacuum expectation value of a single current insertion vanishes, ⟨ ∂S ( λ )/ ∂λ ⟩ =
0, such as is the case for the electromagnetic current. It is clear that the second term in Equation (28)only acts to modify the unperturbed correlator, and hence cannot generate the temporal enhancement associated withthe energy shift. Focusing purely on the first term, and inserting an explicit form for the electromagnetic externalfield, the corresponding second derivative of the correlator becomes ∂ G ( ) λ ( p ; t ) ∂λ RRRRRRRRRRR λ = = ∫ d x e − i p ⋅ x Γ ∫ d yd z ( e i q ⋅ y + e − i q ⋅ y )( e i q ⋅ z + e − i q ⋅ z )⟨ χ ( x , t )J µ ( z )J µ ( y ) χ ( )⟩ . (29)The correlator defined here involves a four-point correlation function with nucleon interpolating operators held atfixed temporal separation t , with the currents inserted across the entire four-volume. Importantly, this expressionis evaluated in the absence of the external field, and hence momentum conservation is exact. It is then possible toperform a spectral decomposition of this correlator in terms of a transfer matrix that is diagonal in the momenta.It is a rather straightforward calculation to perform the standard procedure of inserting a complete sets of states,and then exploit translational invariance to complete the spatial integrals. Since the temporal integrals over thecurrents extend over all time, each distinct time ordering of the 4-point function must be treated separately. Howeverthe contribution to the energy shift can only come from the contribution where the two current operators both appearbetween the nucleon creation and annihilation operators. Isolating the contributions that give rise to the dominant te − E N ( p ) t behavior at asymptotic times gives: ∂ G ( ) λ ( p ; t ) ∂λ RRRRRRRRRRR λ = = tA ( p ) e − E N ( p ) t E N ( p ) ⟨ N ( p ) ∣∫ d z ( e iq ⋅ z + e − iq ⋅ z ) J µ ( z )J µ ( )∣ N ( p )⟩ + . . . , (30)where the subleading terms are suppressed by the el-lipsis. Note that the spin indices have been suppressedhere, however a detailed presentation is provided in Ap-pendix B. Finally, a comparison of this form with Equa-tion (26) and the Compton amplitude, Equation (3),leads to the result quoted in Equation (23).In principle the derivation presented in this section(and Appendix B) can be generalized to mixed cur-rents by adding an additional perturbation (or current)to Equation (22) with a different coupling strength, λ ′ ,and current momentum, q ′ , which can, in general, betaken to be different from λ and q . This would allow access to interference terms, allowing one to study u – d flavor interference effects, spin-dependent amplitudes orthe off-forward Compton amplitude and generalized par-ton distributions with q ≠ q ′ . Details of a prescriptionfor the off-forward Compton tensor will be presented ina forthcoming paper [43]. Table I. Details of the gauge ensembles used in this work. N f c SW κ l κ s L × T a m π m N m π L Z V N cfg [fm] [GeV] [GeV]2 + ×
64 0.074(2) 0 . ( ) . ( ) . IV. SIMULATION DETAILSA. Gauge ensembles
We use a single gauge ensemble generated by theQCDSF/UKQCD Collaborations employing a stout-smeared non-perturbatively O( a ) -improved Wilson ac-tion for the dynamical up/down and strange quarks anda tree-level Symanzik improved gauge action [44]. Wework on a volume of L × T = ×
64, the bare cou-pling parameter is β = .
5, and the lattice spacing, a = . ( ) fm, is set using a number of flavor-singletquantities [45–48]. We are working on the SU ( ) -flavorsymmetric point where the masses of all three quark fla-vors are set to approximately the physical flavor-singletmass, m = ( m s + m l )/ ≃
470 MeV and m π L = .
6. Further details and advan-tages of this choice are discussed in [25, 46]. The renor-malization constant for the local vector current is deter-mined to be Z V = . ( ) by imposing the chargeconservation on the Sachs electric form factor calculatedat Q =
0. This value is in agreement within statisticalprecision with the value determined in the chiral limitusing the RI ′ -MOM scheme [49]. We tabulate the detailsin Table I for the Reader’s convenience. B. Feynman-Hellmann implementation
We implement the second-order Feynman-Hellmanntheorem through the valance quarks. It can clearly beimplemented at the hybrid Monte Carlo level, whichwould relay the effects of the perturbations to the sea-quarks [29], however this would lead to a significant in-crease in required computing resources, so in this workwe focus on the quark-line connected contributions to theCompton amplitude. To this end, we add the perturba-tion given in Equation (22), S ( λ ) = S + λ ∫ d z ( e i q ⋅ z + e − i q ⋅ z )J ( z ) , (31)to the valence quark action only, where the renormalizedlocal vector current, J ( x ) = Z V ¯ q ( x ) iγ q ( x ) , is chosento be along the z -direction, µ =
3. The second exponen-tial term symmetrizes the Fourier transform and ensuresthe Hermiticity of the action. In order to evaluate thesecond-order energy shift with respect to λ at λ =
0, onehas to compute additional quark propagators at severalchoices of λ . This added cost of computation is coun-tered by optimizing the inversion of the perturbed Dirac X X X XX X XXXXXX uud
Figure 1. The six possible ways of inserting two currents toa nucleon. Upper three correspond to the uu flavor contribu-tions while the lower leftmost is for dd . Remaining two arefor the ud , which we omit in this work. matrix. We adopt an approach where we feed the unper-turbed propagator as an initial guess to the inversion ofthe perturbed one, which results in roughly a factor of10 gain in inversion time.In order to improve the stability of estimating theenergy-shifts at λ =
0, one aims to introduce the small-est possible perturbations by choosing a suitable λ . Theobjective is to keep λ sufficiently small to minimize thecontamination from λ effects, yet large enough to ensurethat the perturbation is not lost within the numericalprecision of the calculation. Our tests indicate that anychoice in the range 10 − > ∣ λ ∣ ≥ − leads to meaningfulresults. Note that the upper bound is sensitive to thequark mass, where a too large λ might lead to increasedinstabilities in the Dirac matrix inversion, particularly asone approaches the physical point. C. Flavor decomposition
The implementation of the Feynman-Hellmann theo-rem described above effectively inserts an external cur-rent on to a quark line by computing its propagator withthe perturbed quark action Equation (31). When bothcurrents are inserted onto the u -quarks or the d -quark, weevaluate the “ uu ” or “ dd ” contributions to the Comptonstructure functions, respectively. By employing positiveand negative pairs of λ ’s (see Section V), one can form u + d and u − d type insertions leading to the possibilityfor isolating a “ ud ” insertion where one current hits a u -quark and the other the d quark. The six different waysof inserting the currents are shown in Figure 1. The ud contribution is particularly interesting since it directlycorresponds to a higher-twist contribution [50], i.e. thetwist-4 cat’s ears diagram. An investigation of these con-tributions is left for future work. D. Isolating the energy shift
The energy of the ground state, N λ , in a weakly cou-pled external field can be expanded as a Taylor series in λ , E N λ ( p ) = E N ( p ) + λ ∂E N λ ( p ) ∂λ ∣ λ = + λ ∂ E N λ ( p ) ∂λ ∣ λ = + O( λ ) . (32)Collecting terms that are even and odd in λ to all orders,we may rewrite the expansion as, E N λ ( p ) = E N ( p ) + ∆ E eN λ ( p ) + ∆ E oN λ ( p ) , (33)where E N ( p ) in the above two expressions correspondsto the unperturbed ( λ =
0) energy. In order to extractthe second order energy shift from the lattice correlationfunctions, we construct a ratio which isolates the even- λ energy shift, ∆ E eN λ ( p ) , R eλ ( p , t ) ≡ G ( )+ λ ( p , t ) G ( )− λ ( p , t )( G ( ) ( p , t )) (34) t ≫ ——→ A λ ( p ) e − E eNλ ( p ) t , (35)where the perturbed two-point functions, G ( )± λ ( p , t ) , aredefined in Equations (24) and (25) and G ( ) ( p , t ) isthe unperturbed one. The large t behavior given inEquation (35) is arrived at by combining Equations (25)and (33). While not necessary for our discussion, forcompleteness we note that the overlap factor is A λ ( p ) = ∣⟨ Ω λ ∣ χ ( )∣ N λ ( p )⟩∣ E N λ ( p ) ×∣⟨ Ω λ ∣ χ ( )∣ N − λ ( p )⟩∣ E N − λ ( p ) ( ∣⟨ Ω ∣ χ ( )∣ N ( p )⟩∣ E N ( p ) ) − . (36)Extraction of the even- λ energy shift ∆ E eN λ then followsstandard spectroscopy methods by fitting R eλ ( p , t ) de-fined in Equation (34) with a single exponential at suffi-ciently large times. Details follow in the next section. V. RESULTS AND DISCUSSIONA. Extracting the structure functions and themoments
In order to illustrate the feasibility and the versatilityof the method, we carry out simulations with several val-ues of current momentum, Q , in the range 3 ≲ Q ≲ . Utilizing up to six randomly placed quark sourcesper configuration, we perform up to O( ) measure-ments for each pair of λ and q . Quark fields are smearedin a gauge-invariant manner by Jacobi smearing [51], Table II. Multiple ω values that we can access with severalcombinations of p = ( p x , p y , p z ) and q = ( q x , q y , q z ) in latticeunits. Note that the ω ≥ values are omitted in the analysis—values of ω outside the allowed range are indicated by italics. ω = p ⋅ q / Q p /( π / L ) q /( π / L )( , , ) ( , , ) ( , , ) ( , , ) ( , , )( , , ) .
00 0 .
00 0 .
00 0 .
00 0 . ( , , ) .
20 0 .
31 0 .
12 0 .
20 0 . ( , , ) .
40 0 .
62 0 .
24 0 .
40 0 . ( , − , ) .
20 0 .
16 0 .
24 0 .
00 0 . ( , − , ) .
40 0 .
16 0 .
35 0 .
20 0 . ( , , ) .
60 0 .
46 0 .
47 0 .
40 0 . ( , , ) .
80 0 .
77 0 .
59 0 .
60 0 . ( , , ) .
71 0 .
80 0 . ( , − , ) .
62 0 .
82 0 .
60 0 . ( , , ) .
92 0 .
94 0 .
80 0 . ( , , ) . where the smearing parameters are tuned to produce arms radius of ≃ . ω values at eachsimulated value of q by varying the nucleon momentum p as shown in Table II.Setting p z = q z =
0, along with our choice of µ =
3, simplifies the Compton tensor such that the secondorder energy shift in Equation (23) corresponds to the F Compton structure function directly, ∂ E N λ ( p ) ∂ λ ∣ λ = = − T ( p, q ) + T ( p, − q ) E N ( p )= − F ( ω, Q ) E N ( p ) , (37)where the energy of the nucleon is calculated via the con-tinuum dispersion relation, E N ( p ) = √ m N + p .We compute the perturbed two-point correlationfunctions (Equation (24)) with four values of λ =[± . , ± . ] . Even- λ energy shifts are extractedfrom the ratio of correlation functions given in Equa-tion (34) following a covariance-matrix based χ analy-sis to pick the best available range for each ratio. Weshow a representative case for q = ( , , ) ( πL ) and p = ( , , ) ( πL ) in Figure 2. Having even- λ energy shiftsfor two λ values, we perform polynomial fits of the form,∆ E eN λ ( p ) = λ ∂ E N λ ( p ) ∂λ ∣ λ = + O( λ ) , (38)to determine the second order energy shift. Given thesmallness of our λ values, higher order O( λ ) terms areheavily suppressed, hence the fit form reduces to a simple t/a − . . . . . . . . l og R e λ ( p , t + ) R e λ ( p , t ) q =(4 , ,
0) (2 π/L ) p =(1 , ,
0) (2 π/L ) λ = 0 . λ = 0 . Figure 2. Effective mass plot of the ratio given in Equa-tion (34) for q = ( , , ) ( πL ) and p = ( , , ) ( πL ) . Shadedhorizontal regions indicate the fit windows along with the ex-tracted λ values with their 1 σ error margins. Data points areshifted for clarity. λ × − ∆ E e N λ ( p ) × − p =(0, 0, 0) b = 2.764 ± .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 . λ . . . . . ∆ E e N λ ( p ) × − q =(4, 1, 0) 2 π/Lp =(0, 0, 0) 2 π/L λ × − ∆ E e N λ ( p ) × − p =(0, 1, 0) b = 2.727 ± .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 . λ . . . . . ∆ E e N λ ( p ) × − q =(4, 1, 0) 2 π/Lp =(0, 1, 0) 2 π/L λ × − ∆ E e N λ ( p ) × − p =(0, 2, 0) b = 2.985 ± .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 . λ . . . . . ∆ E e N λ ( p ) × − q =(4, 1, 0) 2 π/Lp =(0, 2, 0) 2 π/L λ × − ∆ E e N λ ( p ) × − p =(1, -2, 0) b = 2.544 ± .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 . λ . . . . . ∆ E e N λ ( p ) × − q =(4, 1, 0) 2 π/Lp =(1, -2, 0) 2 π/L λ × − ∆ E e N λ ( p ) × − p =(1, -1, 0) b = 2.503 ± .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 . λ . . . . . ∆ E e N λ ( p ) × − q =(4, 1, 0) 2 π/Lp =(1, -1, 0) 2 π/L λ × − ∆ E e N λ ( p ) × − p =(1, 0, 0) b = 2.560 ± .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 . λ . . . . . ∆ E e N λ ( p ) × − q =(4, 1, 0) 2 π/Lp =(1, 0, 0) 2 π/L λ × − ∆ E e N λ ( p ) × − p =(1, 1, 0) b = 2.326 ± .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 . λ . . . . . ∆ E e N λ ( p ) × − q =(4, 1, 0) 2 π/Lp =(1, 1, 0) 2 π/L λ × − . . . ∆ E e N λ ( p ) × − p =(1, 2, 0) b = 1.459 ± .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 . λ . . . . . ∆ E e N λ ( p ) × − q =(4, 1, 0) 2 π/Lp =(1, 2, 0) 2 π/L λ × − ∆ E e N λ ( p ) × − p =(2, -1, 0) b = 1.596 ± .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 . λ . . . . . ∆ E e N λ ( p ) × − q =(4, 1, 0) 2 π/Lp =(2, -1, 0) 2 π/L λ × − ∆ E e N λ ( p ) × − p =(2, 0, 0) b = 0.907 ± .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 . λ . . . . . ∆ E e N λ ( p ) × − q =(4, 1, 0) 2 π/Lp =(2, 0, 0) 2 π/L λ × − − . . . ∆ E e N λ ( p ) × − p =(2, 1, 0) b = -0.014 ± .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 . λ . . . . . ∆ E e N λ ( p ) × − q =(4, 1, 0) 2 π/Lp =(2, 1, 0) 2 π/L Figure 3. λ dependence of ∆ E eN λ ( p ) given in Equation (38).Fit form is f ( λ ) = bλ . Error bars are smaller than the sym-bols. one parameter polynomial. A representative fit of Equa-tion (38) to the energy shifts for q = ( , , ) ( πL ) and p = ( , , ) ( πL ) is shown in Figure 3. In Figure 4 weshow the subtracted Compton structure function for uu and dd insertions as obtained from the energy shifts viaEquation (37) for each value of ω for q = ( , , ) ( πL ) .With our particular choice of Lorentz indices and mo-menta, we can connect the lattice Compton amplitude(Equation (37)) to the moments of the structure func-tions (Equation (16)) as, F ( ω, Q ) = ( ω M ( ) ( Q ) + ω M ( ) ( Q )+ ⋅ ⋅ ⋅ + ω n M ( ) n ( Q ) + ⋅ ⋅ ⋅) . (39)Since the Compton amplitude is directly related to theexperimental cross-section, it must be positive definitivefor the entire kinematic region. Consequently, this holdsfor the uu and dd moments as well. Hence, the moments M ( ) n are constrained to be monotonically decreasing, M ( ) ( Q ) ≥ M ( ) ( Q ) ≥ ⋅ ⋅ ⋅ ≥ M ( ) n ( Q ) ≥ ⋅ ⋅ ⋅ ≥ . (40)More generally, the sequence of moments satisfy theHausdorff moment criteria [52], yet the simple monotonicdecreasing form of Equation (40) allows an assessmentof the constraint on the moments provided by the data. .
12 0 .
24 0 .
35 0 .
47 0 .
59 0 .
71 0 .
82 0 . ω − . . . . . . F qq ( ω , Q ) q =(4 , ,
0) 2 π/L uudd mean ± stat. ± sys.µ uu . +0 . . − . − . µ uu . +0 . . − . − . µ uu . +0 . . − . − . µ uu . +0 . . − . − . µ dd . +0 . . − . − . µ dd . +0 . . − . − . µ dd . +0 . . − . − . µ dd . +0 . . − . − . Figure 4. ω dependence of the subtracted Compton structurefunction F qq ( ω, Q ) of nucleon at Q = .
66 GeV . Fits de-picting the extraction of the moments via Equation (39) with n = This series is rapidly converging and stable with respectto the truncation order. However, imposing the abovecondition in a least-squares analysis is not so straight-forward, but it can be easily implemented in a Bayesianapproach. The particular Bayesian inference implemen-tation we are employing [53] has the advantage of us-ing adaptive algorithms to optimize the hybrid MonteCarlo parameters [54], which removes the extra effort offine-tuning such parameters and returns the sampling re-sults within mere minutes with convergence checks imple-mented.In the present analysis, we sample the moments fromuniform distributions with bounds M ( ) ( Q ) ∈ [ , ] and M ( ) n ( Q ) ∈ [ , M ( ) n − ( Q )] , for n >
1, to enforcethe monotonic decreasing nature of the moments. Uni-form prior distributions are chosen since these are un-informative distributions and remove a source of bias.The sequences of individual uu or dd moments are se-lected according to a multivariate probability distribu-tion, exp (− χ / ) , where, χ = ∑ i,j [F ,i − F obs1 ( ω i )] C − ij [F ,j − F obs1 ( ω j )] , (41)is the χ function with the covariance matrix C ij , ensur-ing the correlations between the data points are takeninto account. Fits depicting the extraction of the mo-ments are shown in Figure 4. Note that, Equation (40)is not necessarily true for the isovector, uu − dd , mo-ments. Therefore, the Bayesian priors for the uu and dd are treated independently. However, by sampling the uu and the dd datasets within the same trajectory, we en-sure underlying correlations between those datasets areaccounted for. Hence, the indices i , j in Equation (41)run through all the ω values and both flavors.The first few moments extracted are given in Table IIIand the isovector moments are plotted in Figure 5 for Table III. First few moments of the structure function F forseveral values of Q . Contributions of each flavor are givenalong with the isovector quantity. Errors are statistical un-certainties only. q /( π / L ) Q [GeV ] M ( ) n ( Q ) uu dd uu − dd ( , , ) . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . ( , , ) . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . ( , , ) . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . ( , , ) . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . ( , , ) . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . M ( ) . + . − . . + . − . . + . − . n . . . . . . . . M ( ) n , uu − dd ( Q ) Q = 2 . Q = 3 . Q = 4 . Q = 5 . Q = 7 . Figure 5. Isovector moments given in Table III. Q are givenin GeV . each choice of Q . Error margins correspond to the high-est posterior density interval with a 68% credible regionand the asymmetric intervals reflect the shape of the pos-terior distributions. We find that the lower momentshave a negligible dependence on the truncation order ofthe series in Equation (39) for n ≥ u and the d moments, thisleads to rather large second moments for the isovector u − d combination, which are even comparable to that ofthe first moment in some cases. While this is likely due to the limited statistics of the current simulations, it mayalso be a signal of significant power corrections which wediscuss in the next section. B. Power corrections and scaling
The moments M ( ) n ( Q ) in Table III appear to indi-cate that power corrections are present throughout ourrange of photon momenta, 3 ≲ Q ≲ . The lead-ing power corrections to moments of structure functionshave essentially two sources, target mass (together withpossible threshold effects) and mixing with operators ofhigher twist. Target mass and threshold effects can be ac-counted for, to a certain extent, by replacing the Bjorken x scaling variable by a generalized scaling variable, ξ (e.g.[55–57]). For example, a commonly used form proposedby Nachtmann [56] is ξ = x + √ + m N x / Q , (42)where m N is the nucleon mass. Besides being power cor-rected, the various moments defined in terms of thesenew ξ scaling variables are mixtures of the (Cornwall-Norton [58]) moments defined in terms of x , with themixings suppressed by powers of 1 / Q . For a generalizedscaling variable incorporating the analytic structure ofthe forward Compton amplitude see Ref. [55].The second source of power corrections in the structurefunction moments are due to contributions from opera-tors with twist-4 and above, represented by the O( / Q ) terms in Equation (20). We note that while in princi-ple it is possible to compute the Wilson coefficient forthe twist-4 contribution to M ( ) ( Q ) non-perturbativelyon the lattice [13, 59, 60], the hyper-cubic nature of thelattice can only accommodate operators of spin four orless, which thwarts any direct prediction of the Wilsoncoefficients for higher moments.Since the moments computed in this work are deter-mined from a fit to the full Compton amplitude, theynaturally include all possible power corrections. Whilethe present data do not have the precision or range of Q to isolate individual power corrections, we are able toaccount for the observed Q -dependence of each momentby fitting with the functional form M ( ) n ( Q ) = M ( ) n + C n / Q + O( / Q ) . (43)Of particular interest is a fit to the lowest isovector mo-ment M ( ) ,uu − dd ( Q ) , which we show in Figure 6 as a func-tion of Q . Here we clearly see that the current data iswell described by Equation (43), with the fit form sug-gesting large power corrections may be present at low Q .However we add a cautionary note that the behavior atsmall Q is heavily influenced by the presence of the largevalue for the moment at Q = .
74 GeV and neglectingthis point would result in a much softer, although stillnon-trivial, Q -dependence in the small- Q region.0 Q [ GeV ] . . . . . . . . . M ( ) , uu − dd ( Q ) Figure 6. Q dependence of the isovector M ( ) ,uu − dd ( Q ) mo-ment. Black data points (tabulated in Table III) are obtainedfrom independent fits to the ω -dependence of F ( ω, Q ) atfixed Q . Curve shows the fit to Equation (43). The phenomenological values of the moments in thelanguage of the parton model, namely v n ( µ ) in Equa-tion (20), are commonly quoted at the scale µ = M ( ) ,uu − dd ( Q ) , say at Q ≳ , then performing a perturbative rescaling down to µ = Q and a further increase in statistics.At this stage we refrain from extending the presentanalysis to the higher moments, however an immediateavenue of study will be to investigate if the observed en-hancement of the M ( ) isovector moment persists withhigher statistics over a larger Q range. Should this bethe case, it will be interesting to compare with the en-hancement observed in the empirical results at small- Q (see e.g. [61]).In addition to studying the Q -dependence of the uu , dd and uu − dd moments of F ( ω, Q ) , the techniquesdescribed in this work can be easily extended to investi-gations of additional quantities such as the higher-twist ud moments (see the final two diagrams in Figure 1) or atest of the Callen-Gross relation at small Q . For a firstattempt see [50]. We emphasize that this work clearlydemonstrates that a study of the Q -dependence of suchobservable on the lattice is now possible. VI. SUMMARY AND CONCLUSIONS
We have presented a derivation of the second-orderFeynman-Hellmann theorem and its relationship to theforward Compton amplitude. In particular, the Comp-ton amplitude can be computed directly on the latticewith a simple extension of the already established latticeimplementations of the Feynman-Hellmann theorem, de-void of operator mixing and related complications of the conventional approach. In order to illustrate the feasi-bility and the versatility of this method in directly prob-ing nucleon structure functions, we have performed high-statistics simulations for several photon momenta, Q , onthe 2 + ×
64 QCDSF/UKQCD lattices at the SU ( ) flavor symmetric point corresponding to a pionmass of ≃
470 MeV. By studying the Compton ampli-tude across a range of kinematics, we have presented non-trivial signals for the first few moments of the nucleonstructure functions. By revealing the Q dependence ofthe low moments, there is a clear opportunity to directlystudy the evolution to the partonic regime — more de-tailed investigations of the power corrections will be pur-sued in future work. Beyond studying the approach tothe partonic regime, the method could also be applied atsmaller Q and probe the dynamics of low-energy Comp-ton scattering processes [62]. While moments of structurefunctions are relatively straightforward, there is also aprospect to invert the Compton amplitude to extract the x -dependence of the structure functions directly — seeRef. [63] for some first attempts. ACKNOWLEDGMENTS
The numerical configuration generation (using theBQCD lattice QCD program [64])) and data analysis (us-ing the Chroma software library [65]) was carried out onthe IBM BlueGene/Q and HP Tesseract using DIRAC2 resources (EPCC, Edinburgh, UK), the IBM Blue-Gene/Q (NIC, J¨ulich, Germany) and the Cray XC40at HLRN (The North-German Supercomputer Alliance),the NCI National Facility in Canberra, Australia (sup-ported by the Australian Commonwealth Government)and Phoenix (University of Adelaide). RH is supportedby STFC through grant ST/P000630/1. HP is supportedby DFG Grant No. PE 2792/2-1. PELR is supported inpart by the STFC under contract ST/G00062X/1. GSis supported by DFG Grant No. SCHI 179/8-1. KUC,RDY and JMZ are supported by the Australian ResearchCouncil grant DP190100297.
Appendix A: Operator product expansion
For completeness, we summarize the relevant notationto complement the OPE discussion in Section II C. The(leading-order) Wilson coefficients are given by: C ( j ) f, n = Q f + O( g ) , j = , ⟨ p, s ∣ [O { µ ...µ n } f − Tr ] ∣ p, s ⟩ = v fn [ p µ . . . p µ n − Tr ] , (A2)in terms of the traceless and symmetric parts of the localquark bilinears: O { µ ...µ n } q = i n − ψ q γ µ ←→ D µ ⋯←→ D µ n ψ q , (A3)1and similarly for the gluons O { µ ...µ n } g = i n − Tr F µ ν ←→ D µ ⋯←→ D µ n F µ n ν , (A4)where ←→ D = (—→ D − ←— D ) . Appendix B: Second order Feynman-HellmannTheorem
For the interested reader, we provide some of the keyintermediate steps to produce the principal derivationpresented in the main text, Equation (23).The 2-point nucleon correlator in an external field,Equation (22), is given by: G ( ) λ ( p ; t ) ≡ ∫ d xe − i p ⋅ x Γ ⟨ Ω λ ∣ χ ( x , t ) χ ( )∣ Ω λ ⟩ , (B1)where Γ the spin-parity projection matrix, with trace im-plied, and ∣ Ω λ ⟩ is the vacuum in the presence of the ex- ternal field. A nucleon interpolating operator is assumedfor χ .The strategy to derive the second-order Feynman-Hellmann relation is to consider the general form of thespectral decomposition of the Euclidean correlator, andmatch the energy shift against the explicit decompo-sition of the correlator in the presence of a weak ex-ternal field. Following the usual procedure of insert-ing a complete set of states in between the operators, ⨋ X d k ( π ) E Xλ ( k ) ∣ X λ ( k )⟩⟨ X λ ( k )∣ , and carrying out themomentum integral, the spectral decomposition of Equa-tion (B1) in the large (Euclidean) time limit, where theground state dominance is realized, is given as, G ( ) λ ( p ; t ) ≃ A λ ( p ) e − E Nλ ( p ) t , (B2)where E N λ ( p ) and A λ ( p ) are the energy of the groundstate nucleon and overlap factor, respectively, in thebackground field. We note that in the presence of thebackground field, the Hamiltonian of the system will mixmomentum states p ± n q —with that p chosen to corre-spond to the lowest kinetic energy, ∣ p ∣ < ∣ p + n q ∣ ∀ n ∈ Z .The second order derivative of Equation (B2) with re-spect to λ , evaluated at λ =
0, is given by: ∂ G ( ) λ ( p ; t ) ∂λ RRRRRRRRRRR λ = = e − E N ( p ) t ⎡⎢⎢⎢⎢⎣ ∂ A λ ( p ) ∂λ − t ( ∂A λ ( p ) ∂λ ∂E N λ ( p ) ∂λ + A ( p ) ∂ E N λ ∂λ ) + t A ( p ) ( ∂E N λ ( p ) ∂λ ) ⎤⎥⎥⎥⎥⎦ . (B3)The derivatives of A λ ( p ) and E N λ ( p ) are understood to be evaluated at λ =
0. The first-order energy shifts vanish, ∂E N / ∂λ =
0, provided we restrict ourselves to the non-Breit-frame kinematics, i.e. ∣ p ∣ ≠ ∣ p ± q ∣ [27, 30]. In this case,the above equation thus reduces to ∂ G ( ) λ ( p ; t ) ∂λ RRRRRRRRRRR λ = = e − E N ( p ) t [ ∂ A λ ( p ) ∂λ − tA ( p ) ∂ E N λ ( p ) ∂λ ] . (B4)where the first term corresponds to the shift in the overlap factor and the second order energy shift is identified inthe t -enhanced or the time-enhanced term. The familiar overlap factor is given by: A ( p ) = ∑ s E N ( p ) Γ ⟨ Ω ∣ χ ( )∣ N ( p , s )⟩⟨ N ( p , s )∣ χ ( )∣ Ω ⟩ . (B5)We now directly evaluate the second-order derivative within the path integral formalism. The 2-point correlationfunction takes the form: λ ⟨ χ ( x , t ) χ ( )⟩ λ = Z( λ ) ∫ D ψ D ψ D U χ ( x , t ) χ ( ) e − S ( λ ) , (B6)where S ( λ ) is the perturbed action given in Equation (22), and Z( λ ) is the corresponding partition function. Pro-jecting the 2-point function to definite momenta and spin gives the standard correlator, G ( ) λ ( p ; t ) = ∫ d x e − i p ⋅ x Γ λ ⟨ χ ( x , t ) χ ( )⟩ λ (B7)To simplify the following expressions, we use the shorthand notation to describe the product of interpolatingoperators, G = ∫ d x e − i p ⋅ x Γ χ ( x , t ) χ ( ) . (B8)2 𝜒𝜒𝜒 𝜒𝜒𝜒 JJ 𝜒 A JJ 𝜒 B JJ 𝜒 C JJ 𝜒 D JJ 𝜒 E JJ 𝜒 F Figure 7. Distinct time orderings of the current insertions, with increasing time assumed from left to right.
The first-order derivative of the correlator is then given by ∂ ⟨G⟩ λ ∂λ = ⟨G⟩ λ ⟨ ∂S ( λ ) ∂λ ⟩ λ − ⟨G ∂S ( λ ) ∂λ ⟩ λ . (B9)The first term corresponds to a vacuum shift and the second term encodes a the three-point correlation functionthat is related to the first-order energy shift. This term has been discussed in detail and applied to the calculationof forward matrix elements [27] and form factors [30]. For the Compton amplitude, the second order derivative isrequired, which is straightforward to evaluate, ∂ ⟨G⟩ λ ∂λ = ⟨G⟩ λ ⟨ ∂ S ( λ ) ∂λ ⟩ λ + ⟨G ∂ S ( λ ) ∂λ ⟩ λ + ⟨G⟩ λ ⟨( ∂S ( λ ) ∂λ ) ⟩ λ + ⟨G⟩ λ ⟨ ∂S ( λ ) ∂λ ⟩ λ ⟨ ∂S ( λ ) ∂λ ⟩ λ − ⟨G ∂S ( λ ) ∂λ ⟩ λ ⟨ ∂S ( λ ) ∂λ ⟩ λ + ⟨G ( ∂S ( λ ) ∂λ ) ⟩ λ . (B10)The first two terms vanish when the external perturbation is purely linear in λ . In the limit λ →
0, vacuum matrixelements of the external fields vanish, ⟨ ∂S ( λ )/ ∂λ ⟩ =
0, assuming the operator doesn’t carry vacuum quantum numbers,such as the electromagnetic current—the scalar current would be an obvious counter example. The term involving ⟨( ∂S ( λ )/ ∂λ ) ⟩ will not in general vanish, however this can only act as a multiplicative factor on the free-field correlatorand hence cannot contribute to the time-enhanced term in Equation (B4). The second-order energy shift can thereforeonly arise from the final term in Equation (B10), ∂ ⟨G⟩ λ ∂λ ∣ λ = = ⟨G ( ∂S ( λ ) ∂λ ) ⟩ + . . . , (B11)where the ellipsis denotes terms that are not time-enhanced. By restoring the explicit form for G , we have ∂ G ( ) λ ( p ; y ) ∂λ RRRRRRRRRRR λ = = ∫ d x e − i p ⋅ x Γ ⟨ χ ( x , t ) χ ( ) ( ∂S ( λ ) ∂λ ) ⟩ , (B12)Using our explicit form for the electromagnetic external field, the corresponding second derivative of the correlatoris given by ∂ G ( ) λ ( p ; t ) ∂λ RRRRRRRRRRR λ = = ∫ d xe − i p ⋅ x Γ ∫ d yd z ( e i q ⋅ y + e − i q ⋅ y )( e i q ⋅ z + e − i q ⋅ z )⟨ χ ( x , t )J µ ( z )J µ ( y ) χ ( )⟩ . (B13)The correlator defined here involves a four-point correlation function with nucleon interpolating operators held atfixed temporal separation t , with the currents inserted across the entire four-volume. Importantly, this expressionis evaluated in the absence of the external field, and hence momentum conservation is exact. It is then possible toperform a spectral decomposition of this correlator in terms of a transfer matrix that is diagonal in the momenta.Given that the Fourier projection of the nucleon sink is at definite momentum p , and ∣ p ∣ < ∣ p ± q ∣ (as discussedabove Equation (B4)), the leading asymptotic behavior of the correlator must have an exponential behavior given by e − E N ( p ) t . By resolving the corresponding t -enhanced coefficient of this exponential, we can identify the second-orderenergy shift, as given in Equation (B4).Assuming that the temporal length is sufficiently large that we can neglect the temporal boundary conditions, thereare six distinct time orderings of where the current insertions can act relative to the nucleon interpolating fields. Theyare shown in Figure 7. Configuration A is the obvious ordering that contains the desired Compton amplitude. This3corresponds to ground-state saturation of the nucleon on either side of the current insertions. The t dependence ofthis particular contribution, including explicit integrals over the current insertion times, will take the form: ∫ t dτ ′ ∫ τ ′ dτ ⟨ χ ( t ) J ( τ ′ ) J ( τ ) χ ( )⟩ ∼ ∫ t dτ ′ ∫ τ ′ dτ e − E N ( p )( t − τ ′ ) e − E X ( p + q )( τ ′ − τ ) e − E N ( p )( τ ) (B14)It is convenient to isolate the current separation time by transforming the coordinates to:∆ = τ ′ − τ (B15)¯ τ = ( τ + τ ′ )/ , (B16)and hence ∫ t dτ ′ ∫ τ ′ dτ e − E N ( p )( t − τ ′ ) e − E X ( p + q )( τ ′ − τ ) e − E N ( p )( τ ) = ∫ t d ∆ ∫ t − ∆ / / d ¯ τ e − E N ( p ) t e −( E X ( p + q )− E N ( p )) ∆ , (B17) = e − E N ( p ) t ∫ t d ∆ e −( E X ( p + q )− E N ( p )) ∆ ( t − ∆ ) . (B18)The term linear in t corresponds to the anticipated time enhancement of Equation (B4)—details of the connectionto the Compton amplitude are given below. Given the condition that E X > E N , the damping ensures that the termproportional to ∆ is independent of t for large times. It is this damping which ensures the current separation remainslocalized in time, and allows the nucleon to saturate to the ground state on either side of the current.Having selected the term of interest, it is necessary to confirm that none of the other possible configurations canscale as te − E N ( p ) t at large times. One potential example would be to consider the nucleon at the source to carrymomentum p + q . This case gives a temporal behavior according to ∫ t dτ ′ ∫ τ ′ dτ ⟨ χ ( t ) J ( τ ′ ) J ( τ ) χ ( )⟩ ∼ ∫ t dτ ′ ∫ τ ′ dτ e − E N ( p )( t − τ ′ ) e − E X ( p + q )( τ ′ − τ ) e − E N ( p + q )( τ ) , (B19) = ∫ t d ∆ ( e − E N ( p + q ) t e −( E X ( p + q )− E N ( p + q )) ∆ − e − E N ( p ) t e −( E X ( p + q )− E N ( p )) ∆ E N ( p + q ) − E N ( p ) ) . (B20)The second term clearly contains a damped exponential and hence the integral over ∆ converges for large t . In thefirst term, the ordering of the levels E X ( p + q ) and E N ( p + q ) will govern which contribution dominates at large t .However in either case, this term is exponentially-suppressed relative to e − E N ( p ) t . This example, that does not exhibitthe desired te − E N ( p ) t behavior, makes it clear that in order to generate the coefficient linear in t , one must have 2intermediate propagators of the lowest energy nucleon, such as in Equation (B14). With three available time windowsand the momentum transfer through the current insertion, it is only possible to achieve this with the lowest-energynucleons separated by an intermediate, energetic state.It is then straightforward to conclude that it is not possible for any of the temporal configurations B to F togenerate a contribution te − E N ( p ) t . To highlight how these other terms contribute, we consider the behavior of theB-type ordering. One of the contributions would take the form ∫ ∞ t dτ ′ ∫ t dτ ⟨ J ( τ ′ ) χ ( t ) J ( τ ) χ ( )⟩ ∼ ∫ ∞ t dτ ′ ∫ t dτ e − E V ( q )( τ ′ − t ) e − E X ( p + q )( t − τ ) e − E N ( p ) τ . (B21)Although a “light” vector meson propagates outside the nucleon interpolators, the τ ′ integral is convergent and noremnant of this mass scale can appear in the t -dependent exponent. And even though the momentum states werechosen to highlight a e − E N ( p ) t contribution, there cannot be a temporal enhancement since the kinematics are chosento ensure E X ( p + q ) > E N ( p ) .Given that the contribution to the second-order energy shift must come from the temporal orientation of type A,we demonstrate how this relates to the Compton amplitude. Explicitly written out, configuration A gives rise to the4-point function: ∂ G ( ) λ ( p ; t ) ∂λ RRRRRRRRRRR Aλ = = ∫ d x e − i p ⋅ x ∫ d yd z ∫ t dτ ′ ∫ τ ′ dτ ( e i q ⋅ y + e − i q ⋅ y )( e i q ⋅ z + e − i q ⋅ z )× Γ ⟨ χ ( x )∣J µ ( z , τ ′ )J µ ( y , τ )∣ χ ( )⟩ . (B22)4We insert complete sets of states next to the nucleon interpolating operators, and translate the operator expressionsaccording to the standard form, χ ( x ) = e − i ˆ P .x χ ( ) e i ˆ P .x and J µ ( z )J µ ( y ) = e − i ˆ P .y J µ ( z − y )J µ ( ) e i ˆ P .y , which leads to ∂ G ( ) λ ( p ; t ) ∂λ RRRRRRRRRRR Aλ = = ∫ d yd z ∫ t dτ ′ ∫ τ ′ dτ ∑ X,Y ∫ d k ( π ) e − E X ( p ) t e −( E Y ( k )− E X ( p )) τ E X ( p ) E Y ( k )× e i ( k − p )⋅ y ( e i q ⋅ y + e − i q ⋅ y )( e i q ⋅ z + e − i q ⋅ z )× Γ ⟨ Ω ∣ χ ( )∣ X ( p )⟩⟨ X ( p )∣J µ ( z − y , τ ′ − τ )J µ ( , )∣ Y ( k )⟩⟨ Y ( k )∣ χ ( )∣ Ω ⟩ . (B23)By adopting the transformation, z ′ = z − y , y ′ = y , the Fourier integral over y ′ can be eliminated, and henceeliminate the k integral: ∂ G ( ) λ ( p ; t ) ∂λ RRRRRRRRRRR Aλ = = ∫ d z ′ ∫ t dτ ′ ∫ τ ′ dτ ∑ X,Y e i q ⋅ z ′ e − E X ( p )( t − τ ) E X ( p )× [ e − E Y ( p − q ) τ E Y ( p − q ) Γ ⟨ Ω ∣ χ ( )∣ X ( p )⟩⟨ X ( p )∣J µ ( z ′ , τ ′ − τ )J µ ( , )∣ Y ( p − q )⟩⟨ Y ( p − q )∣ χ ( )∣ Ω ⟩+ e − E Y ( p ) τ E Y ( p ) Γ ⟨ Ω ∣ χ ( )∣ X ( p )⟩⟨ X ( p )∣J µ ( z ′ , τ ′ − τ )J µ ( , )∣ Y ( p )⟩⟨ Y ( p )∣ χ ( )∣ Ω ⟩] + ( q → − q ) (B24)As described above in Equation (B20), the term involving the momentum transfer between in and out states cannotcontribute to the energy shift, it is only the term involving a p → p matrix element that is of interest. By applying theresult of Equation (B17), and noting that at large t , the correlator must be dominated by the state E X = E Y = E N : ∂ G ( ) λ ( p ; t ) ∂λ RRRRRRRRRRR Aλ = = ∫ d z e i q ⋅ z e − E N ( p ) t ( E N ( p )) ∑ s,s ′ ∫ t d ∆ ( t − ∆ )× Γ ⟨ Ω ∣ χ ( )∣ N ( p , s )⟩⟨ N ( p , s )∣J µ ( z , ∆ )J µ ( , )∣ N ( p , s ′ )⟩⟨ N ( p , s ′ )∣ χ ( )∣ Ω ⟩ + ( q → − q ) + . . . , (B25)where the ellipsis represents terms that are suppressed at large t , relative to te − E N ( p ) t . Here, in identifying the ground-state nucleon, the spin sums implied by the ∑ X,Y have been restored. Because the matrix element ⟨ N ∣ J ( ∆ ) J ( )∣ N ⟩ is exponentially damped at large ∆, the term in the integrand proportional to ∆ also cannot generate a contributionto the second-order energy shift. Hence the only remaining term contributing to the energy shift is: ∂ G ( ) λ ( p ; t ) ∂λ RRRRRRRRRRR Aλ = = ∑ ss ′ A ss ′ ( p ) t e − E N ( p ) t E N ( p ) [∫ t d ∆ ∫ d z e i q ⋅ z ⟨ N ( p , s )∣J µ ( z , ∆ )J µ ( , )∣ N ( p , s ′ )⟩ + ( q → − q )] + . . . , (B26)where a spin-density overlap is used: A ss ′ ( p ) = E N ( p ) Γ ⟨ Ω ∣ χ ( )∣ N ( p , s )⟩⟨ N ( p , s ′ )∣ χ ( )∣ Ω ⟩ = δ ss ′ A ( p ) . (B27)A comparison of the form presented in Equation (B26) with Equation (B4), together with the Compton amplitude inEquations. (15) and (11) with q =
0, yields our result quoted in Equation (23). [1] G. Martinelli and C. T. Sachrajda, Nucl. Phys. B ,660 (1996), arXiv:hep-ph/9605336.[2] M. Beneke, Phys. Rept. , 1 (1999), arXiv:hep-ph/9807443. [3] G. Martinelli and C. T. Sachrajda, Nucl. Phys. B ,355 (1989).[4] M. G¨ockeler, R. Horsley, E.-M. Ilgenfritz, H. Perlt, P. E.Rakow, G. Schierholz, and A. Schiller, Phys. Rev. D , , 107 (2018), arXiv:1711.07916 [hep-ph].[6] X. Ji, Phys. Rev. Lett. , 262002 (2013),arXiv:1305.1539 [hep-ph].[7] H.-W. Lin, J.-W. Chen, S. D. Cohen, and X. Ji, Phys.Rev. D , 054510 (2015), arXiv:1402.1462 [hep-ph].[8] C. Alexandrou, K. Cichy, V. Drach, E. Garcia-Ramos, K. Hadjiyiannakou, K. Jansen, F. Steffens,and C. Wiese, Phys. Rev. D , 014502 (2015),arXiv:1504.07455 [hep-lat].[9] K. Orginos, A. Radyushkin, J. Karpie, andS. Zafeiropoulos, Phys. Rev. D , 094503 (2017),arXiv:1706.05373 [hep-ph].[10] K. Cichy and M. Constantinou, Adv. High Energy Phys. , 3036904 (2019), arXiv:1811.07248 [hep-lat].[11] K.-F. Liu and S.-J. Dong, Phys. Rev. Lett. , 1790(1994), arXiv:hep-ph/9306299.[12] M. G¨ockeler, R. Horsley, E.-M. Ilgenfritz, H. Perlt, P. E.Rakow, G. Schierholz, and A. Schiller, Phys. Rev. D ,5705 (1996), arXiv:hep-lat/9602029.[13] S. Capitani, M. G¨ockeler, R. Horsley, H. Oelrich, D. Pet-ters, P. E. Rakow, and G. Schierholz, Nucl. Phys. BProc. Suppl. , 288 (1999), arXiv:hep-lat/9809171.[14] W. Detmold and C. J. D. Lin, Phys. Rev. D73 , 014501(2006), arXiv:hep-lat/0507007 [hep-lat].[15] Z. Davoudi and M. J. Savage, Phys. Rev. D , 054505(2012), arXiv:1204.4146 [hep-lat].[16] Y.-Q. Ma and J.-W. Qiu, Phys. Rev. Lett. , 022003(2018), arXiv:1709.03018 [hep-ph].[17] X. Ji and C. Jung, Phys. Rev. Lett. , 208 (2001),arXiv:hep-lat/0101014.[18] S. Hashimoto, PTEP , 053B03 (2017),arXiv:1703.01881 [hep-lat].[19] M. T. Hansen, H. B. Meyer, and D. Robaina, Phys. Rev.D , 094513 (2017), arXiv:1704.08993 [hep-lat].[20] X. Feng, M. Gorchtein, L.-C. Jin, P.-X. Ma, andC.-Y. Seng, Phys. Rev. Lett. , 192002 (2020),arXiv:2003.09798 [hep-lat].[21] R. A. Brice˜no, Z. Davoudi, M. T. Hansen, M. R.Schindler, and A. Baroni, Phys. Rev. D , 014509(2020), arXiv:1911.04036 [hep-lat].[22] A. Chambers, R. Horsley, Y. Nakamura, H. Perlt,P. Rakow, G. Schierholz, A. Schiller, K. Somfleth,R. Young, and J. Zanotti, Phys. Rev. Lett. , 242001(2017), arXiv:1703.01153 [hep-lat].[23] S. J. Brodsky, Acta Phys. Polon. B , 635 (2005),arXiv:hep-ph/0411028.[24] F. Hautmann and D. E. Soper, Phys. Rev. D , 074020(2007), arXiv:hep-ph/0702077.[25] W. Bietenholz, V. Bornyakov, M. G¨ockeler, R. Hors-ley, W. G. Lockhart, Y. Nakamura, H. Perlt, D. Pleiter,P. E. L. Rakow, G. Schierholz, A. Schiller, T. Streuer,H. St¨uben, F. Winter, and J. M. Zanotti, Phys. Rev. D , 054509 (2011), arXiv:1102.5300 [hep-lat].[26] R. Horsley, R. Millo, Y. Nakamura, H. Perlt, D. Pleiter,P. Rakow, G. Schierholz, A. Schiller, F. Winter, and J. Zanotti (QCDSF, UKQCD), Phys. Lett. B , 312(2012), arXiv:1205.6410 [hep-lat].[27] A. J. Chambers, R. Horsley, Y. Nakamura, H. Perlt,D. Pleiter, P. E. L. Rakow, G. Schierholz, A. Schiller,H. St¨uben, R. D. Young, and J. M. Zanotti (CSSM,QCDSF/UKQCD), Phys. Rev. D , 014510 (2014),arXiv:1405.3019 [hep-lat].[28] A. Chambers, R. Horsley, Y. Nakamura, H. Perlt,P. Rakow, G. Schierholz, A. Schiller, and J. Zanotti(QCDSF), Phys. Lett. B , 30 (2015), arXiv:1410.3078[hep-lat].[29] A. Chambers, R. Horsley, Y. Nakamura, H. Perlt,D. Pleiter, P. Rakow, G. Schierholz, A. Schiller,H. St¨uben, R. Young, , and J. Zanotti, Phys. Rev. D , 114517 (2015), arXiv:1508.06856 [hep-lat].[30] A. J. Chambers, J. Dragos, R. Horsley, Y. Nakamura,H. Perlt, D. Pleiter, P. E. L. Rakow, G. Schierholz,A. Schiller, K. Somfleth, H. St¨uben, R. D. Young, andJ. M. Zanotti (QCDSF, UKQCD, CSSM), Phys. Rev. D , 114509 (2017), arXiv:1702.01513 [hep-lat].[31] W. Detmold, Phys. Rev. D , 054506 (2005), arXiv:hep-lat/0410011.[32] T. Primer, W. Kamleh, D. Leinweber, and M. Burkardt,Phys. Rev. D , 034508 (2014), arXiv:1307.1509 [hep-lat].[33] Z. Davoudi and W. Detmold, Phys. Rev. D , 074506(2015), arXiv:1507.01908 [hep-lat].[34] M. J. Savage, P. E. Shanahan, B. C. Tiburzi, M. L. Wag-man, F. Winter, S. R. Beane, E. Chang, Z. Davoudi,W. Detmold, and K. Orginos, Phys. Rev. Lett. ,062002 (2017), arXiv:1610.04545 [hep-lat].[35] C. Bouchard, C. C. Chang, T. Kurth, K. Orginos,and A. Walker-Loud, Phys. Rev. D , 014504 (2017),arXiv:1612.06963 [hep-lat].[36] A. Agadjanov, U.-G. Meißner, and A. Rusetsky, Phys.Rev. D , 031502 (2017), arXiv:1610.05545 [hep-lat].[37] R. Devenish, Deep inelastic scattering (Oxford UniversityPress, Oxford, 2004).[38] A. V. Manohar, in
Lake Louise Winter Institute: Sym-metry and Spin in the Standard Model Lake Louise, Al-berta, Canada, February 23-29, 1992 (1992) pp. 1–46,arXiv:hep-ph/9204208 [hep-ph].[39] K.-F. Liu, Phys. Rev.
D62 , 074501 (2000), arXiv:hep-ph/9910306 [hep-ph].[40] D. Drechsel, B. Pasquini, and M. Vanderhaeghen, Phys.Rept. , 99 (2003), arXiv:hep-ph/0212124.[41] R. A. Brice˜no, M. T. Hansen, and C. J. Monahan, Phys.Rev. D , 014502 (2017), arXiv:1703.06072 [hep-lat].[42] J. Liang, T. Draper, K.-F. Liu, A. Rothkopf, and Y.-B. Yang (XQCD), Phys. Rev. D , 114503 (2020),arXiv:1906.05312 [hep-ph].[43] A. H. Gunn et al. , in preperation (2020).[44] N. Cundy, M. G¨ockeler, R. Horsley, T. Kaltenbrunner,A. D. Kennedy, Y. Nakamura, H. Perlt, D. Pleiter,P. E. L. Rakow, A. Sch¨afer, G. Schierholz, A. Schiller,H. St¨uben, and J. M. Zanotti, Phys. Rev. D , 094507(2009), arXiv:0901.3302 [hep-lat].[45] W. Bietenholz, V. Bornyakov, M. G¨ockeler, R. Hors-ley, W. G. Lockhart, Y. Nakamura, H. Perlt, D. Pleiter,P. E. L. Rakow, G. Schierholz, A. Schiller, T. Streuer,H. St¨uben, F. Winter, and J. M. Zanotti (QCDSF-UKQCD Collaboration), Phys. Rev. D , 054509 (2011).[46] W. Bietenholz, V. Bornyakov, N. Cundy, M. G¨ockeler,R. Horsley, A. D. Kennedy, W. G. Lockhart, Y. Naka- mura, H. Perlt, D. Pleiter, P. E. L. Rakow, A. Sch¨afer,G. Schierholz, A. Schiller, H. St¨uben, and J. M. Zanotti,Phys. Lett. B690 , 436 (2010), arXiv:1003.1114 [hep-lat].[47] R. Horsley, J. Najjar, Y. Nakamura, H. Perlt, D. Pleiter,P. E. L. Rakow, G. Schierholz, A. Schiller, H. St¨uben,and J. M. Zanotti (QCDSF-UKQCD),
Proceedings, 31stInternational Symposium on Lattice Field Theory (Lat-tice 2013): Mainz, Germany, July 29-August 3, 2013 ,PoS
LATTICE2013 , 249 (2014), arXiv:1311.5010 [hep-lat].[48] V. Bornyakov, R. Horsley, R. Hudspith, Y. Nakamura,H. Perlt, D. Pleiter, P. Rakow, G. Schierholz, A. Schiller,H. St¨uben, and J. Zanotti, (2015), arXiv:1508.05916[hep-lat].[49] M. Constantinou, R. Horsley, H. Panagopoulos, H. Perlt,P. Rakow, G. Schierholz, A. Schiller, and J. Zanotti,Phys. Rev. D , 014502 (2015), arXiv:1408.6047 [hep-lat].[50] A. Hannaford-Gunn, R. Horsley, Y. Nakamura, H. Perlt,P. E. L. Rakow, G. Schierholz, K. Somfleth, H. St¨uben,R. D. Young, and J. M. Zanotti, in (2020)arXiv:2001.05090 [hep-lat].[51] C. R. Allton, C. T. Sachrajda, R. M. Baxter, S. P. Booth,K. C. Bowler, S. Collins, D. S. Henty, R. D. Kenway, B. J.Pendleton, D. G. Richards, J. N. Simone, A. D. Simpson,B. E. Wilkes, and C. Michael (UKQCD), Phys. Rev. D , 5128 (1993), arXiv:hep-lat/9303009.[52] F. Hausdorff, Mathematische Zeitschrift , 74 (1921).[53] J. Salvatier, T. V. Wiecki, and C. Fonnesbeck, PeerJComputer Science (2016).[54] M. D. Hoffman and A. Gelman, Journal of Machine Learning Research , 1593 (2014).[55] E. D. Bloom and F. J. Gilman, Phys. Rev. Lett. , 1140(1970).[56] O. Nachtmann, Nucl. Phys. B , 237 (1973).[57] M. Kuroda and G. Schierholz, Nucl. Phys. B , 148(1979).[58] J. M. Cornwall and R. E. Norton, Phys. Rev. , 2584(1969).[59] S. Capitani, M. G¨ockeler, R. Horsley, D. Petters,D. Pleiter, P. Rakow, and G. Schierholz, Nuclear PhysicsB - Proceedings Supplements , 173 (1999), proceed-ings of the 7th International Workshop on Deep InelasticScattering and QCD.[60] W. Bietenholz, N. Cundy, M. Goeckeler, R. Horsley,H. Perlt, D. Pleiter, P. Rakow, G. Schierholz, A. Schiller,T. Streuer, and J. Zanotti, PoS LAT2009 , 138 (2009),arXiv:0910.2437 [hep-lat].[61] W. Melnitchouk, R. Ent, and C. Keppel, Phys. Rept. , 127 (2005), arXiv:hep-ph/0501217.[62] V. Lensky, F. Hagelstein, V. Pascalutsa, andM. Vanderhaeghen, Phys. Rev. D , 074012 (2018),arXiv:1712.03886 [hep-ph].[63] R. Horsley, Y. Nakamura, H. Perlt, P. E. L. Rakow,G. Schierholz, K. Somfleth, R. D. Young, andJ. M. Zanotti (QCDSF-UKQCD-CSSM), in (2020)arXiv:2001.05366 [hep-lat].[64] T. R. Haar, Y. Nakamura, and H. Stuben, EPJ WebConf. , 14011 (2018), arXiv:1711.03836 [hep-lat].[65] R. G. Edwards and B. Joo (SciDAC Collabora-tion, LHPC Collaboration, UKQCD Collaboration),Nucl.Phys.Proc.Suppl.140