# A Preliminary Determination of the Second Mellin Moment of the Pion's Distribution Amplitude Using the Heavy Quark Operator Product Expansion

William Detmold, Anthony V. Grebe, Issaku Kanamori, C.-J. David Lin, Santanu Mondal, Robert J. Perry, Yong Zhao

AA Preliminary Determination of the Second Mellin Moment of the Pion’s DistributionAmplitude Using the Heavy Quark Operator Product Expansion

William Detmold and Anthony V. Grebe ∗ Center for Theoretical Physics, Massachusetts Institute of Technology,77 Massachusetts Avenue, Cambridge, MA 02139, USA

Issaku Kanamori

RIKEN Center for Computational Science, Kobe 650-0047, Japan

C.-J. David Lin

Institute of Physics, National Chiao-Tung University, 1001 Ta-Hsueh Road, Hsinchu 30010, Taiwan andCentre for High Energy Physics, Chung-Yuan Christian University, Chung-Li, 32032, Taiwan

Santanu Mondal

Los Alamos National Laboratory, Theoretical Division T-2, Los Alamos, NM 87545, USA

Robert J. Perry † Institute of Physics, National Chiao-Tung University, 1001 Ta-Hsueh Road, Hsinchu 30010, Taiwan

Yong Zhao

Physics Department, Brookhaven National Laboratory, Bldg. 510A, Upton, NY 11973, USA (HOPE Collaboration)

We explore the feasibility of determining Mellin moments of the pion’s light cone distributionamplitude using the heavy quark operator product expansion (HOPE) method. As the ﬁrst step ofa proof of principle study we pursue a determination of the second Mellin moment. We discuss ourchoice of kinematics which allows us to successfully extract the moment at low pion momentum. Wedescribe the numerical simulation, and describe the data analysis, which leads us to a preliminarydetermination of the second Mellin moment in the continuum limit in the quenched approximationas (cid:10) ξ (cid:11) = 0 . I. INTRODUCTION

At high energies, exclusive processes in quantum chro-modynamics (QCD) may be described with the aid ofthe so-called light cone distribution amplitudes (LCDAs),convolved with a short distance perturbative kernel [1,2]. These distribution amplitudes contain the non-perturbative information about the process, and are theresult of a Fock space truncation where only the lowest,valence Fock state is retained. In this sense, one may con-sider the pion’s light cone distribution amplitude φ ( x, µ )as giving the probability amplitude for converting a pioninto a pair of collinear quark and antiquark with longi-tudinal momentum fraction x and 1 − x , respectively. Itis deﬁned via the matrix element (cid:104) | ψ ( z n ) /nγ W [ z n, z n ] ψ ( z n ) (cid:12)(cid:12) π + ( p ) (cid:11) = if π ( p · n ) (cid:90) dx e − i ( z x + z (1 − x )) p · n φ π ( x, µ ) , (1)where p µ is the momentum of the pion, n is a light-like vector ( n = 0), z and z are real numbers, f π = ∗ [email protected], Speaker at APLAT 2020 † [email protected], Speaker at APLAT 2020 .

132 GeV is the pion decay constant, µ is the renor-malization scale and W [ a, b ] is a Wilson line required toensure gauge invariance of the matrix element. We take x to be the longitudinal momentum fraction of the u -quarkin the Fock state (cid:12)(cid:12) ud (cid:11) . Momentum conservation then im-plies that the d quark has momentum fraction 1 − x . Inthe isospin limit, where the masses of the up and downquarks are degenerate, the light cone distribution ampli-tude is symmetric under the interchange x → − x , thatis φ π ( x, µ ) = φ π (1 − x, µ ) . (2)We shall assume isospin symmetry in this work.The value of precise determinations of the LCDA liesin the object’s process independence . This allows one todescribe many exclusive processes in QCD with the samedistribution amplitude convolved with a process depen-dent perturbative kernel.Currently, the only ab-initio approach to the determi-nation of this object is a numerical computation usinglattice QCD. Unfortunately, direct calculation of suchlight-cone objects is impossible on a Euclidean latticesince the light cone, deﬁned by z = 0, is a reduced toa single point ( z µ E = 0). Despite this diﬃculty informa-tion about the LCDA may be determined from the latticein a number of indirect ways. The traditional approach a r X i v : . [ h e p - l a t ] O c t involves a determination of the Mellin moments of theLCDA [3, 4]. These are deﬁned by (cid:104) ξ n (cid:105) µ = (cid:90) − dξ ξ n φ ( ξ, µ ) , (3)where ξ = 2 x − x is the momentum fraction of thecollinear quark anti-quark pair. Noting again the isospinsymmetry x ↔ (1 − x ) we see that only the even momentsmay be non-zero for the pion. These moments may berelated to local matrix elements which are immediatelyamenable to a lattice calculation. It is possible to writethe full distribution amplitude with the knowledge of theMellin moments alone: φ ( ξ, µ ) = 12 π (cid:90) ∞−∞ ds (cid:32) ∞ (cid:88) n =0 ( is ) n n ! (cid:104) ξ n (cid:105) µ (cid:33) e − isξ . (4)Unfortunately, the breaking of the full rotation groupon the lattice leads to operator mixing and thus powerdivergences appear in twist-2 operators with spin higherthan four [5]. These power divergences make the determi-nation of the higher moments more diﬃcult. Neverthe-less, this approach has been well studied and has yieldedresults for the ﬁrst non-trivial moment of the pion andkaon [6, 7]. A number of other proposals in the literatureseek to overcome this diﬃculty [8–28].While much good work has been done in the extrac-tion of the pion LCDA, it is clear that more must stillbe done to acquire precise predictions of this object.With this view, it is clearly of interest to explore otherproposals for the calculation of the distribution ampli-tude. One such approach, which we pursue in this workis the so-called heavy quark operator product expansion(HOPE) [29, 30]. The HOPE method builds on the con-ventional operator product expansion (OPE) approach,by performing the numerical simulation with a ﬁctitiousheavy quark species, which leads to a number of advan-tages over the standard treatment [29]. This method al-lows the extraction of the Mellin moments of the LCDA,and thus in principle allows the reconstruction of the am-plitude within a wide range of x . In this paper, we dis-cuss the application of the HOPE method to the pion’sLCDA. In particular, we discuss kinematic choices whichlead to an eﬃcient extraction of the second Mellin mo-ment, and discuss the resulting preliminary extraction ofthe second Mellin moment. II. SUMMARY OF THE CALCULATION

The HOPE method is a multi-step procedure. Thus,before beginning our discussion of the kinematics usedand the numerical study we performed, we provide anoverview of the calculation. The starting point of thiswork is a study of the anti-symmetric version of the ma- trix element in Minkowski space, T µν ( p, q ) = (cid:90) d z e iq · z (cid:104) | T [ J µA ( z/ J νA ( − z/ | π ( p ) (cid:105) , (5)given by U µν ( p, q ) = 12 (cid:18) T µν ( p, q ) − T νµ ( p, q ) (cid:19) = (cid:90) d z e iq · z (cid:104) | T [ J [ µA ( z/ J ν ] A ( − z/ | π ( p ) (cid:105) , (6)where the axial-vector current is replaced by the heavy-light ﬂavour changing current: J µA = ¯Ψ γ µ γ ψ + ¯ ψγ µ γ Ψ , (7)with ψ being the light quark ﬁeld, and Ψ being the heavyquark ﬁeld. We note that it is also possible to study theLCDA Mellin moments using the corresponding heavy-light vector current. By applying the OPE to the abovematrix element, we can show [30] that to leading twist,the antisymmetric tensor U µν ( p, q ) may be written in theisospin limit as U µν ( p, q ) = 2 if π (cid:15) µναβ q α p β ˜ Q ∞ (cid:88) n even C n ( η )2 n ( n + 1) C ( n ) W ( ˜ Q ) × (cid:104) ξ n (cid:105) ζ n + O (1 / ˜ Q ) , (8)where f π ≈ .

132 GeV is the pion decay constant, ˜ Q , η and ζ are kinematic variables given by˜ Q = − q − m ψ , (9) η = p · q (cid:112) p q , (10) ζ = (cid:112) p q ˜ Q . (11) C ( n ) W are the Wilson coeﬃcients, and C n ( η ) are the Gegen-bauer polynomials, which arise as a result of resummingtarget mass eﬀects [31, 32].In order to accurately extract the Mellin moments, oneneeds to determine the Wilson coeﬃcients beyond zerothorder. Since these Wilson coeﬃcients only account forthe ultraviolet eﬀects of QCD, they may be calculatedusing perturbation theory. The Wilson coeﬃcients maybe written C ( n ) W ( ˜ Q ) = 1 + α s c (1) n + . . . . (12) Note that Ref. [30], uses a normalization for the Mellin mo-ments which diﬀers by a factor of 2 n from our convention. Oursagrees with the ‘standard’ normalization which allows us to di-rectly compare our result with other determinations of the secondMellin moment. By calculating the matrix element T µν ( p, q ) on the lat-tice, one may then perform a ﬁt to the form of the heavyquark OPE, and thus obtain the Mellin moments of thedistribution amplitude. We note that when computingthe hadronic tensor, kinematics should be chosen suchthat we remain in the unphysical region. This requireschoosing ( p + q ) < m hl ≈ ( m Ψ + Λ QCD ) , (13)where m hl is the mass of the lightest heavy-light me-son. This ensures that the analytic continuuation toMinkowski space is may be straightforwardly obtainedby the replacement q → iq . Thus this method is con-strained to work in the windowΛ QCD (cid:28) (cid:112) q < m Ψ (cid:28) a (14)By performing the calculation at a number of lattice spac-ings, one may then extrapolate to the continuum in theusual way [33]. A. Eﬃcient Kinematics for an Extraction of theSecond Mellin Moment

In this work, we are primarily interested in an extrac-tion of the second Mellin moment of the pion’s LCDA.From Eqn. 8, it is possible to see that the n th moment isweighted by the kinematical factor C n ( η )2 n ( n + 1) C ( n ) W ( ˜ Q ) ζ n . (15)For this section, we shall assume that the Wilson coeﬃ-cients are unity. This results in O ( α S ) errors, but will noteﬀect the features discussed here. We deﬁne the weightfunction W ( n ) = C n ( η )2 n ( n + 1) ζ n . (16)This weighting factor is the origin of the diﬃculty in ex-tracting the higher moments in OPE approaches. Forexample in our numerical work we ﬁx the physical size ofthe system to be L × a = 1 .

92 fm for all choices of thelattice spacing, a . Thus the smallest unit of momentumis ∆ p = 2 πLa = 0 .

64 GeV , (17)with the pion at rest p = (0 , , × .

64 GeV and thecurrent insertion momentum q = (0 , , × .

64 GeVwith m π = 0 .

56 GeV and m Ψ = 2 . q max[ W (0) , q ] = 1 (18)max[ W (2) , q ] = 0 . , (19) p =( ) , q =( ) p =( ) , q =( ,0,1 ) p =( ) , q =( ) - - R e [ W ( n ) ] p =( ) , q =( ,0,1 ) p =( ) , q =( ) I m [ W ( n ) ] FIG. 1. We demonstrate the suppression factor obtained fromthe kinematic weight function W ( n ) for several diﬀerent kine-matic choices. We examine the combinations p = (0 , , q = (0 , ,

1) (blue circles), p = (1 , , q = (1 / , ,

1) (earthsquares) and p = (4 , , q = (2 , ,

1) (garnet diamonds),where all momenta should be multiplied by 0.64 GeV to deter-mine their physical values. As we explain, in the case where p · q (cid:54) = 0, the weighting function will be complex. Since p · q = 0 for the kinematic choice described by the blue cir-cles, there is no imaginary part, and we thus exclude thosepoints from the lower plot for clarity. with higher moments further suppressed. Under nor-mal circumstances, an extraction of even the ﬁrst non-trivial (ie, the second) moment with this particular choiceof kinematics would be a challenging task. Note how-ever that by changing the kinematics, one may reducethe kinematic suppression. This fact is demonstrated inFig. 1, where a number of diﬀerent choices of kinemat-ics are shown. We note that in general the extraction ofhigher Mellin moments requires higher pion momentum,which poses a challenge for numerical determinations.Since we wish to numerically simulate the Comptontensor so that we may determine the second Mellin mo-ment, it is advantageous to explore our kinematic op-tions to best optimize the desired signal. In particular,by studying the properties of this weight function, wecan determine kinematics which allow us direct accessto the second Mellin moment, somewhat bypassing thekinematical suppression.To begin, we note that while this function is real inMinkowski space, it is in general complex in Euclideanspace. This is because p = iE π ( p ), and we take q real.Noting again the deﬁnitions of the kinematical variables η = p · q (cid:112) p q , (20) ζ = (cid:112) p q ˜ Q . (21)We see that while ζ n is always real, under certain kine-matical choices, η is complex: η = p · q (cid:112) p q + i E π ( p ) q (cid:112) p q . (22)Note that only the even moments are non-zero due toour assumption of isospin symmetry. The correspondingGegenbauer polynomials are also even. Thus we onlyhave even factors of η , and so we see that if we have thespatial inner product p · q (cid:54) = 0, the coeﬃcient of (cid:10) ξ (cid:11) is complex. Note that since the kinematic factors areabsent from the zeroth moment, this allows one separatethe contribution from the lowest moment, and thus gaindirect access to the second moment.There are several caveats to this. Firstly, the overallnormalization of the HOPE can spoil this result. In par-ticular consider the term (cid:15) µναβ q α p β . In this work, westudy the combination µ = 1, ν = 2. We thus have (cid:15) αβ q α p β = q iE π ( p ) − q p . (23)Since this is an overall multiplicative factor, it will in gen-eral imbue all the moments (including the zeroth), witha complex kinematical factor. We can ensure this doesnot occur by taking kinematics where either p = 0 or q = 0. In either case, the kinematic factor will again beeither purely real or purely imaginary, and thus the spe-cial kinematics may be used to directly access the secondmoment. Secondly, in this discussion we have neglectedthe role of the Wilson coeﬃcients, however, we note thatthese can only give corrections which are O ( α S ), andwill be numerically small. Thus as we shall see the ‘spe-cial kinematics’ are still eﬀective in isolating the secondMellin moment. A demonstration of the special kinemat-ics is shown in Fig. 2. To summarize, in this work, weuse the conditions p · q (cid:54) = 0 , (24) p = 0 . (25)In particular, we choose to perform the simulations withthe momentum p = (1 , , × .

64 GeV , (26) q = (1 / , , × .

64 GeV , (27)The reason for the apparent fractional lattice momentumis that as we shall see the ‘physical’ momenta are linear ξ ξ ξ + ξ - - - - - q ( GeV ) R e [ U ( p , q ) ] ( M e V ) ξ ξ ξ + ξ - - - q ( GeV ) I m [ U ( p , q ) ] ( M e V ) FIG. 2. Examining the special kinematics. By choosing thekinematics p = (1 , , × .

64 GeV, q = (1 / , , × .

64 GeVand considering the real and imaginary parts independently,it is possible to see that while the imaginary part is satu-rated with the contribution from the lowest moment, the realpart allows one to directly access the second Mellin momentdirectly. combinations of the inserted momenta, and in particu-lar, we will see that we must include a factor of half inthe deﬁnition of q . This kinematic choice leads to lesskinematical suppression;max[ W (0) , q ] = 1 (28)max[ W (2) , q ] = 0 . , (29)but most importantly allows one direct access to the sec-ond Mellin moment. Having now optimized our kinemat-ical choice, we proceed to discuss the numerical simula-tion, and resulting extraction of the Mellin moment. III. LATTICE COMPUTATION

The hadronic tensor is the Fourier transform of acurrent-current correlator, so it can be written in termsof 2- and 3-point functions of the form C ( τ π , p ) = (cid:90) d x e i p · x (cid:104) |O π ( x , τ π ) O † π ( , | (cid:105) (30) C µν ( τ e , τ m ; p e , p m ) = (cid:90) d x e d x m e i p e · x e e i p m · x m (cid:104) |T (cid:2) J µA ( x e , τ e ) J νA ( x m , τ m ) O † π ( , (cid:3) | (cid:105) (31)where T is the time-ordering operator, O π is the pioninterpolating operator, and J µA ≡ ¯Ψ γ µ γ ψ + ¯ ψγ µ γ Ψ (32)is the ﬂavor-changing axial current insertion operatorthat converts the pion’s light quarks ψ into valence heavyquarks Ψ and vice versa. In the large-time limit, the two-and three-point functions asymptote to C ( τ π , p ) ∼ (cid:12)(cid:12) (cid:104) π ( p ) |O † π ( , | (cid:105) (cid:12)(cid:12) E π e − E π τ π (33) C ( τ e , τ m ; p e , p m ) ∼ (cid:104) π ( p ) |O † π ( , | (cid:105) E π e − E π ( τ e + τ m ) / (cid:90) d z e i q · z (cid:104) |T (cid:104) J µA (cid:16) z (cid:17) J νA (cid:16) − z (cid:17)(cid:105) | π ( p ) (cid:105) (34)with z = x e − x m , (35) p = p e + p m , (36) q = 12 ( p m − p e ) . (37)The 3-point correlator is shown graphically in Fig. 3 andcan be computed via a sequential propagator through theoperator. The source and sink of the 2-point functionand the source of the 3-point function are constructedusing both Gaussian and link smearing to suppress ex-cited state contamination. Fitting C , C µν at large τ π , τ e , and τ m lets us extract R µν ( τ ; p , q ) = (cid:90) d z e i q · z (cid:104) |T (cid:104) J µ (cid:16) z (cid:17) J ν (cid:16) − z (cid:17)(cid:105) | π ( p ) (cid:105) (38)The hadronic tensor is then deﬁned as the Fourier trans-form of R µν in the τ = z direction: U µν ( p, q ) ≡ (cid:90) dτ e iq τ R [ µν ] ( τ ; p , q ) (39) A. Lattice Parameters

In this study, we used Chroma [34] to measure cor-relators at two heavy quark masses, with bare quarkmasses of about 1.6 GeV and 2.5 GeV. (The renor-malized quark masses found by the ﬁts are substantiallyheavier than the bare masses.) To accomodate such largequark masses, we use ﬁne lattice spacings, ranging from0.041 fm to 0.060 fm (and with future plans to includean additional ensemble with a = 0 .

030 fm). The physicalvolumes are tuned to about 1.9 fm on all ensembles.

FIG. 3. The three-point correlation function used to com-pute the hadronic tensor of the pion. The pion is created atthe origin and ﬂavor-changing axial currents are inserted attimes τ e , τ m . The quark propagating between the currentsis artiﬁcally heavy due to the ﬂavor-changing nature of thecurrents. Due to critical slowing down, using dynamical fermionswould be prohibitively expensive, especially for a pre-liminary simulation. Therefore, this calculation is per-formed in the quenched approximation using latticesfrom Ref. [35]. The ensembles used here and the mea-surements performed on them are summarized in Table I.We use Wilson clover fermions with the clover coeﬃ-cient set non-perturbatively to the value in Ref. [36]. Wetuned the pion mass to about 560 MeV across the en-sembles. In addition to reducing the computational cost,this unphysically heavy pion mass ensures that m π L > ×

10 momentum insertions × B. Reducing Noise

At the kinematics used, to O ( α s ), the second moment (cid:104) ξ (cid:105) is proportional to the real part of the hadronic tensor(and the imaginary part of the hadronic tensor is mostlyindependent of (cid:104) ξ (cid:105) ). Thus, measuring Re[ U µν ] gives aclean probe of (cid:104) ξ (cid:105) without much contamination fromhigher-twist eﬀects. However, while this is a clean signal,it is also a small one: At the kinematics used, the realpart of U µν is 2–3 orders of magnitude smaller than theimaginary part.The 3-point correlator (and therefore the ratio R µν ) ispure imaginary , correspond to the antisymmetric and For p, q in Minkowski space, the hadronic tensor is purely imag-inary, as can be seen from the operator product expansion. The

TABLE I. The ensembles used in this study and the number of measurements performed on each. β a (fm) L × T κ light κ heavy, 1 κ heavy, 2 N cfg N src Light Props Heavy Props6.30168 0.060 32 ×

64 0.135146 0.119867 0.112779 450 7 3150 126,0006.43306 0.048 40 ×

80 0.135170 0.122604 0.116599 250 2 500 20,0006.59773 0.041 48 ×

96 0.135028 0.124420 0.119228 341 3 1023 40,920 symmetric parts of R µν . Speciﬁcally,Re[ U µν ( p , q )] = Re (cid:20)(cid:90) ∞−∞ dτ R µν ( τ ; p , q ) e − iq τ (cid:21) ∝ (cid:90) ∞ dτ [ R µν ( τ ; p , q ) − R µν ( − τ ; p , q )] sin( q τ )(40)Thus, we need to measure a diﬀerence R µν ( τ ; p , q ) − R µν ( − τ ; p , q ) that is two orders of magnitude smallerthan each of the terms constituting the diﬀerence. Theprecision to which we can measure this diﬀerence dependson how well the two terms are correlated (which wouldcause correlated errors to cancel). However, for moder-ately large τ , the correlators C µν ( τ e , τ e ± τ ; p e , p m ) usedto compute R µν ( ± τ ; p , q ) have poorly correlated uncer-tainties since the sinks are temporally separated on thelattice.We could obtain better correlations – and there-fore better error cancellation – if we could com-pute R µν ( − τ ; p , q ) using C µν ( τ m , τ e ; p e , p m ), sincethen the correlators used to compute R µν ( τ ; p , q ) and R µν ( − τ ; p , q ) would be at the same timeslices (up to in-terchange of the two current insertions). However, withthe current setup where the sequential propagator passesthrough the ﬁrst current inserted, this would require cur-rent insertions at all desired τ m , which would be pro-hibitively expensive. Instead, we use γ -hermiticity towrite C µν ( τ e , τ m ; p e , p m ) ∗ = C νµ ( τ m , τ e ; − p m , − p e ) (41)where p e and p m are related to p and q via (36), (37).This lets us compute both terms in the right-hand sideof (40) in terms of correlators with τ m ≥ τ e , since R µν ( τ ; p , q ) − R µν ( − τ ; p , q ) = R µν ( τ ; p , q )+ R µν ( τ ; − p , q )(42)Now, the terms in the right-hand side of (42) are morehighly correlated, so we would expect larger cancellationof correlated errors. This eﬀect is shown in Fig. 4, whereuncertainties are reduced by a factor of about 10 by usingthe right-hand side of (42) rather than the left-hand side. Euclidean-space data R µν are related to the Minkowski-spacehadronic tensor via Laplace transform U µν ( q, p ) = (cid:90) ∞−∞ dτ e − q τ R µν ( τ ; q , p )whose kernel is purely real for real q . Thus, if U µν with q ∈ R is imaginary, R µν ( τ ; p , q ) must be too. R μν ( τ ; p , q )+ R μν ( τ ; - p , q ) R μν ( τ ; p , q )- R μν (- τ ; p , q ) τ ( fm ) ⅈ ℱ - [ R e ( U μ ν ) ] ( M e V ) FIG. 4. Comparing both sides of the equality in (42) R µν ( τ ; p , q ) + R µν ( τ ; − p , q ) (blue) and R µν ( τ ; p , q ) − R µν ( − τ ; p , q ) (earth), both measured with two sources on 450conﬁgurations. These quantities agree in expectation, but theformer has uncertainties an order of magnitude smaller thanthe latter. IV. FITTING TO THE HOPE

At the kinematics used here, the n th moment picks upa factor of (cid:16) p · q ˜ Q (cid:17) n (cid:46) . n , so the contribution of fourthmoment is suppressed by a factor of about 50 relative tothat of the second moment. As a result, in this work,we will neglect higher-moment contributions, so we canwrite the operator product expansion as U µν = 2 if π ε µνρλ q ρ p λ ˜ Q (cid:20) C (0) W + (cid:104) ξ (cid:105) p · q ) − p q

6( ˜ Q ) C (2) W + · · · + O (cid:18) Λ QCD ˜ Q (cid:19)(cid:21) (43)where ˜ Q = − m − q , m Ψ is the renormalized heavyquark mass, and C ( n ) W are perturbatively calculable Wil-son coeﬃcients. For this analysis, we have calculated theWilson coeﬃcients to 1-loop order, and we will publishthe results in forthcoming work [33]. The remaining pa-rameters ( f π , m Ψ , (cid:104) ξ (cid:105) ) will be ﬁt to the data.In principle, one could measure f π separately usingthe pion-axial current. However, measurements involv-ing heavy quarks are known to involve additional nor-malization factors, which have been approximated by El-Khadra, Kronfeld, and Mackenzie [37]. If we ﬁt f π fromthe hadronic tensor, any errors in this overall normaliza-tion factor will be absorbed into f π (a nuisance parameterfor us) rather than the parameter of interest, (cid:104) ξ (cid:105) .It should be noted that the operator product expan-sion is a statement about continuum physics. We will ap-proximate the left-hand side of (43) from an observabledeﬁned on the lattice and equate this to the right-handside deﬁned on the continuum, so our equation — andtherefore all parameters extracted from it — are only ac-curate to O ( a ). We will have to extrapolate away thislattice spacing dependence at the end of the calculation.To a good approximation, we can neglect the contribu-tion of (cid:10) ξ (cid:11) to the imaginary part of the hadronic tensor(this is a percent-level contribution) and use Im( U µν ) toextract f π and m Ψ . We can then use these ﬁt values asinputs to the determination of (cid:10) ξ (cid:11) using Re( U µν ). Fig. 5shows the results of this ﬁtting procedure at a = 0 .

06 fm. τ ( fm ) ℱ - [ I m ( U μ ν ) ] ( G e V ) τ ( fm ) ⅈ ℱ - [ R e ( U μ ν ) ] ( M e V ) FIG. 5. The real and imaginary parts of the hadronic ten-sor data (at a = 0 .

06 fm) ﬁtted to the continuum operatorproduct expansion.

V. LATTICE ARTIFACTS AND SYSTEMATICERRORS

Our results (including those shown in Fig. 5) are con-taminated with a variety of lattice artifacts and other sys-tematic errors, including excited-state contamination anddiscretization eﬀects. We are also working at an unphysi- cal pion mass of m π = 560 MeV in the quenched approx-imation. Finally, the operator product expansion we useis truncated, neglecting contributions both from highermoments (suppressed by the weight function W ( n )) andfrom higher-twist eﬀects that scale as Λ QCD / ˜ Q or m π / ˜ Q .Since this is a preliminary study, we will for now ignorethe eﬀects of quenching, unphysical m π , and higher-twisteﬀects. We also do not perform a chiral extrapolation,athough we note that others have studied the chiral be-haviour of the Mellin moments, and found the chiral de-pendence to be small. For example, Refs. [38, 39] showedthat at next-to-leading order in chiral perturbation the-ory, all possible non-analytic corrections to the matrixelements are contained in f π . This has also been previ-ously observed empirically, for example in Ref. [40]. Inthis analysis, we will focus on excited-state contamina-tion and lattice spacing dependence. A. Excited-State Contamination

For equation (34) to be valid, the pion creation oper-ator must be far from both current insertions. For ourmeasurements, τ m ≥ τ e , so excited states should be sup-pressed exponentially in τ e . Varying τ e at a variety ofcurrent separations τ = τ m − τ e shows substantial con-tamination at small τ e (see Fig. 6).In Fig. 7, we make a preliminary attempt to quantifythis error. A multi-state ﬁt to the ratio R µν used todeﬁne the hadronic tensor shows that excited-state con-tamination eﬀects are percent-level by τ e ∼ . τ e on the ﬁner ensembles, we estimate that excited-state contamination would be percent level there as well,which is smaller than the other errors (both statisticaland systematic) in our analysis. B. Continuum Extrapolation

Since our ﬁt variables f π and (cid:104) ξ (cid:105) were obtainedby matching lattice data to the continuum OPE, theyare contaminated with discretization errors. While thefermions used have a non-perturbatively set clover coeﬃ-cient, the axial current operators are unimproved, so theleading-order discretization artifacts are O ( a ).At each of the two heavy quark masses, the continuumlimit is found by linear extrapolation of the form (cid:104) ξ (cid:105) latt = (cid:104) ξ (cid:105) cont + Ca (44)as shown in Fig. 8. Since this study is only prelimi-nary, the range of lattice spacings used is currently small,leading to relatively large uncertainties in the continuumlimit: the more precise of our measurements gives (cid:104) ξ (cid:105) MS µ =2 GeV = 0 . ± .

07 (45) τ e = τ e = τ e = τ e = τ e = τ e = τ ( fm ) ℱ - [ I m ( U μ ν ) ] ( G e V ) FIG. 6. Excited state contamination is controlled by the sepa-ration between the pion source and the ﬁrst current insertion( τ e ). The inverse Fourier transform of the hadronic tensor( R µν ) is measured at a range of τ e values on the coarsest lat-tice ( a = 0 .

060 fm). Excited state contamination is clearlyvisible at τ e (cid:46) . τ e . τ e ( fm ) ℱ - [ I m ( U μ ν ) ] ( G e V ) FIG. 7. To make a quantitative estimate of excited-state con-tamination, we ﬁt the τ = 0 data shown in Fig. 6 with atwo-state ﬁt in τ e . Excited-state contamination is about 1%at τ e = 0 . Future work will increase statistics on the three existingensembles and add a fourth, ﬁner ensemble (at a = 0 . (cid:104) ξ (cid:105) in the continuum should be indepen-dent of the heavy quark mass up to higher twist eﬀects,which scale as inverse powers of ˜ Q and therefore vanish as m Ψ → ∞ . The level of precision of current measurementsis not suﬃcient to resolve these eﬀects, so the ﬁt valuesof (cid:104) ξ (cid:105) are compatible at the two heavy quark masses inthe continuum limit. VI. CONCLUSION

This preliminary study explores the potential of ex-tracting the moments of the pion light-cone distribution ( fm ) ξ ( fm ) ξ FIG. 8. The continuum extrapolation of (cid:104) ξ (cid:105) at heavy quarkbare masses of about 1.6 GeV (blue) and 2.5 GeV (earth),corresponding to renormalized masses of about 2.6 GeV and3.2 GeV, respectively. The uncertainties in the continuumlimit are large but should be reduced in future work. amplitude using the operator product expansion with aheavy valence intermediate quark. In the quenched ap-proximation at an unphysical pion mass, this methoddoes allow determination of (cid:104) ξ (cid:105) in the continuum limit,albeit with a large uncertainty (0 . ± . (cid:104) ξ (cid:105) using this method. ACKNOWLEDGEMENTS

CJL and RJP were supported by the Taiwanese MoSTGrant No. 109-2112-M-009-006-MY3 and MoST GrantNo. 109-2811-M-009-516. The work of IK is partiallysupported by the MEXT as “Program for Promoting Re-searches on the Supercomputer Fugaku” (Simulation forbasic science: from fundamental laws of particles to cre-ation of nuclei) and JICFuS. YZ is supported in part bythe U.S. Department of Energy, Oﬃce of Science, Of-ﬁce of Nuclear Physics, under grant Contract NumberDE-SC0012704 and within the framework of the TMDTopical Collaboration. AVG is supported by the U.S. De-partment of Energy, Oﬃce of Science, Oﬃce of NuclearPhysics under grant Contract Numbers DE-SC0011090and DE-SC0313021. WD is supported by the U.S. De-partment of Energy, Oﬃce of Science, Oﬃce of Nuclear Physics under grant Contract Numbers DE-SC0011090and DE-SC0018121. Numerical calculations are per-formed at the HPC facilities at National Chiao-Tung Uni-versity. We acknowledge support of computing resourcesfrom ASRock Rack Inc. Finally, we all acknowledge thehospitality of National Chiao Tung University and Mas-sachusetts Institute of Technology. [1] V. Chernyak, A. Zhitnitsky, and V. Serbo, JETP Lett. , 594 (1977).[2] G. Lepage and S. J. Brodsky, Phys. Rev. D , 2157(1980).[3] A. S. Kronfeld and D. M. Photiadis, Phys. Rev. D ,2939 (1985).[4] G. Martinelli and C. T. Sachrajda, Phys. Lett. B ,151 (1987).[5] M. Gockeler, R. Horsley, E.-M. Ilgenfritz, H. Perlt, P. E.Rakow, G. Schierholz, and A. Schiller, Phys. Rev. D ,5705 (1996), arXiv:hep-lat/9602029.[6] V. Braun, S. Collins, M. G¨ockeler, P. P´erez-Rubio,A. Sch¨afer, R. Schiel, and A. Sternbeck, Phys. Rev. D , 014504 (2015), arXiv:1503.03656 [hep-lat].[7] G. S. Bali, V. M. Braun, S. B¨urger, M. G¨ockeler, M. Gru-ber, F. Hutzler, P. Korcyl, A. Sch¨afer, A. Sternbeck, andP. Wein, JHEP , 065 (2019), arXiv:1903.08038 [hep-lat].[8] V. Braun and D. M¨uller, Eur. Phys. J. C , 349 (2008),arXiv:0709.1348 [hep-ph].[9] G. S. Bali et al. , Eur. Phys. J. C , 217 (2018),arXiv:1709.04325 [hep-lat].[10] Y.-Q. Ma and J.-W. Qiu, Phys. Rev. Lett. , 022003(2018), arXiv:1709.03018 [hep-ph].[11] R. S. Suﬁan, C. Egerer, J. Karpie, R. G. Edwards, B. Jo´o,Y.-Q. Ma, K. Orginos, J.-W. Qiu, and D. G. Richards,(2020), arXiv:2001.04960 [hep-lat].[12] Z. Davoudi and M. J. Savage, Phys. Rev. D , 054505(2012), arXiv:1204.4146 [hep-lat].[13] X. Ji, Phys. Rev. Lett. , 262002 (2013),arXiv:1305.1539 [hep-ph].[14] X. Ji, Sci. China Phys. Mech. Astron. , 1407 (2014),arXiv:1404.6680 [hep-ph].[15] X. Ji, Y.-S. Liu, Y. Liu, J.-H. Zhang, and Y. Zhao,(2020), arXiv:2004.03543 [hep-ph].[16] C. Alexandrou, K. Cichy, M. Constantinou, K. Jansen,A. Scapellato, and F. Steﬀens, Phys. Rev. Lett. ,112001 (2018), arXiv:1803.02685 [hep-lat].[17] H.-W. Lin, J.-W. Chen, X. Ji, L. Jin, R. Li, Y.-S. Liu,Y.-B. Yang, J.-H. Zhang, and Y. Zhao, Phys. Rev. Lett. , 242003 (2018), arXiv:1807.07431 [hep-lat].[18] J.-H. Zhang, J.-W. Chen, X. Ji, L. Jin, and H.-W. Lin,Phys. Rev. D , 094514 (2017), arXiv:1702.00008 [hep-lat].[19] J.-H. Zhang, L. Jin, H.-W. Lin, A. Sch¨afer, P. Sun, Y.-B.Yang, R. Zhang, Y. Zhao, and J.-W. Chen (LP3), Nucl.Phys. B , 429 (2019), arXiv:1712.10025 [hep-ph]. [20] R. Zhang, C. Honkala, H.-W. Lin, and J.-W. Chen,(2020), arXiv:2005.13955 [hep-lat].[21] A. Radyushkin, Phys. Rev. D , 034025 (2017),arXiv:1705.01488 [hep-ph].[22] K. Orginos, A. Radyushkin, J. Karpie, andS. Zafeiropoulos, Phys. Rev. D , 094503 (2017),arXiv:1706.05373 [hep-ph].[23] B. Jo´o, J. Karpie, K. Orginos, A. V. Radyushkin,D. G. Richards, and S. Zafeiropoulos, (2020),arXiv:2004.01687 [hep-lat].[24] K.-F. Liu and S.-J. Dong, Phys. Rev. Lett. , 1790(1994), arXiv:hep-ph/9306299.[25] J. Liang, T. Draper, K.-F. Liu, A. Rothkopf, and Y.-B. Yang (XQCD), Phys. Rev. D , 114503 (2020),arXiv:1906.05312 [hep-ph].[26] M. Oehm, C. Alexandrou, M. Constantinou, K. Jansen,G. Koutsou, B. Kostrzewa, F. Steﬀens, C. Urbach,and S. Zafeiropoulos, Phys. Rev. D , 014508 (2019),arXiv:1810.09743 [hep-lat].[27] K. Can et al. , (2020), arXiv:2007.01523 [hep-lat].[28] A. Chambers, R. Horsley, Y. Nakamura, H. Perlt,P. Rakow, G. Schierholz, A. Schiller, K. Somﬂeth,R. Young, and J. Zanotti, Phys. Rev. Lett. , 242001(2017), arXiv:1703.01153 [hep-lat].[29] W. Detmold and C. Lin, Phys. Rev. D , 014501 (2006),arXiv:hep-lat/0507007.[30] W. Detmold, I. Kanamori, C. D. Lin, S. Mon-dal, and Y. Zhao, PoS LATTICE2018 , 106 (2018),arXiv:1810.12194 [hep-lat].[31] H. Georgi and H. Politzer, Phys. Rev. Lett. , 1281(1976), [Erratum: Phys.Rev.Lett. 37, 68 (1976)].[32] O. Nachtmann, Nucl. Phys. B , 237 (1973).[33] W. Detmold, A. V. Grebe, I. Kanamori, C.-J. D. Lin,S. Mondal, R. J. Perry, and Y. Zhao, “The Heavy QuarkOperator Product Expansion I: Wilson Coeﬃcients,” Inpreparation.[34] R. G. L. C. Edwards and B. U. C. Jo’o, Nucl. Phys. Proc.Suppl. , 832 (2005), arXiv:hep-lat/0409003.[35] W. Detmold and M. Endres, Phys. Rev. D , 074507(2018), arXiv:1801.06132 [hep-lat].[36] R. G. Edwards, U. M. Heller, and T. R. Klassen, Phys.Rev. Lett. , 3448 (1998), arXiv:hep-lat/9711052v2.[37] A. El-Khadra, A. Kronfeld, and P. Mackenzie, Phys.Rev. D , 3933 (1997).[38] J.-W. Chen and I. W. Stewart, Phys. Rev. Lett. ,202001 (2004), arXiv:hep-ph/0311285.[39] J.-W. Chen, H.-M. Tsai, and K.-C. Weng, Phys. Rev. D , 054010 (2006), arXiv:hep-ph/0511036.[40] V. Braun et al. , Phys. Rev. D74