Parton physics from a heavy-quark operator product expansion: Lattice QCD calculation of the second moment of the pion distribution amplitude

The pion light-cone distribution amplitude (LCDA) is a central non-perturbative object of interest for high-energy exclusive processes in quantum chromodynamics. In this article, the second Mellin moment of the pion LCDA is determined as a proof-of-concept calculation for the first numerical implementation of the heavy-quark operator product expansion (HOPE) method. The resulting value for the second Mellin moment, determined in quenched QCD at a pion mass of $m_\pi=550$ MeV at a factorization scale of 2 GeV is $ \langle \xi^2 \rangle = 0.210 \pm 0.013\text{ (stat.)} \pm 0.034\text{ (sys.)}$. This result is compatible with those from previous determinations of this quantity.


I. INTRODUCTION
The pion light-cone distribution amplitude (LCDA) plays an important role in parton physics. It is central to a description of a range of exclusive processes in high energy quantum chromodynamics [1]. This LCDA, denoted as φ π , is defined via its relation to the matrix element for the transition between the vacuum and the (charged 1 ) pion state, where µ is the renormalization scale and W[−z, z] is a light-like Wilson line between light-cone coordinates −z to z (z 2 = 0). In the above equation, f π , p, and p µ are the decay constant, the three-momentum, and the fourmomentum of the pion. In the light-cone gauge, φ π (ξ, µ) can be interpreted as the probability amplitude for the pion to transition to the state of a quark and an antiquark that carry (1 + ξ)/2 and (1 − ξ)/2 fractions of the pion momentum, respectively [1]. This LCDA plays an important role in understanding exclusive processes in QCD [1,2], in addition to being a crucial ingredient for extracting information regarding the Cabibbo-Kobayashi-Maskawa matrix in flavour physics via hadronic decays of the B meson [3][4][5][6].
Since φ π (ξ, µ) encodes non-perturbative physics from strong interactions, it is natural to attempt to compute this quantity with lattice QCD (LQCD). However, LQCD calculations are normally carried out in Euclidean space because of the need to employ Monte-Carlo methods in evaluating the path integrals. Therefore it is challenging to adopt LQCD for investigating parton physics which involves dynamics on the light cone. This leads to the conventional LQCD approach of determining the Mellin moments for various parton distribution functions (PDFs) and LCDAs. In general these moments can be extracted by computing matrix elements of local operators that result from an operator product expansion (OPE) [7][8][9][10][11]. The analytic continuation between Euclidean and Minkowski spaces for these matrix elements is straightforward. In principle, knowledge of the relevant Mellin moments enables the construction of PDFs and LCDAs. Nevertheless, such a strategy has been limited by the possibility of having reliable LQCD results for only the first few moments because the breaking of O(4) Euclidean space-time symmetry by the lattice regularization requires the subtraction of power divergences in the renormalization of the relevant local operators [7]. These power divergences already appear in the computation of the first non-trivial (the second) Mellin moment for φ π (ξ, µ). They can be evaded by using a method proposed in Ref. [8], where one chooses particular combinations of the Lorentz indices in the local operator for computing the second moment. However, this method is not applicable for the extraction of higher moments, making it desirable to have other approaches with which to extract information of LCDAs and other light-cone quantities, including PDFs.
Alternative strategies for extracting information about LCDAs and PDFs using LQCD have been proposed in the last two decades [12][13][14][15][16][17][18][19][20], and their implementations are being intensively pursued, see Refs. [21][22][23][24] for reviews. To access information contained in higher Mellin moments while bypassing the above-mentioned complication in renormalizing the local operators, a generic character of these strategies is the computation for hadronic matrix elements of nonlocal operators. For instance, one much studied procedure, the quasi-PDF approach [17], features the calculation of matrix elements involving a space-like Wilson line. This work follows the method suggested in Ref. [14]. It relies on investigating the OPE analysis for hadronic amplitudes in Euclidean space with the insertion of two local quark bilinears that contain a fictitious, valence heavy quark. For this reason, this approach is termed the heavy-quark OPE (HOPE) method. The introduction of this heavy quark affords several advantages, such the removal of certain highertwist effects, the inclusion of an extra heavy scale to control the OPE, and simple properties of analytic continuation to Minkowski space. As pointed out in recent work [25], this method can be used to extract the ξ-dependence in φ π (ξ, µ) through a perturbative matching procedure. The current article presents its application to the numerical determination of the (second) Mellin and Gegenbauer moments 2 . Since this is the first numerical implementation of the HOPE strategy, this calculation bears the character of a feasibility study, which investigates the relatively well-known second Mellin moment. This allows a better quantification of the efficacy of the method. The success of the current work paves the way for further studies using this technique. For this reason, the quenched approximation and an unphysical pion mass of ∼ 550 MeV are used. Since the continuum extrapolation is an important ingredient in this approach, the calculations are performed at four choices of the lattice spacing, a, ranging from 0.04 fm to 0.08 fm.
Definitions relevant to this method and the conventions adopted in this article are declared in Sec. II. The HOPE method and its particular features for the current work are explained in Sec. III, while Sec. IV describes the details of the numerical implementation. The analysis of the data and the results are discussed in Sec. V. The conclusion of this work and the outlook in this research direction are then given in Sec. VI.

II. DEFINITIONS AND CONVENTIONS
Before proceeding to describe the numerical lattice determination for the second Mellin and Gegenbauer moments of φ π (ξ, µ), a review of the definition of these moments is in order. The discussion presented here also serves the purpose of setting the conventions and notation that will be used in subsequent calculations.
Conformal symmetry in QCD [29] implies that the pion LCDA can be conveniently studied using the Gegenbauer 2 Preliminary results from this calculation have been presented in contributions to the proceedings of recent lattice conferences [26][27][28].
To this order, their renormalization scale dependence is [30] where α s is the strong coupling, β 0 = 11 − 2N f /3 (N f being the number of flavours) is the coefficient of the leadingorder (LO) QCD β-function, and γ n is the anomalous dimension, Since γ n increases monotonically with n, one expects that a truncated version of the OPE in Eq. (2) can be a good approximation to φ π (ξ, µ) at large enough renormalization scales. In the regime where µ Λ QCD (with Λ QCD being the QCD dynamical scale), the pion LCDA is dominated by the zeroth Gegenbauer moment. Equations (4) and (5) lead to the asymptotic form the pion LCDA, where the normalization φ 0 = 1 has been imposed.
The Gegenbauer moments in Eqs. (2) and (3) can be expressed as linear combinations of the Mellin moments ξ n , which are defined as For instance, from Eqs. (2) and (7) it is straightforward to obtain In general, knowledge of ξ 0 , ξ 2 , ..., ξ n is equivalent to that of φ 0 , φ 2 , ..., φ n . These Mellin moments, ξ n , can be related to matrix elements of local, twist-two, operators, where the Lorentz indices are totally symmetrized, with and the traces are taken in all possible pairs amongst the Lorentz indices, µ 0 , µ 1 , . . . , µ n . As discussed in the last paragraph, from the leading-order result of QCD perturbation theory, it is natural to expect that knowledge of the first few Gegenbauer moments allows one to construct φ π (ξ, µ) reliably at sufficiently large µ. This also implies that π (ξ, µ), defined in Eq. (11), on the values of ξ 2 . The range of ξ 2 values shown in the plot covers typical results for this second Mellin moment at µ ∼ 2 GeV from modern lattice computations. They lead to consistency with a single-humped or double-humped structure of the DA, and more precise measurements would resolve this.
obtaining important information about the LCDA at µ Λ QCD is possible from the first few Mellin moments 3 . This point can be illustrated by investigating an extreme scenario where one truncates the Gegenbauer OPE in Eq. (2) at n = 2. The pion LCDA constructed with this truncated OPE is denoted φ (2) π (ξ, µ). Using Eq. (8), one obtains Figure 1 shows the result with φ (2) π (ξ, µ) at ξ 2 (µ) = 0.2, 0.25 and 0.3. Note that these are typical values for this second Mellin moment at µ ∼ 2 GeV from modern lattice computations [9-11, 32, 33]. This figure demonstrates that the shape of the pion LCDA can depend strongly upon ξ 2 (µ). Naturally, the inclusion of higher moments will likely reduce the sensitivity to the second moment, but nevertheless, this exercise shows that the second Mellin moment is a phenomenologically interesting quantity.

III. STRATEGY AND CORRELATION FUNCTIONS
To present the calculation for the second Mellin and Gegenbauer moments of φ π (ξ, µ), it is first necessary to describe the strategy and the correlators that have to be computed using LQCD. To extract moments for the pion LCDA employing the HOPE method, the hadronic amplitudes are computed, where with J µ A being the axial-vector current involving a light quark ψ and a fictitious valence heavy quark Ψ with mass m Ψ , The antisymmetrization in Eq. (12) is performed explicitly to reduce statistical noise, and the axial-axial correlator in Eq. (13) was empirically found to have less excited state contamination than the analogous correlator with two vector current insertions.

A. Relevant results from the HOPE strategy
Conformal symmetry in QCD [29] implies that it is natural to proceed with the Gegenbauer OPE for studying the twist-two contribution to the hadronic amplitude defined in Eq. (12). This OPE allows one to express the twist-two component of V [µν] (q, p) in terms of the same Gegenbauer moments, φ n (µ), defined in Eq. (3). These Gegenbauer moments do not mix under renormalization group (RG) evolution at one-loop [34,35]. At this order, generically the Gegenbauer OPE leads to [25] where µ is the renormalization scale, m Ψ is the mass of the fictitious valence heavy quark, and F n are coefficients that can be computed in QCD perturbation theory and can be expressed as functions of the kinematic variables with Q 2 = −q 2 . For simplicity, higher-twist contributions to V [µν] (q, p) will be discarded below in this section. They will be discussed in detail in Sec. V.
Employing results in Ref. [25], it can be demonstrated that As an example, the leading-order (tree-level) result of F n is independent of µ and gives, up to O(ω 2 ), where the relation between the leading two Mellin and Gegenbauer moments given in Eq. (8) is used. In this work, the goal is to compute the second Mellin and Gegenbauer moments of φ π (ξ, µ), working in the kinematic regime wherẽ ω 1, such that the OPE for V [µν] (p, q) can be truncated at the order ofω 2 .
Beyond the leading order (LO) in QCD perturbation theory, adopting Eqs. (15), (17) and (8), one obtains where the Wilson coefficients, C W (Q 2 , µ, m Ψ ), are linear combinations of F n,m (Q 2 , µ, m Ψ ) [25]. Since this work uses a relatively heavy pion (m π ∼ 550 MeV), it is beneficial to resum higher-twist target mass effects proportional to m π . The resummation prescription given in [36] is to replaceω n by ζ n C 2 n (η)/(n + 1)Q 2 , where ζ = p 2 q 2 /Q 2 , η = p · q/ p 2 q 2 , and C 2 n (η) is a Gegenbauer polynomial. In other words, Truncating at the order ofω 2 , where the explicit one-loop expressions for C (0) W (Q 2 , µ, m Ψ ) are given in Ref. [25]. Equation (21) is used in this analysis to extract ξ 2 . As described in Refs. [14,25], in addition to ξ 2 , f π and m Ψ are also fit parameters in the analysis procedure that will be presented in detail in Sec. V. Note that while the hadronic matrix element is renormalization scheme and scale independent, the factorization of this matrix element into short-distance Wilson coefficients and long-range Mellin moments are dependent on the renormalization scheme and scale. The calculation of the Wilson coefficients was determined in the MS scheme and thus the fitted heavy-quark masses and Mellin moments are directly extracted in this scheme.

B. The correlation functions
The power of the hadronic tensor lies in its amenability to lattice QCD calculations. The pion LCDA defined in Eq. (1) cannot be computed directly in Euclidean-space LQCD due to the light-like separation vector z. In contrast, the hadronic tensor V µν can be written in terms of quantities calculable on the lattice. Defining then the hadronic tensor is the Fourier transform of R µν in the temporal direction: Using lattice methods, one can compute two-point and three-point correlation functions and The three-point correlator is shown diagrammatically in Figure 2.
For 0 τ T , the 2-point correlator is saturated with the contribution of the lowest-lying hadronic state and can be written as which allows determination of the overlap factor 4 Z π (p) = 0|O π |π(p) and the pion energy E π (p).
Similarly, for 0 τ e , τ m T /2, the 3-point correlation function takes the form with p = p e + p m and τ = τ e − τ m , allowing one to extract from the two-and three-point correlators. From this and Eq. (21) and (23), one can extract the second Mellin moment on the lattice.

A. Lattice Action and O(a) Improvements
Order-a corrections to correlation functions arise from both the action and the interpolating operators [37][38][39]. Thus in order to remove these effects, one must in general improve both. The ensembles were generated with the standard Wilson gauge action where β = 6/g 2 is inverse coupling and P µν is the Wilson plaquette. The gauge action is automatically O(a) improved.
The Wilson fermion action, which is omitted from ensemble generation in the quenched approximation but is used in propagator construction, is given by where κ (f ) is the hopping parameter for flavour f [40]. The fermion action can be improved by addition of the clover term where σ µν = [γ µ , γ ν ] /2i andF µν is a discretized version of the field strength tensor corresponding to a sum over plaquettes. The clover coefficient c sw is taken from the non-perturbative tuning in Ref. [41].
The O(a)-improved, renormalized axial current operator is given by where Z

(0)
A is the axial-vector renormalization constant calculated in the chiral limit,m ij = (m i +m j )/2 is the average value of the masses of the two quark fields andb A , c A , c A are couplings which must be tuned to remove the O(a) corrections [42]. However, only b A is required for the O(a) improvement of the hadronic matrix element considered in this work.
This work studies the antisymmetric correlator R [µν] (τ, p, q), which may be obtained by taking the antisymmetric combination of Eq. (22). For this specific matrix element, it is possible to show that the term proportional to ac A vanishes by symmetry (see Appendix A 2 for details). The terms proportional to c A can also be shown to vanish; details of this argument are given in Appendix A 3. This means that both renormalization and O(a) improvement of the current are, for the purposes of this work, purely multiplicative effects. Z

(0)
A is given in Ref. [41] andb A is given in Ref. [42]. However, in the time-momentum representation analysis procedure described below, any multiplicative factors only affect the fitted pion decay constant f π and not the second moment ξ 2 , and thus the second moment is independent of any uncertainties in these O(a) improvement parameters. As a result, lattice artifacts in ξ 2 only enter at O(a 2 ) or higher order.

B. Lattice Parameters
Despite the O(a)-improvement scheme described above, the method in this work requires very fine lattice spacings because of the large mass used for the intermediate heavy quark propagator. With renormalized heavy quark masses ranging from about 2 GeV to 4.5 GeV, lattice spacings between 0.08 fm and 0.04 fm are needed to keep am Ψ 1. At larger values of am Ψ , lattice artifacts were uncontrolled and could not be reliably removed.
Due to critical slowing down, generating dynamical configurations at such fine spacings is expensive and beyond the scope of this preliminary study. As a result, this analysis uses gauge configurations generated in the quenched approximation following the multiscale procedure of Ref. [43]. The lattice spacings for the ensembles with L/a = 24, 32, and 48 had previously been determined in Ref. [43] using Wilson flow with a reference scale of w 0.4 = 0.193 fm [44]. This scale-setting procedure was repeated for the ensemble with L/a = 40. The lattice geometries were tuned to a constant physical volume of 1.92 fm, which was kept small to reduce computational costs.
Finite volume effects were suppressed by using light quark masses tuned to give m π ≈ 550 MeV so that m π L ≈ 5.3. The heavy quark masses were chosen to give approximately constant masses of the heavy-heavy pseudoscalar meson across the four lattices. Details on the lattices and quark masses used are listed in Table I. The required 2-and 3-point functions were generated using the software package Chroma with the QPhiX inverters [45,46].

C. Choice of Heavy-Quark Masses
The operator product expansion in Eq. (21) will require higher-twist corrections that scale like Λ QCD /Q or m π /Q. Unlike light-quark operator product expansions that rely on large momenta or small distances to suppress higher-twist effects [47], this work relies on the heavy intermediate mass for this suppression. With m Ψ q 4 |q|, higher-twist  effects will scale as Λ QCD /m Ψ or m π /m Ψ , so this analysis requires Λ QCD , m π m Ψ . Separately, lattice artifacts enter as powers of either aΛ QCD or am Ψ . Suppressing these will require am Ψ 1, so altogether it is required that With Λ QCD ∼ 250 MeV and m π ≈ 550 MeV, m Ψ > 1.8 GeV provides some suppression of higher twist effects. To fit the residual higher-twist effects away will require a range of heavy quark masses, and therefore this analysis will consider m Ψ as large as 4.5 GeV. The lattice spacing is small enough for the finest discretization to accommodate such a heavy mass while maintaining am Ψ < 1; coarser lattices out of necessity have a smaller range of heavy quark masses. Figure 4 shows the quark masses used on each of the four ensembles considered here. The plot shows the trade-off between discretization effects (which can depend on (amΨ) 2 ) and higher-twist effects (which scale as 1/mΨ). At fixed lattice spacing, one can decrease the higher twist effects at the cost of increasing discretization errors, and the available trade-offs at the four lattice spacings studied here are shown by the blue curves. The coloured points show the masses actually used in this study. The black dashed line at amΨ = 1.05 shows the cutoff beyond which discretization effects were no longer found to be well controlled, and the gray line at amΨ = 0.7 shows a more conservative threshold used for analysing systematic errors arising from lattice artifacts.

D. Choice of Kinematics
The heavy-quark OPE is given by where p is the momentum of the incoming pion and q is the difference in momenta between the two outgoing currents, and where the ellipsis represents the contributions of higher moments that are negligibly small in this analysis. The exact form of the higher-twist effects suppressed by Λ/Q is not known, but symmetries (see Appendix A 1) constrain it to be proportional to ε µνρλ q ρ p λ .
In order to enhance the contribution of the second moment, one would like its prefactor to be as large as possible. However,Q 2 must be large to suppress higher twist effects, and p is limited by noise that grows exponentially with the pion energy on the lattice. In this work,p ≡ p 2π/L was constrained to be one unit of momentum, which for the volumes used in this work corresponds to |p| = 640 MeV. At these kinematics, the second moment is a small contribution to the hadronic tensor. As such, it is desirable to isolate its effect from the much larger contribution of the zeroth moment. In this study, the axial current indices are fixed to be µ = 1, ν = 2, the prefactor on the right-hand side of Eq. (34) becomes If the kinematics are chosen such that p 3 = 0, then this prefactor is purely imaginary. At tree level, the entire contribution of the zeroth moment will be pure imaginary as well. However, p · q = iE π q 4 − p · q is generically complex (as long as p · q = 0), so the contribution of the second moment to the hadronic tensor will have nonzero real part. The effect of these special kinematics is shown in Fig. 5. This work met these criteria by choosinĝ in units of 2π/L, as well as all combinations ofp,q that are equivalent under lattice symmetries. With these choices, ξ 2 can be extracted as the leading contribution in the real part of V µν . (Note thatq = (p m −p e )/2, so it is quantized in half-integers rather than integers.) There are two caveats to this argument: Investigation of optimal kinematics for the HOPE procedure. Feasible lattice simulations are restricted to small lattice momentum, and thus for the HOPE method, smallω. Generally speaking, this results in moments beyond the leading zeroth moment being highly suppressed. This behaviour may be seen in the imaginary part of the amplitude where the variation of ξ 2 has minimal effect on the matrix element. By choosingp = (1, 0, 0) andq = (1/2, 0, 1), it is possible to show that the leading contribution to the real component of the amplitude arises from ξ 2 . Thus this choice of kinematics offers improved access to the second Mellin moment.
1. The Wilson coefficient of the zeroth moment C (0) W , while real at tree level, becomes complex at 1-loop order and also contributes to the real part of V µν . However, this contribution is suppressed by α s , and it is known analytically, so this small correction can be subtracted.
2. Higher-twist contributions may also be complex, but their contribution to Re[V µν ] must also contain powers of (p · q) 2 /(Q 2 ) 2 , just like the second moment contribution. They are further suppressed by Λ QCD /Q, so it is expected that they are smaller than the second moment contribution, particularly for largeQ 2 .
However, assuming the second Mellin moment is ξ 2 ∼ 0.25, both of the above contributions are subdominant at the special kinematics chosen here.

E. Computing Real and Imaginary Parts of V µν
To compute the hadronic tensor, values of R µν (τ ) ∝ C µν 3 (τ e , τ m = τ e + τ ) are needed at both positive and negative values of τ . In particular, using the fact that R µν (τ ) is pure imaginary 5 , the imaginary and real parts of V µν can be written in terms of symmetric and antisymmetric combinations of R µν (±τ ): At the kinematics of interest, the real part of V µν is about two orders of magnitude smaller than the imaginary part, so computing the difference in Eq. (37) requires a delicate cancellation between R µν (±τ ). The computation becomes more tractable if the two terms are highly correlated, as this increases the statistical power of the correlated difference. These correlations are substantially enhanced if values of C µν 3 for τ < 0 are obtained using the identity 6 C µν Then Eq. (37) and (38) can be written as Consequently, one can obtain both τ > 0 and τ < 0 at the same sets of current insertion times, which will enhance the correlations. A demonstration of this reduction in statistical error is shown in Figure 6.

F. Excited State Contamination and Choice of τe
The three point correlator C µν is computed by creating a pion source, propagating one of the quarks forward to τ e , creating a sequential source, and then tying together the sequential heavy-quark propagator and the other light quark propagator at the sink. Since τ e ≤ τ m is chosen in this work, excited state effects arise from the fact that the combination of states created by the pion interpolator have not fully relaxed to the ground state before τ e , so they are suppressed exponentially in τ e . Excited-state effects are reduced by using a Gaussian-smeared pion source [48] with smearing radius equal to the inverse pion mass (aw smear = {4.5, 6.0, 8.0, 9.0} for L/a = {24, 32, 40, 48}, respectively). With this smearing, numerical study on the L/a = 32 lattices showed that excited state contamination is estimated to be about 1% for a source-operator separation τ e of about 0.7 fm, as shown in Figure 7. 6 This identity can be proven by writing C µν 3 in terms of the quark propagators  Since τ e must be fixed at runtime, τ e ∼ 0.7 fm is chosen, leading to a ∼ 1% systematic error due to excited-state contamination. Excited state contamination in the 2-point function is better controlled since one does not need to choose the source-sink separation at runtime, and one can afford the very conservative fit range of [T /4, 3T /4] since the statistical errors on the 2-point function are smaller than those on the 3-point function (for an example of the goodness of fit, see Fig. 8).

V. ANALYSIS, RESULTS AND DISCUSSION
Extraction of the second moment from the 2-and 3-point correlation functions measured here is nontrivial due to signal contamination by both lattice artifacts and higher-twist effects. The extraction of Z π (p), E π (p) from fitting the 2-point correlation function C 2 (τ, p) at late τ and the construction of R µν (τ ) in Eq. (22) are relatively straightforward. However, comparing the 3-point data to the OPE of the hadronic amplitude can be done in multiple ways, which can lead to somewhat different systematics. As this is the first numerical study of the HOPE method, two analysis methods, called the time-momentum analysis and the momentum-space analysis, are performed. This enables a crosscheck of the results and ensures that they are robust against systematics in the analysis procedure. These analysis methods are as follows: 1. Time-momentum analysis (iii) Fit f π , m Ψ and ξ 2 to the hadronic amplitude in the continuum limit using the momentum-space, continuum HOPE formula presented in Sec. III A.
These alternative procedures are shown diagrammatically in Fig. 9. The following sections will detail both analysis strategies. Construct ratio from bare correlators: Perform Fourier transform to momentum space: V µν (p, q; a) = a τ e iτ q 4 R µν (τ, p, q; a) Extrapolate data to continuum: Fit continuum matrix element to HOPE formula: Fit lattice matrix element to inverse Fourier transform of HOPE formula: Extrapolate fixed lattice spacing parameters fπ(a, mΨ), and ξ 2 (a, mΨ) to continuum, mΨ → ∞ limit.
Momentum-space analysis Time-momentum analysis In order to avoid contamination with UV divergences near τ = 0, points with τ < 3a (grayed out in the plots) are excluded from the fit.

A. Time-Momentum Analysis
The ratio R µν (τ ; p, q) was constructed for 0 < τ ≤ τ max ≈ 0.6 fm. The statistical quality of the signal deteriorates with time, and large-τ data may be more susceptible to higher-twist contamination, motivating the cut at τ max . The symmetric and anti-symmetric components of R µν (τ ) are constructed as described in Section II-B.
An example fit to V µν for a single heavy quark mass at a single lattice spacing is shown in Fig. 10. At the chosen kinematics, the second moment provides a negligible contribution to the imaginary part of the hadronic tensor (see Fig. 5), so the fitting procedure can be split into two steps: one in which f π and m Ψ are fit to the imaginary part of V µν and a second step that consists of a single-parameter fit of ξ 2 to the real part of V µν , where f π and m Ψ are used as inputs. At values of τ comparable to the lattice spacing, uncontrolled discretization effects are to be expected. Additionally, if the two current insertions are close in space-time relative to the lattice spacing, they may mix with lower-dimensional operators and lead to UV divergences. Both of these effects suggest that small-τ data should be removed from the fits. Empirically, the χ 2 values for the fits to the various heavy-quark masses became reasonable if the τ ≤ 2a data are excluded, so all fits will only use data with τ ≥ 3a.
This fitting procedure compares lattice data to a continuum, twist-2 OPE. As a result, the extracted second moment ξ 2 (a, m Ψ ) will be contaminated by both lattice artifacts and higher-twist corrections. The lattice artifacts enter at O(a 2 ) (see Appendices A 2, A 3 for details), and by dimensional analysis, a 2 must be accompanied by two factors of a mass scale, either the typical momentum scale of Λ QCD or the heavy quark mass m Ψ , so there may be discretization effects proportional to a 2 , a 2 m Ψ , or a 2 m 2 Ψ . With am Ψ < 1.05, these terms were sufficient to describe lattice artifacts without need for additional O(a 3 ) terms. Higher-twist effects scale as powers of Λ QCD /Q or m π /Q, and Λ QCD ∼ m π in this analysis. The fitting procedure effectively integrates over the q 4 dependence, and m Ψ |q|, so the twist-3 contribution can be approximated by a Λ/m Ψ term. Therefore, to extract ξ 2 in the continuum limit without higher-twist contamination, ξ 2 (a, m Ψ ) is fit to the formula  Table I, plotted as a function of a 2 with the heavier masses at each lattice spacing displaced slightly to the right for visual clarity. The black star at a 2 = 0 represents the extrapolated value from the fit to Eq. (42).
spacing of a = 0.06 fm, the magnitudes of the various terms are where the renormalization scale for ξ 2 is taken to be µ = 2 GeV and the error bars are purely statistical. Neither the higher-twist nor the discretization effects can be neglected at the precision considered in this work. The fit result for ξ 2 is shown in Figure 11.

B. Estimation of Systematic Uncertainties for Time-Momentum Analysis
The analysis procedure described in the previous subsection contains several systematic errors. Excited state contamination in the 3-point function was estimated in Sec. IV F to be a ∼ 1% effect (and contamination in the 2-point function is much smaller). Finite-volume effects are expected to scale as 1 mπL e −mπL < 10 −3 and are negligible compared to both statistical and other systematic errors.
This work uses an unphysically heavy pion mass of m π ∼ 550 MeV. Previous studies [49] have indicated that ξ 2 at such a pion mass differs from its physical value by about 5%. Therefore, this is taken as a systematic effect arising from the unphysical pion mass.
Other systematic errors can be estimated by studying the effects of changing input parameters or varying the analysis procedure.
• The global fit described in Eq. (42) restricted the heavy quark masses to those satisfying am Ψ < 1.05. To test whether this cut is sufficient to exclude lattice artifacts of O(a 3 ) or higher, one could choose a more conservative cut, using only data satisfying am Ψ < 0.7. Refitting with this more limited data set results in a fit value of ξ 2 = 0.226 ± 0.043. Although these two results are compatible within one standard deviation, the difference between the central values (0.016) is taken as the estimate of the systematic uncertainty from the continuum extrapolation. • The global fit contains a Λ QCD /m Ψ term to account for the twist-3 contribution. In principle, higher-twist contributions are also present. To estimate such systematic effects, one could add a Λ QCD /m 2 Ψ term to the global fit in Eq. (42). This changes the fit result to ξ 2 = 0.185 ± 0.017 which has a central value differing from that of the primary procedure by 0.025. This is taken to be the systematic uncertainty from higher twist effects.
• As explained in Sec. V A, at small values of τ = τ m − τ e , the data are contaminated with uncontrolled lattice artifacts. The primary fit omits the τ /a = 0, 1, and 2 points, where such effects are the most significant and result in unacceptable χ 2 values in the fits. To analyse errors arising from the placement of this cut, one can exclude τ /a = 3 from the fits, which gives a modified result of ξ 2 = 0.208 ± 0.014 and therefore a small systematic uncertainty from the difference in central values of 0.002.
• The Wilson coefficients C W are calculated in perturbation theory, and in this analysis, they are only computed to 1-loop order. As an estimate of the magnitude of higher-loop corrections, one can perform this analysis at a larger renormalization scale of µ = 4 GeV and then run back to µ = 2 GeV using Eq. (4). Such a procedure results in ξ 2 (µ = 4 GeV) = 0.216 ± 0.012, which evolves to ξ 2 (µ = 2 GeV) = 0.218 ± 0.014, giving a systematic uncertainty of 0.008 from the change in central value.
The dominant sources of uncertainty are from the continuum and higher-twist extrapolations. In principle, both these extrapolations can be better controlled by including finer lattice spacings, which would also allow the inclusion of larger heavy-quark masses. However, computations at finer lattices are expensive and therefore beyond the scope of this preliminary work. The error from quenching is formally uncontrollable, although empirically it is a 10-20% effect in many calculations. To perform a precise comparison of this result to dynamical calculations would require redoing these calculations on dynamical ensembles.

C. Determination of fπ
The previous two subsections describe the determination of the second moment of the pion LCDA using the timemomentum analysis procedure. To check the validity of the HOPE strategy, it is worth noting that the pion decay constant f π is computed as a byproduct of this analysis. As is clear in the OPE formula, Eq. (21), f π is an overall normalization factor for the hadronic amplitude V µν .
One can extrapolate the f π values computed at various heavy quark masses on the four ensembles to the continuum using the same procedure as for the extrapolation of ξ 2 , giving a global fit value of 161 ± 2 MeV after removal of lattice discretization and higher-twist effects, where the error reflects statistical uncertainties only. It should be noted that this measurement suffers from not only the systematic errors mentioned in the previous subsection but also additional uncertainties from the normalization constants Z A andb A to which the second moment is completely insensitive in the time-momentum analysis. 7 On the other hand, f π can be directly measurable on the lattice via the axial-axial correlator at 0 τ T , where A 4 is the local unsmeared axial current ψγ 4 γ 5 ψ and the convention where f π ∼ 130 MeV at the physical pion mass is used. As written, Eq. (48) contains O(a) corrections, so the continuum extrapolation must include a term linear in a rather than in a 2 . A low-statistics (one source per configuration) computation of f π at the four lattices gives a continuum value of f π = 157 ± 6 MeV.
Despite the systematics that could affect the value extracted from the hadronic tensor measurement, the two determinations of f π are in good agreement (see Fig. 12). While f π is not directly relevant to the calculation, this serves as a useful cross-check of the validity of the operator product expansion and, more generally, of this calculational method.

D. Momentum-Space Analysis
A further check of the validity of the time-momentum representation method can be obtained by analysing the same data using a momentum-space analysis. The starting point for the momentum-space analysis is the time-momentum representation ratio R [µν] (τ, p, q; a) constructed in Eq. (28). The lattice-regularized data are converted to momentumspace via 7 Since this normalization factor is only used in the determination of fπ rather than in the computation of ξ 2 that is the main focus of the paper, this work uses the approximationb A am ij ≈ b A am ij , which is correct up to a mass-independent O(a) term [42]. The value of b A was taken to be a constant value of 1.25 across all lattice spacings, which is consistent with the values quoted in Ref. [42] at all lattice spacings measured in that work.
where τ max was taken to be approximately 1 fm = 5 GeV −1 . By Fourier transforming the tree-level HOPE equation, it is possible to show that the numerical data decay exponentially in τ as where E q = m 2 Ψ + q 2 . Since R [µν] (τ, p, q; a) exhibits exponential decay in τ with an exponent with magnitude greater than approximately 2 GeV, the truncation in the sum in Eq. (49) is expected to be well controlled. While the discrete Fourier transform formally only produces a discrete set of Fourier modes, in this work, interpolation between these Fourier modes is achieved by evaluating Eq. (49) for arbitrary q 4 . Note that the largest Fourier mode, q max 4 must be taken sufficiently small to remain below the Nyquist frequency which corresponds to requiring q max 4 < π/a. For the ensemble with the coarsest lattice spacing (a = 0.0813 fm) this results in the constraint that q max 4 < 7.5 GeV. In practice, data at momenta close to the Nyquist frequency may possess large lattice artifacts, and in this analysis q max 4 = 5 GeV is chosen.
While data at non-zero τ are guaranteed to have a well-defined continuum limit, τ = 0 data contain additional UV divergences arising from the mixing of the current-current operator with lower-dimensional operators. After performing the Fourier transform, this divergence will appear as an additive shift in the numerical data. Thus in order to ensure that the hadronic amplitude considered in this work has a well-defined continuum limit (after finite, multiplicative renormalization), a single subtraction is first performed at fixed lattice spacing: where q 4,sub is chosen to be q 4,sub = q max 4 = 5 GeV. This choice is informed by the desire to minimize the statistical noise introduced in this process. As a result of this subtraction, the matrix element may be expressed generically as where V [µν] sub (p, q) is the continuum hadronic matrix element. As emphasized previously, the simplicity of this continuum limit is one of the advantages of considering a current-current correlator. As argued in Sec. IV A, it is possible to show that the use of the the Sheikholeslami-Wohlert (clover) improved action leads to the removal of O(a) corrections for the matrix element studied here (see Appendices A 2, A 3 for details), so (2),sub (p, q) + O(a 3 ) .
It is important to note that the continuum extrapolation requires taking a → 0 along a line of constant physics. This is achieved by tuning the bare parameters to ensure certain physical quantities remain fixed as lattice spacing is varied, but there are some mistunings resulting from percent level inaccuracies in the tuning process. While such a mistuning appears as a relatively mild effect in the time-momentum representation analysis where it results in a variation of f π , the momentum-space approach is more sensitive to any such mistuning, since it affects the continuum extrapolation of the hadronic matrix element V [µν] (p, q). Evidence of this mistuning can be seen in Fig. 3. It is important to note that the tuning of the light-quark sector parameters to their physical values is an issue for all methods, and is not unique to this approach. One can reduce the systematic error associated with this mistuning by considering a ratio which is less sensitive to the volume and pion mass dependence. Examining the HOPE, the leading volume and pion mass dependent quantities are found in the prefactor (for µ = 1, ν = 2) for the special kinematics considered. In order to reduce the volume and pion mass dependence, the ratio V [12] sub (p, q)/(E π q 3 ) is formed and then continuum extrapolated. In contrast to the time-momentum representation analysis, a further cut must be made on the data included in the analysis. All data presented in the time-momentum representation analysis satisfy the constraint am Ψ < 1.05. This restriction is placed to ensure control over lattice artifacts. However, in addition to this constraint, the momentumspace analysis requires data at a sufficient number of lattice spacings for the continuum extrapolation to be performed. Re[V [12] sub (p, q)/E π q 3 ] (GeV Re[V [12] sub (p, q)/E π q 3 ] (GeV Re[V [12] sub (p, q)/E π q 3 ] (GeV In particular, the model for the continuum extrapolation contains two free parameters, and thus the analysis must be limited to the subset of data where the the heavy quark mass satisfies the constraint am Ψ < 1.05 for at least three lattice spacings. From Fig. 4, this criterion constrains this analysis to make use of only the lightest three heavy quark masses. The continuum extrapolations of the real and imaginary parts of the hadronic amplitude are shown in Figs. 13, 14.

E. Extraction of Second Moment
Having extrapolated the real and imaginary parts of the hadronic matrix element to the continuum, a global fit is performed to the order-α S HOPE formula given in Eq. Im[V [12] sub (p, q)/E π q 3 ] (GeV Im[V [12] sub (p, q)/E π q 3 ] (GeV from the continuum data is It is important to note that while the above parameterization of the twist-three piece is the most natural, other terms like m Ψ /Q 4 are also possible. The resulting fit is shown in Fig. 15. As a result of this global analysis, the first two moments of the pion LCDA are found to be f π = 0.173 ± 0.001 GeV, and ξ 2 (µ = 2 GeV 2 ) = 0.210 ± 0.013, where the statistical uncertainty is determined from a bootstrap analysis of the numerical data. Since this a global fit, the three heavy quark masses are also determined to be m Systematic uncertainties in ξ 2 for the momentum-space analysis may also be estimated following the same procedure as described in the time-momentum representation analysis. Since both analyses use the same lattice data they share some sources of systematic error. Thus the excited state contamination is taken to be a ∼ 1% effect, and as with the time-momentum representation analysis, finite-volume effects are taken to be negligible. Finally, as in the timemomentum representation analysis, the unphysically heavy pion mass is assumed to contribute a 5% systematic error. In addition to these shared sources of systematic uncertainty, the momentum-space method has several additional sources of systematic error. These arise from the difference in the order of steps of the analysis, and are discussed below: • Following the procedure employed above, O(a 3 ) effects are studied by making a more conservative cut on am Ψ .
In particular, the cut is chosen to be am Ψ < 0.7, which is consistent with the time-momentum representation. The resulting fit leads to a value for the second moment of ξ 2 = 0.222 ± 0.068. Taking the differences of central values leads to a systematic error of 0.012.
• The effects of higher-twist contributions may be studied by adding a twist-four ansatz to the continuum HOPE formula in Eq. (55). The form chosen is V [12] higher-twist (q, p, m Ψ )/(E π q 3 ) = V [12] (q, p, m Ψ )/(E π q 3 ) + mirroring the choice for the twist-three term. The resulting value for the second moment is ξ 2 = 0.245 ± 0.019. This results in a systematic uncertainty of 0.035.
• The higher-loop corrections to the Wilson coefficients are studied by repeating the above analysis at a renormalization scale of µ = 4 GeV, and then evolving back to µ = 2 GeV using the renormalization-group evolution of the Gegenbauer moments given by Eq. (4). This gives ξ 2 (µ = 4 GeV) = 0.236 ± 0.016. Running this to 2 GeV results in ξ 2 (µ = 2 GeV) = 0.239 ± 0.017. Taking the difference between the evolved ξ 2 and the original fitted ξ 2 (µ = 2 GeV) gives a systematic uncertainty of 0.029.
A breakdown of the sources of systematic error described here is given in Table III. This analysis of systematic errors leads to a final value for the second moment of ξ 2 (µ = 2 GeV) = 0.210 ± 0.013(stat) ± 0.044(syst). These two sources of error may be added in quadrature to obtain the final result ξ 2 (µ = 2 GeV) = 0.210 ± 0.046. Similarly to the time-momentum space analysis, the dominant source of systematic uncertainty arises from the higher-twist terms. This issue is made worse in momentum-space due to the additional cut on lattice data required for the continuum extrapolation of the hadronic amplitude.

F. Discussion of Results
The ratio R [µν] (τ, p, q; a) was analysed using two alternative approaches, termed the time-momentum representation (TMR) analysis and the momentum-space (Mom) analysis. The results of the second moment from these approaches are The central values and statistical errors are the same in both approaches. The agreement of central values is the result of statistical coincidence; with a different choice of fit parameters this extrapolated central value is expected to vary. The equivalence of the statistical error is relatively unsurprising, since both approaches mostly share the same raw lattice data.
As a cross-check, the pion decay constant f 2pt π = 0.158 ± 0.005 GeV was extracted from a conventional analysis of the axial-vector 2-point correlation function. This is to be compared with the HOPE-derived values f TMR π = 0.161 ± 0.002 (stat.) GeV and f Mom π = 0.173 ± 0.001 (stat.) GeV, with systematic uncertainties in f π . The systematic uncertainties in these determinations f π are likely comparable to the systematic uncertainty in ξ 2 (about 10-20%), or perhaps slightly larger due to the added uncertainty in the normalization factors.
Examining both procedures allows the study of the advantages and shortcomings of both approaches and serves as a further cross-check of the analysis. The above equations show that the time-momentum representation approach results in a smaller systematic error than that of the momentum-space analysis. While the systematic uncertainty incurred from the truncation of the twist expansion is the largest systematic error in both analyses, the additional cut placed on the data in the momentum-space analysis results in the removal of data with the heaviest heavy quark masses. Since higher-twist corrections are suppressed by factors of 1/Q ∼ 1/m Ψ , this results in less control over the higher-twist effects.
Given the above considerations, the central value is chosen to be the more precise time-momentum representation analysis value of This corresponds to a Gegenbauer moment of Most previous lattice calculations have used local operators to compute ξ 2 . In the quenched approximation, ξ 2 was previously computed to be 0.280 ± 0.051 at a renormalization scale of µ = 2.67 GeV [9]. Running this down to 2 GeV gives ξ 2 = 0.285 ± 0.054, which agrees with this quenched calculation, albeit with a large error bar.
More recent calculations with the local operator method have been performed in dynamical QCD, giving ξ 2 = 0.28 ± 0.02 [10] and ξ 2 = 0.235 ± 0.008 [50], both at µ = 2 GeV. A separate approach is to proceed via the quasidistribution amplitude (the distribution amplitude analogue of the quasi-PDF), which was used to give a result of ξ 2 = 0.244 ± 0.030 [51]. These results are compared to the second moment determined in this work in Fig. 16. Formally, the uncertainty from quenching cannot be controlled, so a precise comparison of the results in this work to these dynamical calculations is not possible. However, in practice, quenching errors are usually at the order of 10-20%, and the calculation presented here agrees with the dynamical results within the listed uncertainties combined with this approximate quenching uncertainty.

VI. CONCLUSION AND OUTLOOK
Factorization theorems in QCD imply that the LCDA of the pion is relevant to a variety of experimental processes. Since the LCDA is a non-perturbative object, it is a quantity of importance for LQCD calculations. While direct computation of φ π (ξ, µ) is impossible in a Euclidean field theory, a range of different theoretical approaches which allow one to indirectly study the LCDA have been proposed and pursued. These methods include direct calculation of the local matrix elements [7][8][9][10][11], factorization apporaches like the pseudo-DA approach [17], and a light-quark operator product expansion [15], have been used or proposed to this effect.
Knowledge of the first non-trivial Mellin moment of the pion ξ 2 provides an important constraint on the shape of the LCDA. Due to the one-loop running behaviour of the Gegenbauer moments, one expects that the second moment is especially important for the shape of the LCDA at large enough renormalization scale. This quantity has previously been studied with the conventional approach of calculating the relevant matrix element of the local twist-two operator. As a result, this quantity is relatively well known and therefore provides a good test of the validity and applicability of the new methods.
This paper presents the first numerical study of the HOPE method to extract the second Mellin moment of the pion LCDA. This approach utilizes a quenched fictitious heavy-quark species which enables more control over higher-twist effects. After a discussion of the numerical calculation of the hadronic matrix element, two alternative approaches were explored for extracting the second Mellin moment of the pion LCDA from the numerical data. These two approaches were termed the time-momentum analysis and the momentum-space analysis. Central to both analyses is the fact that the matrix element has a well-defined continuum limit after multiplicative renormalization.
In the time-momentum analysis, the O(α s ) formula of the HOPE is fit to the lattice data, and the resulting fit parameters are then extrapolated to the continuum. In the momentum-space analysis, the order of operations is reversed, and instead after Fourier transforming the lattice data, the correlators are extrapolated to the continuum before being compared with the O(α s ) continuum HOPE. Both analyses produce results in good agreement with each other, and with other calculations in the literature. Due to the order of operations, more lattice data are included in the time-momentum analysis. This leads to an improved estimate of the statistical error in the second Mellin moment. As a result, the final value of ξ 2 determined in this work is ξ 2 (µ = 2 GeV) = 0.210 ± 0.036 .
The uncertainty in this result is dominated by systematic effects, especially from higher-twist terms and the continuum extrapolation. From Table II, it is clear that reducing the systematic effect arising from higher-twist contributions is the most important task for improving these calculations. For this purpose, it could be helpful to adopt other heavy-quark formulations, such as that in Ref. [52], for future lattice numerical calculations.
The results shown here demonstrate the viability of the HOPE approach to determine moments of light-cone quantities with comparable statistical precision to results from other methods. This paves the way for further investigations of the pion LCDA, including dynamical studies of the second moment using the HOPE method and a determination of higher Mellin moments. Early numerical studies have commenced, and preliminary results for the fourth moment are discussed in Ref. [28]. The success of this approach for the LCDA also suggests that the HOPE method can be applied to the study of other light cone quantities. Key objects of interest are the kaon LCDA which would allow the study of Mellin moments of a system with non-zero strangeness and pion PDF and helicity PDF, for which the Wilson coefficients have already been calculated to one loop order [25].