Pion distribution amplitude from Euclidean correlation functions: Exploring universality and higher-twist effects

Building upon our recent study arXiv:1709.04325, we investigate the feasibility of calculating the pion distribution amplitude from suitably chosen Euclidean correlation functions at large momentum. We demonstrate in this work the advantage of analyzing several correlation functions simultaneously and extracting the pion distribution amplitude from a global fit. This approach also allows us to study higher-twist corrections, which are a major source of systematic error. Our result for the higher-twist parameter $\delta^\pi_2$ is in good agreement with estimates from QCD sum rules. Another novel element is the use of all-to-all propagators, calculated using stochastic estimators, which enables an additional volume average of the correlation functions, thereby reducing statistical errors.


I. INTRODUCTION
The lattice approach to QCD enables the computation of a multitude of hadronic parameters with high precision from first principles. Since the inception of this method, the list of quantities amenable to lattice simulation has been ever increasing. As the scientific focus moves on to ever larger classes of quark-gluon correlations the need for high precision lattice simulations to complement experimental data becomes ever more urgent. Hadronic contributions to the muon anomalous magnetic moment, which is on the verge of becoming a sensitive probe of physics beyond the Standard Model, constitute one such prominent example. In particular lattice calculations of the hadronic "light-by-light" scattering contribution [2,3] are set to become more precise than inferring this quantity from experimental measurements, see, e.g., [4,5] and references therein. Another venue which currently attracts a lot of attention is how lattice QCD may contribute to the determination of parton (i.e., quark and gluon) distributions in hadrons [6], which are scale-dependent nonperturbative quantities that enter the description of "hard" processes via QCD factorization theorems.
The possibility to calculate parton distributions from Euclidean correlation functions has been discussed for decades. For early work see, e.g., [7][8][9]. Recently, with the work by Ji [10] where it was strongly emphasized that nothing prevents one from accessing light-cone dynamics starting from Euclidean space, such approaches gained prominence. Several proposals exist, which differ in detail but share the same general strategy: The parton distributions are not calculated directly but extracted from suitable Euclidean correlation functions ("lattice * philipp.wein@physik.uni-regensburg.de cross sections" in the terminology of [11]; we prefer to use the term "Euclidean correlation functions" in this context because cross sections, in general, do not have a simple path-integral representation). After taking the continuum and other appropriate limits these can be expressed in terms of parton distributions in the framework of QCD factorization in continuum theory, in analogy to the extraction of parton distributions from fits to experimentally measured structure functions. In other words, the role of lattice QCD can be to provide a complementary set of observables from which parton distributions can be extracted, ideally, employing global fits combining lattice input with experimental data on hard reactions.
Such calculations are at an exploratory stage. At present the main task is to develop specific techniques that will eventually allow one to control all systematic errors. The pion light-cone distribution amplitude (DA) is the simplest parton function of this kind and offers itself as a laboratory where many of the relevant issues can be investigated. It also allows one to compare the strengths and weaknesses of the existing methods. Moreover, the pion DA is interesting in its own right as the main nonperturbative input to studies of hard exclusive reactions with energetic pions in the final state, e.g., the γ * γ → π transition form factor and weak B-meson decays B → π ν , B → ππ, etc.
In our recent publication [1], we have showcased the position space approach proposed in [12] and we will be using the same framework here. The new contribution of this work is to illustrate the advantages of considering several correlation functions simultaneously. Such a multi-channel approach not only leads to better statistics, but, most importantly, allows one to control and estimate higher twist corrections which otherwise lead to large systematic errors. The possibility to constrain higher twist corrections from the studies of lattice correlation functions is interesting within a much more gen- eral context and can have important applications. For the case at hand we find that the higher twist corrections extracted from lattice simulations agree very well with earlier estimates based on QCD sum rules and the phenomenology of hard exclusive reactions.
This article is organized as follows: Starting in Sec. II with a brief discussion of our approach and relating it to other methods used in the literature, we proceed in Sec. III to formulate the collinear factorization of correlation functions in position space, including one-loop results for the investigated current combinations. In Sec. IV we detail the methods used in our lattice computation. We present our results in Sec. V, before we conclude.

II. HEURISTIC DISCUSSION
Here we discuss a simple example for how the information on parton distributions at light-like separations can be extracted from the study of Euclidean correlation functions. We start from the pion transition form factor F πγγ (q 2 1 , q 2 2 ) of the reaction π 0 (p) → γ * (q 1 ) + γ * (q 2 ), which can be obtained from the matrix element of the product of two electromagnetic currents d 4 z e i(q1−q2)·z 2 ⟨0 T {j µ ( z 2 )j ν (− z 2 )} π 0 (p)⟩ = ie 2 µναβ q α 1 q β 2 F πγγ (q 2 1 , q 2 2 ) , where e is the electric charge and p = q 1 + q 2 is the pion momentum. The form factor F πγγ (q 2 1 , q 2 2 ) can be measured experimentally, at least in principle. If at least one of the photon virtualities is large, the form factor can also be calculated in QCD in terms of a single nonperturbative function describing the quark momentum fraction distribution u in the pion at small transverse separation, the pion DA. For the heuristic discussion in this section we consider the leading contribution shown in Fig. 1; the corrections are discussed in the next section. To this accuracy one obtains [13] F πγγ (q 2 1 , q 2 2 ) = − where F π ≃ 92 MeV is the pion decay constant. If the form factor is measured for a wide range of photon virtualities, the pion DA φ π (u) can be extracted from this relation (up to various higher order correction terms).
In practice such measurements are very difficult and experimental information is only available for kinematical situations where one virtuality is large and the second is close to zero [14,15], which is not sufficient to map out the complete u-dependence.
The integral ∫ d 4 z of (1) receives contributions from both space-like and time-like separations. Space-like correlation functions can readily be accessed in lattice simulations. However, addressing time-like distances is not at all straight-forward. The central observation at the root of the recent development is that time-like contributions are not needed (in the present context) as the complete information on the pion DA in principle is already contained in the space-like correlator.
Indeed, to the same accuracy as above, where is the pion DA in longitudinal position space, which is analogous to the Ioffe-time parton distribution in deep inelastic lepton-hadron scattering [16,17]. The correlation function in Eq. (3) can be calculated on the lattice for space-like separations z 2 < 0 and in principle arbitrarily large values of the scalar product p · z. In this way Φ π (p · z) can be directly measured [12]. Before going into details, we discuss the structure of the position space DA at a qualitative level to understand what kind of information can be obtained from such a measurement. Note that in the limit of exact isospin symmetry the equality φ π (u) = φ π (1 − u) holds. As a consequence Φ π (p · z) is a real function, Φ π (p · z) = Φ π (−p · z), with the normalization condition Φ π (0) = 1. The second derivative at the origin, Φ ′′ π (0), is related to the first nontrivial moment of φ π (u) which is usually denoted as ⟨ξ 2 ⟩ and referred to as the second Mellin moment in the DA literature: where ξ = 2u−1. This moment can be obtained on the lattice using conventional techniques [18][19][20][21] as the matrix element of a local operator that contains two covariant derivatives. Higher derivatives of Φ π at the origin are sensitive to higher moments. It has become standard to write the pion DA as a series expansion in orthogonal (Gegenbauer) polynomials where a π 0 = 1. Note that to one-loop accuracy the coefficients a π n (µ) do not mix under evolution of the scale µ. Moments of the DA can be written in terms of the coefficients in the Gegenbauer expansion, e.g., The corresponding expansion of the DA in position space is in terms of Bessel functions (conformal partial waves [12]) where The first few conformal partial waves F n (p · z 2), n = 0, 2, 4, are shown in Fig. 2. Since F n (ρ) ∼ ρ n for ρ → 0, the sum in (8) 3. Three models for the pion distribution amplitude (11) in momentum fraction (upper panel) and position space (lower panel). Note that Φπ(−p · z) = Φπ(p · z). extract the information on the pion DA beyond the first few moments, one has to include measurements at large p · z [12].
To view this from a somewhat different perspective, consider, for illustrative purposes, the one-parameter class of models at the reference scale µ 0 = 2 GeV. Three particular choices, cover a wide range of shapes that appear to be phenomenologically acceptable. These longitudinal momentum fraction space DAs and the corresponding position space DAs Φ π (p · z) are plotted in Fig. 3. The differences between the models increase with p · z. However, as demonstrated in Fig. 2, in the range accessible with present-day lattice calculations ( p · z ≲ 5) the differences are almost entirely due to the variation of the second Gegenbauer moment: a π 2 (µ 0 ) = 0.389, 0.146, 0 for the three above models, respectively.
So far we discussed the situation at tree level. Taking into account QCD corrections, the position space pion DA Φ π (p · z) in Eq. (3) will be substituted by a function of both scalar invariants, z 2 and p · z, of the form where φ (2) π ≡ φ π is the twist-two DA. The C VV n are coefficient functions that depend at most logarithmically on z 2 and are calculable in perturbation theory, while µ F is the factorization scale. We will tacitly assume using dimensional regularization and the modified minimal subtraction (MS) scheme. The superscript "VV" indicates the dependence of the coefficient functions on the choice of the correlation function used to define the position space DA -two vector currents for the present example Eq. (3). The leading (and higher) twist pion DAs are universal nonperturbative functions and independent of this choice. The power-suppressed O(z 2 ) correction terms correspond to higher twist pion DAs, like φ (4) π . The factorization scale µ F should be chosen similar in size to 2 √ −z 2 to prevent large logarithms from appearing in the coefficient functions.
The function Φ VV π (p · z, z 2 ), and/or similar correlation functions with different choices of currents, can be calculated on the lattice within certain ranges of the two arguments. Different strategies have been suggested how useful information can be extracted from such lattice data. In this work we follow the proposal of [12] as well as our work [1], carrying out the complete analysis in position space. We keep the distance between the currents sufficiently small to suppress higher twist effects and to enable the perturbative evaluation of the coefficient functions at the scale µ F ∼ 2 √ −z 2 ≥ 1 GeV, i.e., √ −z 2 ≲ 0.4 fm. At the same time √ −z 2 should be much larger than the lattice spacing, in this work a ≈ 0.071 fm, to tame discretization effects.
In the literature it has also been suggested to carry out a one-dimensional Fourier transformation of the lattice data in order to define new observables that are closer in spirit to the initial DA in longitudinal momentum fraction space, e.g., a quasi-distribution [22][23][24][25][26][27][28][29][30][31][32][33] or a pseudo-distribution [34][35][36] φ ps Both expressions are designed in such a way that, to leading twist accuracy, they reproduce the pion DA at tree level. In the existing calculations which employ the above methods two spatially separated quark fields are connected with a Wilson line. Equivalently, this construction can be viewed as a correlation function involving two bilinear currents with an auxiliary "heavy" quark field [37][38][39] rather than the light quark field we use in Eq. (3). Apart from employing a Wilson line [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36]40] or an auxiliary light quark propagator [1,12,41], other obvious choices for connecting the two positions include a scalar propagator [8,9], a heavy quark propagator [42], or just employing Coulomb gauge [43]. Another technical difference of the quasi-distribution work relative to our approach is the use of the largemomentum factorization scheme at an intermediate step ("Large momentum effective theory" (LaMET) [44,45]) to emphasize that, for a large pion momentum and at a fixed quark momentum fraction, large distance (i.e., higher twist) contributions are suppressed.

A. Collinear factorization in position space
A general approach to implement collinear factorization of QCD amplitudes in position space is provided by the light-ray operator product expansion (OPE) [46][47][48][49][50][51]. For a generic current product one writes where while Z k are the renormalization factors for the currents, Π µ F l.t. [. . .] is the leading twist projection operator, C 12 is the coefficient function, and the ellipses stand for higher twist contributions. For simplicity we disregard the flavour structure, showing only the contribution of flavour-nonsinglet axialvector operators that will be important for this work. The corresponding expression for the product of quark and antiquark fields connected by a Wilson line is exactly the same, with Z 1 , Z 2 substituted by the quark field renormalization factors in space-like axial gauge.
The leading twist projection of a nonlocal quarkantiquark operator is defined as the generating function of renormalized local leading twist operators (traceless and symmetrized over all indices), e.g., where z = (z 1 + z 2 ) 2 and Here and below we indicate trace subtraction and symmetrization by enclosing the involved Lorentz indices in parentheses, for example O (µν) = 1 2 (O µν +O νµ )− 1 4 g µν O λ λ . The light-ray OPE differs from the usual Wilson expansion in local operators by imposing a different power counting. In the latter case one assumes that the distance between the currents is small, z 12 ∼ ηΛ −1 QCD with η → 0, and the operator matrix elements are of order unity in this limit, ⟨O n,k µ1...µn ⟩ ∼ Λ n QCD . In this case only a finite number of local operators on the r.h.s. of Eq. (17) has to be kept, and also the higher twist operators must be added progressing to higher powers of η: The relevant expansion parameter is the operator dimension, not the twist. The light-ray OPE assumes instead that ⟨O n,k µ1...µn ⟩ ∼ η −n Λ n QCD so that z µ1 12 . . . z µn 12 ⟨O n,k µ1...µn ⟩ = O(1) and in this case the series (17) must be resummed to all orders. Such a situation occurs if the hadron has large momentum, p = O(η −1 ) and hence p · z 12 = O(1), since for generic hadronic matrix elements where the reduced matrix element ⟨⟨O n,k ⟩⟩ = O(1).
Higher twist operators of the same dimension have smaller spin (by definition). As a consequence their matrix elements involve less powers of the large momentum and are suppressed. At the amplitude level, expanding in powers of the large momentum corresponds to the classification in terms of the so-called collinear twist, see, e.g., [51,52]. Note that the above power counting is applicable both in Minkowski and Euclidean space. In Minkowski space one can employ a reference frame where all components of the momentum are small and simultaneously the separation between the currents is almost light-like, z µ = O(1), z 2 = O(η 2 ) → 0. In this way the usual interpretation as the light-cone expansion arises.
The light-ray OPE provides a technique to deal with leading twist projected operators (17) as a whole, avoiding the local expansion. These can be viewed as analytic operator functions of the separation between the currents (all short-distance and light-cone singularities are subtracted) and satisfy the equation [48] Explicit expressions for the projection operator Π µ F l.t. can be found in [48,[51][52][53]. This technique combined with the background field method has proven to be very efficient and has found many applications, e.g., in light-cone sum rules [54] for the calculation of higher twist contributions, and for the derivation of the evolution equations for off-forward parton distributions [55,56].
Hadronic matrix elements of the operator (17) define leading twist parton distributions. Specializing to our case, the pion DA is defined via where [50] The second term in the last line is the (twist-four) pion mass correction, which is analogous to the Nachtmann target mass correction in deep-inelastic scattering.

B. Choice of currents and one-loop results
In this work we perform a lattice study of the set of correlation functions where the currents J X ≡q Γ X u are defined as and contain an up quark u and an auxiliary quark field q.
In this study we assume that the auxiliary quark has different flavour, q ≠ u, d, but same mass (m q = m u ) as the light quarks. For convenience and better readability we invoke the obvious notation T µν VA ≡ T V µ A ν , etc. We do not consider the correlation functions of S(P) with V(A) currents because they are dominated by (chiral odd) higher twist DAs. For the correlators with two Lorentz indices the most general invariant decomposition reads where the prefactors are by construction invariant under rescaling of z and all invariant functions T XY ≡ T XY (p·z, z 2 ) have the same mass dimension. The Lorentz decomposition for T µν AV is obtained from the one for T µν VA by replacing V ↔ A. One can show that T VA ≡ T (1) VA is the only invariant function in the VA correlator that receives contributions from the leading twist DA at leading order in perturbation theory, so that we only consider this structure in what follows. The projection needed to isolate it is specified in Appendix A. Finally, T SP and T PS are scalar functions which we write below as T SP and T PS , respectively, to unify the notation.
Separating a common overall prefactor, it is convenient to write the correlation functions in the form where to tree level accuracy and neglecting higher twist is the position space pion DA of Eq. (4). We further separate the leading twist (LT) contribution from the higher twist (HT) part: where the higher twist contributions are of O(z 2 ), cf. Eq. (12). The calculation of the one-loop, i.e., O(α s ), correction at leading twist is relatively straight-forward.
Using the Gegenbauer expansion of the pion DA, Eqs. (6) and (8), the result can be written as Setting the renormalization and factorization scales to the same value µ = µ F , we obtain, to O(α s ) accuracy, where α s = α s (µ), the functions F n are defined in Eq. (9), ρ = p · z 2, C F = 4 3 , η = 1 + 2γ E + ln(−z 2 µ 2 4),s = 1 − s. In the following, we will choose µ ≡ 2 √ −z 2 . The plus prescription is defined as usual, The sum in (28) converges very rapidly since cf. Fig. 2, so that for moderate p · z only the first few Gegenbauer moments give a sizeable contribution [12]. The two-loop corrections O(α 2 s ) are known for the VV correlator [12] but not for other cases, so that we do not include them in this study.
For reasons that will be explained in Section IV B, we choose to analyse linear combinations of correlation functions: PS + SP, VA + AV, and VV + AA, for which the leading quark mass correction originating from the chiral odd part of the auxiliary quark propagator does not contribute. In this sum, e.g., the imaginary parts of the SP and PS correlators cancel each other and drop out. For the VA + AV case the terms linear in the quark mass m q drop out completely after applying the projection onto the invariant function of interest as described in Appendix A. Note that the entire difference between f VV and f AA is due to this quark mass correction, which is converted into an O(m 2 π ) term using the axial Ward identity and m q = m u .
As observed already in [12], the higher twist correction in the VV channel has opposite sign compared to the leading twist term. We find a similar behaviour for the AV channel. In contrast, the twist-four correction for the SP correlation function has the same sign as the leading twist contribution. Numerically, the higher twist corrections turn out to be approximately of the same size as the leading perturbative correction at √ −z 2 2 ∼ 0.2 fm ≃ 1 GeV −1 , and become gradually less important for smaller distances. At the lowest scale considered in this study, 1 GeV, and at p·z = 0 the combined one-loop and higher twist correction yields approximately −40%, −20%, and +50% for the VV, VA, and SP channels, respectively. We will find that these estimates are strongly supported by our lattice data, cf. Section V.

IV. LATTICE CALCULATION
We employ the same gauge ensemble as in [1] (ensemble IV of [60], generated by QCDSF and RQCD), which allows a direct comparison between the sequential source method [61] (used in [1]) and the stochastic method (applied in this work) for the scalar-pseudoscalar channel. We employ the Wilson gluon action with two mass-degenerate flavours of nonperturbatively order a improved Sheikholeslami-Wohlert [62] (= clover) Wilson fermions. The lattice consists of 32 3 × 64 points with periodic boundary conditions (antiperiodic in time for the fermion fields). The inverse gauge coupling parameter reads β = 5.29 and the hopping parameter value is κ = 0.13632. This corresponds to the lattice spacing a ≈ 0.071 fm = (2.76 GeV) −1 [63] and a pion mass m π = 0.10675(59) a ≈ 295 MeV [64]. To reduce autocorrelations we have used a binsize N bin = 20 for the N conf = 2000 configurations we have analysed, cf. Table I. In order to improve the overlap between the interpolating current at the source and the pion state at large momentum, we employ the momentum smearing technique of [65] (see also [21]) with APE smeared spatial gauge links [66].
The operator renormalization is performed as described in [67]: The renormalization factors are calculated nonperturbatively within the RI ′ -MOM scheme [68,69] (along with a subtraction of lattice artifacts in oneloop lattice perturbation theory). These are then converted to the MS scheme using three-loop (continuum) perturbation theory. The corresponding factors for the MS scale µ = 2 GeV can be found in Table III of [64]. To be consistent, we employ the N f = 2 specific running of α s in all perturbative calculations. To this end, we combine the results of [63] and [70] to obtain a value of α s at 1000 a ≈ 2.76 TeV. From there we evolve it downwards using five-loop running [71]. The pseudoscalar and scalar currents are evolved to other scales using the fourloop mass anomalous dimension, which is consistent with the order used in [67]. The numerical values of the N f = 2 coefficients are summarized, e.g., in [72], which also includes the five-loop calculation.
Disconnected quark line diagrams have proven to be notoriously challenging in lattice simulations. They can be avoided by implementing an appropriate flavour structure of our currents: We pretend that the auxiliary quark field q of Eq. (23) is of a different flavour but shares its mass with the light quarks: m q = m u = m d . This does not present any limitation as the perturbative matching is carried out using the same conventions.
In the following sections we use boldface letters for the space components of the distance and momentum, (z µ ) = (0, z) and (p µ ) = (E p , p). In the actual lattice calculation we evaluate the three-point functions using currents positioned at z, relative to our origin 0. These are "shifted" afterwards to the symmetric locations as in Eq. (23) by multiplication with the appropriate phase.

A. Stochastic estimation of correlation functions
We wish to compute the correlation functions Eq. (23). The corresponding three-point functions for different Γ structures are depicted in Fig. 4, where the straight lines correspond to quark propagators. The momentumsmeared, momentum-projected pion source is located at the Euclidean time slice 0. Translational invariance of the correlation function implies that z and 0 can be shifted to the positions z 2 and −z 2, respectively, by multiplication with the phase e ip·z 2 . Previously, in [1] we computed propagators, starting from a point source at the position (t, 0), smeared the resulting propagator at the time slice t = 0 and computed a sequential propagator [61] from there. This propagator and the original propagator were then contracted with the Γ X structure at (t, z), making use of γ 5 -Hermiticity of the propagator G, i.e., G xy = γ 5 G † yx γ 5 . In order to increase the statistics, ideally one would average over different spatial positions y of the current J Y , placing J X at positions y+z, keeping the relative distance vector fixed. It turns out that this is indeed possible, introducing stochastic propagators [73], albeit at the cost of additional (but small) stochastic noise. Therefore, our new approach is to start from a momentum smeared pion source at t = t src and compute stochastic forward propagators from there, using the "one-endtrick" [74]. As with the sequential method adopted by us previously, the external momentum is fixed at the source. Since we need to keep the distance z between the local currents J X and J Y at the sink fixed, volume averaging would not be possible if we created a sequential propagator at the sink. Instead, we use a second stochastic volume source at t sink to connect these two currents. In order to reduce the associated stochastic noise we utilize the hopping parameter expansion in the way suggested in [75,76] (see also [77,78] for related work) to block out the dominant short-distance noise contributions when connecting the two currents with a stochastic propagator. Below we describe our implementation in detail.
We define the momentum smearing operator F p = Φ n (ζp) with n smearing iterations (n = 200 in our cal-culation). This is diagonal in spin and constructed on the time slice t src , iteratively applying the operation where U x,j is an APE smeared [66] spatial gauge link connecting the lattice points (t src , x) and (t src , x + a), for details see [21,65]. In practice this smearing is implemented by multiplying the spatial connectors within the time slice in question by the appropriate phases, U x,j ↦ e −iakj U x,j , where k = ζp. We choose ζ = 0.8 and = 0.25. We write the Wilson-Dirac operator as We also define the time, spin and colour diagonal momentum projection operator ϕ p with the components where i and α denote the colour and spin indices, respectively. We then solve where χ andχ (as well as ξ) are Dirac vectors with colour, spin and spacetime components. Above, we have suppressed these indices for enhanced readability. In our conventions ξ, χ andχ (as well as η and s, which will be introduced below) are dimensionless. We define a momentum smeared interpolator that, when applied to the vacuum, will create states with the quantum numbers of a π 0 carrying the (spatial) momentum p, and a local isovector current J v = ūΓu −dΓd 2 with an arbitrary Dirac structure Γ. In the following we assume, for the sake of readability, that all sources have been shifted to t src = 0 (exploiting translational invariance) and denote the source-sink distance as t. We can now obtain the average over the spatial volume V 3 of smearedlocal two-point functions (39) as an inner product over colour, spin and (threedimensional) space: Here we have suppressed all unnecessary indices. The disconnected contractions drop out since we have exact isospin symmetry. The minus sign in the first line is due to fermion anti-commutation. Within the scalar product (A, B) = A † B we sum over all indices that are not displayed on either side, in this case spatial position, colour and spin. In the second line we used a −4 D −1 0t = q 0qt for q ∈ {u, d}. In the last step we made use of the orthonormality of the noise vectors when aver- where X, Y represent multi-indices, and of the γ 5 -Hermiticity of the propagator.
To construct the desired three-point function, we use additional spin-partitioned [79] (also referred to as spinexplicit or "diluted" in the literature) stochastic sources η (k,α) , k = 1, . . . , n st , and α = 1, . . . , 4, with the components η for each value of k and α to obtain s (k,α) . The lattice propagator G = a −4 D −1 from (t, y) to (t, y + z) can now be estimated as up to a stochastic error that decreases ∝ 1 √ N conf n st , where N conf is the number of gauge configurations and n st = 10 is the number of spin-partitioned stochastic sources.
The operator H in Eq. (35) only couples nearest neighbours for the action we use. Employing the geometric series where we can split up the propagator into the first sum in Eq. (45) that does not contribute at distances z (and distances that are separated by a larger number of hops) and a part that contributes. In the stochastic estimation the first part still adds to the noise. This undesirable effect can be removed, left multiplying the solution with H m(z) [75,76]. Looping over momenta and times, we define temporary scalar fields for m ≤ 10 and all currents of interest. These fields implicitly depend on p and t. Also note that the solutions χ andχ of Eq. (37) depend on the momentum p. The three-point correlation functions can now readily be obtained by replacing J v (t, y) ↦ J † X (t, z 2)J Y (t, −z 2) in Eq. (39) (cf. Fig. 4). The result reads where the value of m ≤ m(z) used within the set of precomputed fields K (m,k,α) X is selected as large as possible for each distance. Note that we have already shifted the above correlation function to the symmetric position. In our study we limit ourselves to the range z i ≤ 5a.
With the previous sequential source method, first one propagator (12 solves) had to be computed. Then for each additional momentum and time separation, two smearing operations were required as well as an additional propagator (12 solves). In our implementation of the new method we vary the distance between the pion source and the sink by changing the time slice where the pion source is placed, enabling us to reuse the stochastic solutions H m s (k,α) and sources η (k,α) of Eqs. pion source is seeded, necessitating only two additional smearing operations and two additional solves. In total, not even taking into account that there is an additional gain from the two possibilities of connecting the valence quark propagators with the stochastic propagator of the auxiliary field (giving us for each momentum p the momentum −p almost for free), the new method does not only allow for a volume average, thereby reducing statistical errors, but turns out to be cheaper by about a factor of two in terms of the total computer time.
The three-point function C 3pt XY admits the same spectral decomposition, Eq. (41), as the two-point function C 2pt . Just the matrix element needs to be replaced: ⟨0 A v 0 (0) π 0 (p)⟩ ↦ ⟨0 J † X (0, z 2)J Y (0, −z 2) π 0 (p)⟩. The overlap factor Z π (p) and the exponential decay cancel when taking the ratio of these two functions. Therefore, in the limit of large Euclidean times, where excited state contributions are exponentially suppressed, the ratio can be related to the matrix element of interest: where Z X is the renormalization factor of the local current J X with respect to the MS scheme [67]. For the scalar and the pseudoscalar currents the renormalization factors acquire a scale dependence due to their anomalous dimension.

B. Reducing discretization effects
In the continuum, the chiral even part (∝ z) of the propagator connecting the two local currents gives the most important contribution, while the chiral odd part (∝ 1) is suppressed by a factor m √ −z 2 and, thus, can be set to zero in a first approximation. However, with Wilson fermions the situation is completely different. We find that the contribution from the chiral odd part, which suppresses the doublers and breaks chiral symmetry, can be of the same order of magnitude as the leading contribution, cf. Fig. 5. The "jumping" of the points nicely demonstrates the strong dependence of the lattice artefacts on the chosen direction. In particular the points along the axes (e.g., (1, 0, 0)), corresponding to the crosses in Fig. 5, exhibit the largest discretization effects, while the points along the diagonal (e.g., (1, 1, 1)) are much better behaved. The large contribution of the chiral odd part of the propagator is a peculiarity of using Wilson fermions, while large discretization effects at short distances are probably a general feature of all position space methods.
The appearance of large contributions from the chiral odd part of the propagator would lead to huge lattice artefacts in the correlator. Therefore, we construct linear combinations of the correlation functions defined in Eq. (23) where the chiral odd part of the propagator drops out to leading order in perturbation theory: For the scalar-pseudoscalar correlator with a pion this is equivalent to taking the real part (cf. [1]). The discretization effects in the chiral even part of the propagator (these correspond to the blue points in Fig. 5) are addressed as follows: We discard data points where the free field discretization effect exceeds 10%. This cut mainly excludes very short distances ( z ≲ 2a) and directions close to the lattice axes. For the remaining data points we define a correction factor c corr (z) such that the corrected propagator satisfies the condition tr zG corr latt (z) where the trace runs over spin and colour indices. To zeroth order accuracy in α s (where G latt = G free latt is the free propagator) this leads to This corresponds to multiplying the blue data points of Fig. 5 by factors so that in the non-interacting case the continuum result is retrieved. One should note that this procedure can only tame distance dependent discretization effects. However, there are also momentum dependent discretization effects, which are not taken into account. It is therefore no surprise that we still find particularly large discretization effects for the high momentum data at small distances. We have therefore decided to include only data points with z ≥ 3a ≈ 0.21 fm, which, setting the scale to µ = 2 z , corresponds to µ ≲ 1.84 GeV. Finally, we remark that the pseudoscalar and scalar currents are (up to small mass-dependent effects) automatically order a improved. In principle we could also have order a improved the axialvector and the vector currents. However, the improvement of γ µ γ 5 and of γ µ would have required us to compute three-point functions with two currents situated at non-equal times (as well as a tensor current in the latter case).

A. Parameter choices and first data survey
Our analysis includes 6 different pion momenta (12, if one counts ±p separately) with absolute values up to p = 2.03 GeV, cf. Table I. For the largest momentum, we have analysed two different directions to increase statistics. Reaching such a large hadron momentum is quite challenging and was achieved by the combination of the momentum smearing technique, which enhances the overlap of the interpolating current with hadrons at large momenta, and the use of stochastic estimators described in Section IV A, which allows us to take a volume average at the cost of additional stochastic noise. The latter trade-off turns out to be very advantageous and yields a significant reduction of the statistical errors compared to the sequential source method used in [1].
Since the lattice data are analysed using QCD factorization in the continuum, we are bound to use sufficiently small separations between the currents to ensure that the coefficient functions are perturbatively calculable. Together with the requirement of controllable discretization effects (see the previous section), this leaves us with the relatively narrow range of possible distances 0.21 fm ≲ z ≲ 0.39 fm, or 3a ≤ z ≤ 5.5a in units of the lattice spacing a ≃ 0.07 fm. Since the direction of z is arbitrary, this constraint still allows for a large data set with ten different values for z . First, however, we should check if the ratios of three-point over two-point functions (50) approach their asymptotic limits. We demonstrate this for the combination (T VV +T AA ) 2 for different momenta and distances in Fig. 6. Clearly, the momentum smearing was extremely successful in removing excited state contributions. Moreover, these seem to affect the two-point function in a similar way as the three-point functions, enabling additional cancellations to take place. The other channels exhibit a very similar behaviour so that we can confidently fit to extended plateaus.
Next, in Fig. 7 we compare our results at two typical distances for two different channels with the expectation obtained using the second Gegenbauer coefficient a π 2 (2 GeV) = 0.1364 that has been determined in [20] with the moment method. The leading twist position space DA (central solid line) is universal for all channels. The dashed lines include our one-loop perturbative corrections while the solid lines also include higher twist effects using the QCD sum rule estimate δ π 2 (2 GeV) = 0.17 GeV 2 [57][58][59]. Unsurprisingly, towards the larger distance z both correction terms become more 1.8 7. Data for the position space DA at two distances compared with expectations obtained using the second Gegenbauer coefficient a π 2 (2 GeV) = 0.1364 determined in [20] with the moment method. The central solid curve corresponds to the (channel independent) tree-level result at leading twist. The dashed lines include one-loop perturbative corrections for the two channels and the outer solid lines also include higher twist contributions (obtained using the QCD sum rule estimate δ π 2 (2 GeV) = 0.17 GeV 2 for the higher-twist normalization constant). The upper data (green) are SP + PS and the lower data (blue) VV + AA.
significant. Sign and magnitude of the predicted splitting are in good agreement with our data. However, there are quantitative differences: Our data still show residual discretization effects, the models for the leading twist and higher twist DAs may not be correct, and there will be two-loop perturbative corrections as well. For the distances shown, the corrections to the leading order leading twist DA are about 25% in size, while even at p · z = 4, the differences between the models plotted in Fig. 3 only amount to about 10%: I.e., within our range of z 2 and p · z values, we are more sensitive to higher twist effects and perturbative corrections than we are to the shape of the leading twist DA. This is also expected from Fig. 2 and the discussion of Section II.
The data points in Fig. 7 as well as in Figs. 8-10 below are obtained by performing a weighted average over all possible combinations of the distance z = (z 1 , z 2 , z 3 ) and momentum p = (p 1 , p 2 , p 3 ) that give the same values for the scalar product p · z and the same z 2 . The markers indicate how many different momenta from Table I contribute to the average: dot= 1, cross= 2, triangle= 3, square= 4, pentagon= 5, hexagon= 6. The VV + AA channel yields the best signal by far, since in this case only one invariant structure exists that is consistent with the symmetries, and one can make use of an additional average over the open Lorentz indices, cf. Appendix A. This averaging is not possible in the VA + AV channel since in this case one needs to project onto the specific leading twist Lorentz structure Eq. (A1b). This projection entails a strong dependence of the signal-to-noise ratio on the momentum direction and for some data points none of the analysed momenta yields a good signal. This explains the outliers with extremely large statistical errors in the figures below. Finally, the SP + PS channel, albeit slightly inferior to VV + AA, also gives small statistical errors.

B. Extraction of distribution amplitude parameters
We are now in a position to analyse the whole data set and attempt to extract the pion DA, carrying out a global fit to all correlation functions using the expressions collected in Sec. III B. In Figs. 8-10 we show our data for four distances, along with such fits. The fits A (turquoise), B (orange) and C (red) correspond to different parametrizations of the leading twist pion DA. Ansatz A corresponds to using the power-law parametrization (10) with a free fit parameter α, while B and C use the Gegenbauer expansion (6) truncated at orders n = 2 and n = 4, respectively. All input parameters are taken at the reference scale µ 0 = 2 GeV and are evolved to µ = 2 z using two-loop evolution equations, apart from the higher twist parameter δ π 2 , where the scale dependence is taken into account at one-loop order. This means fits A and B have two free parameters -α, δ π 2 (A) and a π 2 , δ π 2 (B) -while fit C has three parameters: a π 2 , a π 4 , δ π 2 . The results are shown in Table II for different fit ranges in 2 z . The numbers in parentheses are the statistical errors, which turn out to be surprisingly small for a π 2 and also for δ π 2 . The Gegenbauer coefficient a π 4 cannot be constrained from our data and including this contribution (ansatz C compared to B) does not lead to a distinct improvement of the fit quality. The reason is obvious from Fig. 2, as the n = 4 partial wave gives a negligible contribution to the correlation functions in the p · z range accessible in our study. The small statistical errors for a π 2 and δ π 2 are encouraging and allow us to analyse the (dominant) systematic errors. In order to gain some insight, we have performed the complete analysis for multiple fit ranges in the distance between the currents. A dependence on the lower bound in the distance (corresponding to larger scales) can indicate discretization effects, while a dependence on the upper bound shows the necessity to calculate higher- The SP + PS correlator as a function of p · z for four different separations between the currents. The turquoise, orange and red bands correspond to fits using the parametrizations A, B and C explained in the text, cf. Table II. The dashed lines are obtained by subtracting the higher twist contributions from the parametrizations.   (10), while B and C use the expansion of the DAs in terms of Gegenbauer polynomials, Eq. (6), truncated at n = 2 and n = 4. The numbers in parentheses give the statistical error. As discussed in the main text, a rather generous systematic uncertainty of 30%-50% should be assigned to these results and the values for a π 4 from ansatz A and B are meaningless. The fit range corresponding to the curves plotted in Figs order corrections to the coefficient functions and, possibly, even higher twist corrections. Such effects are clearly visible, cf. Table II. As a second method to estimate the systematic uncertainty, one may assume that not yet calculated higher order perturbative effects are of the size of ∼ 50% of the one-loop correction. Both error estimation methods lead to the conclusion that, for the time being, one has to assign a systematic error of at least 30%-50% to the given numbers for a π 2 and δ π 2 , especially since other systematic uncertainties originating from an unphysically large pion mass as well as finite volume and lattice spacing corrections have not been addressed in this study.

C. Discussion
Within the present range of distances and momenta our data appear to be very sensitive to higher twist corrections. These corrections can be quantified within our approach and the corresponding parameter δ π 2 proves to be only weakly correlated with the shape parameters of the pion DA. This can be explained as follows.
First, it is crucial that perturbative and higher-twist corrections for the VV + AA and PS + SP correlators have similar magnitude and opposite sign, cf. Fig. 7. The higher-twist corrections contribute mostly to the difference of these two correlation functions, and much less to their sum. The effect of adding the a π 2 parameter to the leading-twist pion DA is just the opposite, i.e., it affects both VV + AA and PS + SP correlators in a similar way. Second, writing the correlation functions Φ XY π (p·z, z 2 ) as an expansion in conformal partial waves similar to Eq. (8) for the DA, one can include higher-twist terms as contributions O(z 2 ) to the Gegenbauer coefficients, see [12] for details. It turns out that this correction is largest for the leading term a π 0 ↦ a π 0 (z 2 ) = a π 0 + cδ π 2 z 2 + . . . and affects a π 2 and higher coefficients rather weakly. As a consequence, the higher-twist parameter δ π 2 can be extracted from position space correlators at small values of p · z , which explains its small statistical error.
Note, however, that the obtained value is tied to using first order perturbative corrections O(α s ) to the correlators, and will likely decrease if further terms are taken into account. This ambiguity is conceptual. It is related to the fact that matrix elements of twist-four operators have quadratic power divergences already in the continuum theory and at the same time the perturbative series in leading twist in the minimal subtraction scheme suffers from factorial divergences (renormalons). One can show [80] that these two deficiencies are related and are cured in the sum of perturbative (leading-twist) and nonperturbative (higher-twist) effects. The highertwist contribution, strictly speaking, should be viewed as an effective parametrization of the sum of the uncalculated higher orders of perturbation theory and "genuine" higher-twist effects; their separation requires additional regularization and is not necessary in the present context.
Our result for a π 2 has good statistical accuracy and all parametrizations of the DA lead to similar values that are somewhat larger than the result from the direct calculation of the second moment in Ref. [20], a π 2 = 0.1364(154)(145) (at 2 GeV). This should not be viewed as a contradiction as the systematic errors in the present study are not yet under control. They will decrease significantly in future, especially if one could reach values of p · z ≳ 5, which would also allow us to start probing the next Gegenbauer coefficient, a π 4 . The leading twist DAs obtained from ansätze A and B with fit range II are plotted in Fig. 11. Note that the error bands only show the statistical error and that the systematic uncertainty (cf. fit range variation in Table II) is considerably larger. Both DAs shown in Fig. 11 are in perfect agreement with our data since they yield similar values for a π 2 , which is, as discussed above, the only parameter that is relevant for the description of the data within the range of p · z that is currently available. In order to distinguish these DAs from each other, one would need data at larger p · z values that are sensitive to higher Gegenbauer coefficients. Our results favour DAs that, at a scale of 2 GeV, are considerably broader than the asymptotic DA. The turquoise and orange bands correspond to DAs at a reference scale µ0 = 2 GeV obtained from the fits to parametrizations A and B for fit range II, cf. Table II. Both DAs lead to an equally good description of our data because they have a similar second Gegenbauer coefficient a π 2 , which is the only physically relevant information needed from the DA at the available range of p · z. Note that the error only includes the statistical error for the used fit range and that the systematic uncertainty is considerably larger. For comparison, we have also included a result obtained using the quasi-distribution approach (dashed line) taken from [28].

VI. CONCLUSION AND OUTLOOK
In this work we demonstrate that the method proposed in [12] for the determination of collinear parton distributions does not only lead to qualitatively appealing results (see our first article on the topic [1]) but is indeed capable of producing quantitative results with surprisingly small statistical errors. The latter is possible due to the combination of momentum smearing (improving the signal for hadrons with large momentum) with stochastic estimation. A main characteristic of our approach is that we use an equal-time correlation function of two local currents, connected by a light quark propagator, instead of a nonlocal operator, connected by a Wilson line gauge transporter. This has multiple advantages: 1. We circumvent problems originating from the renormalization of nonlocal operators entirely, since the local currents we use can be renormalized using well-tested standard methods.
2. Using a quark propagator easily allows us to evaluate distances that are not aligned with a lattice axis. While this is also possible when using a smeared Wilson line [31], the latter may interfere with the renormalization. On-axis separations are actually the worst case scenario as we find discretization effects to be largest for these directions, cf. Fig. 5. A restriction to the axes also implies a considerable reduction of the data set: in Figs. 8-10 this would correspond to having only one data point per momentum per plot.
3. We can evaluate multiple channels, which gives us an additional handle on the systematic error. It is crucial that higher twist corrections for different correlation functions are related and can have opposite sign. The channels we have analysed lead to consistent results and can be used in a global analysis to obtain values for the leading Gegenbauer coefficient of the leading twist DA and for the higher twist normalization constant. Note that it should be possible to include data from the Wilson line approach as an additional channel in such a global analysis.
4. Using two local currents instead of one nonlocal current has the nice feature that one can in principle apply the usual operator improvement within the Symanzik improvement program to remove O(a) effects.
5. For the matrix element with two local operators finite volume effects have been calculated in [81]. The results therein show that even for an intermediate lattice size with m π L = 4 (in our case m π L ≈ 3.4) one has to expect large volume effects (∼ 10%) once the distance between the two currents approaches half of the lattice extent, i.e., if z ≈ 0.5L. In this respect, it is helpful that our analysis method is restricted to relatively small distances z ≤ 5.5a ≈ 0.2L where perturbative QCD is applicable, meaning that these volume effects are under control.
Another important feature of our analysis method is that we match the perturbative QCD calculation and the lattice data directly in position space. Note that such a position space analysis is not tied to using a light quark propagator, but can also be performed within the Wilson line approach (see, e.g., [34,35]). The obvious advantage over a quasi-distribution type analysis is that one can directly see on the data level (in Figs. [8][9][10] whether the perturbative matching between the off-lightcone correlation function one calculates and the lightcone quantities one is interested in actually works. From a global fit to our data we obtain values for a π 2 and δ π 2 with unexpectedly small statistical errors. An analysis of the fit range dependence showed that we have reached an accuracy where the systematic uncertainties by far dominate. Nevertheless one can say that the value obtained for a π 2 indicates a DA that, at 2 GeV, is considerably broader than the asymptotic one. The value we obtain for the higher twist matrix element (B10), is only slightly larger than sum rule estimates, which lie at approximately δ π 2 = 0.17 GeV 2 at a scale of 2 GeV. To our knowledge this is the first determination of δ π 2 from lattice QCD.
We find that, restricting the analysis to distances where perturbation theory is applicable, even our largest momentum (with p = 2.03 GeV) is still slightly too small for the data to be sensitive to a π 4 . However, it is clear that this situation will improve dramatically if one could reach values of p > 2.5 GeV.
Having reached small statistical errors only to find a large systematic uncertainty may seem a bit unsettling at first. In fact, the opposite is the case, since all main problems we have identified can be solved by systematically improving the analysis and provide us with some guidance towards the next necessary steps. On the lattice side of the calculation, we find discretization effects to be the gravest issue (despite all efforts to tame them described in Section IV B). We plan to address this problem with a twofold strategy, both by drastically reducing the lattice spacing and, in the long run, by implementing O(a) improvement. To reduce the systematic uncertainty from the perturbative side of the calculation our results clearly call for a two-loop analysis and a more systematic study of higher twist effects.

ACKNOWLEDGMENTS
This work has been supported by the Deutsche Forschungsgemeinschaft (SFB/TRR-55), the Polish National Science Center (NCN grant number UMO-2016/21/B/ST2/01492) and the Studienstiftung des deutschen Volkes. We acknowledge PRACE for awarding us access to Marconi-KNL hosted by CINECA at Bologna, Italy. Part of the analysis was carried out on the QPACE 2 [82] Xeon Phi installation of the SFB/TRR-55 in Regensburg. We used a modified version of the Chroma [83] software package along with the LibHadronAnalysis library and the multigrid solver implementation of [84] (see also [85]).

Appendix A: Lorentz-projection operators
In order to project onto T XY defined in Eq. (25), we can use the projection matrices P VA µν = p · z 2((p · z) 2 − p 2 z 2 ) 2 (2(p · z) 2 + p 2 z 2 )(p µ z ν + z µ p ν ) − 3p · z(z 2 p µ p ν + p 2 z µ z ν ) − p · z((p · z) 2 − p 2 z 2 )g µν , such that T XY = P XY µν T µν XY . For the vector-axialvector channel this projection is the only possibility to obtain T VA . In the case of the vector-vector channel (and the axialvector-axialvector channel), however, one can obtain T VV (or T AA ) from any channel with two fixed indices µ and ν as long as ε µνρσ p ρ z σ ≠ 0: In our final analysis we use a weighted average of the individual channels, where the weight is defined as the inverse standard deviation squared of the respective channel. This yields a much better signal than the projection with (A1a), which basically averages over all vectorvector/axialvector-axialvector channels.

Appendix B: Higher twist corrections
In this appendix we provide some details on the calculation of the higher twist corrections. First of all, for nonvanishing quark masses also the chiral odd twist-three pion DAs have to be taken into account: They enter our calculation multiplied by the quark mass and become part of the pion mass correction. Since these contributions are small we have used the simplest asymptotic expressions, and omitted corrections due to the three-particle quarkantiquark gluon DA [57]. Complete expressions for the twist-three matrix elements can be found in [57,58].