Pion distribution amplitude from Euclidean correlation functions

Following the proposal in [1], we study the feasibility to calculate the pion distribution amplitude (DA) from suitably chosen Euclidean correlation functions at large momentum. In our lattice study we employ the novel momentum smearing technique [2,3]. This approach is complementary to the calculations of the lowest moments of the DA using the Wilson operator product expansion and avoids mixing with lower dimensional local operators on the lattice. The theoretical status of this method is similar to that of quasi-distributions [4], which has recently been applied to the same problem in [5]. The similarities and differences between these two techniques are highlighted.


Introduction
In recent years there has been increasing interest in the possibility to determine parton distribution functions from Euclidean correlation functions, bypassing Wilson's operator product expansion. The general scheme of such calculations is to consider a product of suitable local currents at a spacelike separation z, sandwiched between hadronic states, and match the lattice calculation of this quantity to the perturbative expansion in terms of collinear parton distributions H |q(n)Γ q(−n)|H , The existing concrete proposals differ mainly in the choice of the Q-field. This can be chosen as an auxiliary scalar in a e-mail: philipp.wein@physik.uni-r.de the fundamental representation of the color group [6,7], or as an (auxiliary) heavy [8] or light [1] quark. Another suggestion [4] is to replace the Q-field propagator by a Wilson line connectingq(z/2) and q(−z/2). This last proposal received the most attention, despite added complications due to the renormalization of the Wilson line [9,10], see Refs. [11][12][13][14][15][16][17][18] for recent discussions, the reason being that it allows for a more direct momentum space interpretation in the framework of the large-momentum effective theory [14,19] (LaMET). The corresponding correlation functions, transformed into the longitudinal momentum fraction representation, have become known as quasi-parton-distributions [4]. While quasi parton-distributions are certainly interesting objects, it was already pointed out in Refs. [1,20] that position space correlation functions (or "lattice cross sections", in the terminology of [21][22][23]) contain the complete information on parton distributions. In all cases the functions calculated on the lattice (for early work, see also [24]) are related to parton distributions by means of QCD factorization in the continuum, which can be done both in position and momentum space. A position space analysis naturally leads to the concept of Ioffe-time distributions [16,20,25,26]. We emphasize that all the above suggestions are equivalent, and their relative virtue will be determined by the possibility to control lattice artifacts and other systematic uncertainties.
In this work we study the simplest function of this kind, the pion distribution amplitude (DA), using the technique suggested in Ref. [1], i.e., we use a light quark (Q = q) in Eq. (1) and perform the analysis directly in position space. The same DA has recently been studied using the quasidistribution approach in Ref. [5]. We consider a correlation function of renormalized scalar and pseudoscalar operators at equal times whereq(z) creates a light quark field of hypothetical flavor q = u, d and square brackets [O] denote operator renormalization in the MS scheme. In what follows, we fix the renormalization scale to the "kinematic" scale in the correlator, μ R = 2/ √ −z 2 (cf. [1,27]). The correlation function (3) can be calculated on the lattice as a function of two variables, p · z = −p·z and z 2 = −z 2 . Here and below we use boldface letters for spatial 3-vectors.
We restrict ourselves to sufficiently small distances, |z|/2 < 1 GeV −1 , such that the same correlation function can be calculated in continuum perturbation theory in terms of the pion DA using standard QCD factorization techniques. The result reads where Φ SP π = Φ π + O(α s ) + higher twist (the various corrections will be discussed later), with F π ≈ 93 MeV is the pion decay constant, the variable u corresponds to the quark momentum fraction and φ π (u) is the (leading-twist) pion DA. The integral of the pion DA is normalized to unity, 1 0 du φ π (u) = 1, and its shape has been hotly debated for more than 30 years. This discussion has been reinvigorated by the strong scaling violation in the πγ * γ form factor observed by the BABAR [28] and, to a lesser extent, the BELLE [29] collaboration, which is difficult to explain unless the pion DA exhibits strong enhancements near the end points, see, e.g., Refs. [30][31][32] for a review and further references. For illustrative purposes we consider three models:   Fig. 2 The position space pion DA Φ π ( p · z) [cf. Eq. (5)] for the three models in Eq. (6) φ (2) at the reference scale μ 0 = 1 GeV. These models and the corresponding Fourier-transformed position space DAs Φ π ( p · z), defined in Eq. (5), are plotted in Figs. 1 and 2, respectively. Measuring the correlation function (3) on the lattice for a range of values of p · z = −p · z and z 2 = −z 2 gives access to the pion DA in position space (5) that contains the full information on the longitudinal momentum fraction distribution.
The main difference of our technique [1] to the approach of Ref. [5] is that the smallness of higher twist and perturbative corrections (for arbitrary pion momentum) is guaranteed by keeping the distance |z| between the currents sufficiently small. A large pion momentum is needed not in order to suppress the corrections, but because it provides the necessary lever arm in the dimensionless variable p · z that is mandatory to distinguish between pion DAs of different shape, see Fig. 2.
In contrast, in the LaMET-based approach of [5] formally a Fourier transform over all values of z is taken, and smallness of perturbative and higher twist corrections is achieved indirectly by considering the asymptotic expansion of the amplitude at large values of the pion momentum for a fixed quark momentum fraction u, which is the Fourier conjugate of p · z. Thus |p| → ∞ implies that the integration region in the Fourier integral shrinks to |z| ∼ 1/|p| → 0.
Another difference is that in Ref. [5] a Wilson line is used to connect the quark and the antiquark, whereas in this study we use a light-quark propagator [1]. To tree-level accuracy the difference in the corresponding coordinate space expressions is simply a different coefficient function in Eq. (4). While taking the Wilson line not along a lattice axis introduces additional difficulties, e.g., concerning renormalization, the separation z can be chosen arbitrarily without prob-lems when a light-quark propagator is used. We consider this possibility as an advantage of our calculation because we have found discretization errors to be largest if z lies along a lattice axis. Also the renormalization of the lattice correlator is greatly simplified when one works with a lightquark propagator (for recent progress regarding the Wilson line approach see [16][17][18]).
Note that we suggest to match the lattice matrix element with the pQCD factorization expression directly in coordinate space. This has the advantage that the lattice data can be directly confronted with the theory since perturbative predictions based on model parametrizations of the DAs can easily be transformed to position space.
The whole program naturally splits into two parts-the lattice calculation where all usual extrapolations/limits have to be taken and the pQCD factorization in terms of the pion DA in the continuum. Our presentation is structured accordingly.

QCD factorization
The complete QCD expression for the correlation function (3) can be written as where H (u, z 2 , μ) = 1 + O(α s ) is a short distance coefficient function that can be evaluated perturbatively, μ is the factorization scale, and T HT stands for power-suppressed (in z 2 ) contributions that can be calculated in terms of the pion DAs of higher twist [33,34]. The factorization scale dependence of the pion DA is considerably simplified by using the expansion where the C 3/2 n (x) are Gegenbauer polynomials. The n = 0 coefficient is fixed to unity, a π 0 = 1, by the normalization condition and the remaining ones, n = 2, 4, . . ., encode all relevant nonperturbative information on the DA. They have to be defined at a certain reference scale μ 0 (a common choice is μ 0 = 1 GeV) and evolved to the scale of the process. The corresponding mixing matrices are known in analytic form to two-loop accuracy [35,36] and numerically for the first few moments to three-loop accuracy [37].
Using this expansion, the leading-twist (LT), i.e., twist two, contribution to the correlation function can be written as Setting both the renormalization and factorization scales to μ ≡ 2/ √ −z 2 we obtain, to O(α s ) accuracy, where The plus prescription is defined as usual, The sum in (9) converges very rapidly since so that for finite ρ ∼ 1 2 p · z only the first few Gegenbauer moments give a sizeable contribution, cf. [1].
The leading higher twist contribution O(z 2 ) can be estimated using models for the twist 4 pion DAs discussed in Refs. [33,34]. For the case at hand these corrections are in general complex. We obtain for the real part where δ 2 π 0.2 GeV 2 at the scale μ = 1 GeV [33,34]. The last two terms take into account the pion mass corrections and are rather small.
We find that the higher twist correction for the scalarpseudoscalar correlation function has the same sign as the leading-twist term, in contrast to the vector-vector correlation function considered in Ref. [1], in which case the higher twist correction has the opposite sign. Numerically, this correction turns out to be about 20% of the leading twist contribution at |z|/2 ∼ 0.2 fm 1 GeV −1 .

Generalities
We wish to avoid the calculation of disconnected quark line diagrams, which are challenging in lattice simulations. This becomes possible by implementing an appropriate flavor structure of our currents. One may consider having a π 0 in the final state and q = d in Eq. (3). However, this matrix element vanishes identically due to isospin symmetry. Instead we pretend that the auxiliary quark field q of Eq. (3) is a different, third flavor but for simplicity we keep it at the same mass m q = m u = m d . This corresponds to our continuum QCD calculation.
In the actual lattice calculation we determine the threepoint function using the sequential source method. Therefore, the currents are situated at z and at (the chosen origin) 0 and are afterwards "shifted" to the symmetric locations in Eq. (3) by multiplication with the appropriate phase.

Correlation functions
The remaining nontrivial part of the lattice simulation is the calculation of the connected triangle diagram depicted in Fig. 3. Introducing a phase matrix ϕ t that is diagonal in position space (with diagonal entries (ϕ t ) yy = e −ip·y ) and is nonzero only on time slice t, we can rewrite the three-point function as follows: where z = (z, 0), G stands for the quark propagator, and is a sequential source. The momentum dependent smearing Φ (p) is performed as described in Ref. [2]. We want to stress that in this situation the new momentum smearing technique is even more cost efficient than described in Ref. [3], since one needs a second inversion of the Dirac operator for each additional momentum anyway. The matrix element (3) can be obtained from C 3pt by canceling the normalization factor describing the overlap of the smeared current with the pion state. The latter can be obtained, e.g., from the two-point function C 2pt of a smeared current at the sink (at time t, as in the three-point function) and a local axialvector current at the source. Neglecting excited state contributions, one finds where Z X is the renormalization factor of the local current X with respect to the MS scheme [38], cf. Sect. 4.

Taming discretization effects
In the continuum, the chiral even part of the propagator connecting the two local currents (proportional to / z) gives the most important contribution, while the chiral odd part (proportional to the unit matrix) 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 differ-  Correction of discretization effects for the example of a fixed p · z = 0.39 and |p| = 1.08 GeV. Points with a correction larger than 10% (gray triangles) will be ignored in the analysis ent. We find that the contribution from the chiral odd part, which removes the doublers and breaks chiral symmetry, can be of the same order of magnitude as the leading contribution, cf. Fig. 4. The "jumping" of the points nicely demonstrates the strong dependence of the lattice artifacts on the chosen direction. In particular the points along the axes [e.g., (1, 0, 0)] 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. However, the appearance of large discretization effects is probably a general feature of all coordinate space methods. The appearance of large contributions from the chiral odd part of the propagator would lead to huge lattice artifacts in the correlator. However, the perturbative calculation shows that in the situation where the two currents are located symmetrically with respect to the chosen origin, contributions from the chiral even and the chiral odd parts of the propagator are nicely separated in some channels. For instance, for the scalar-pseudoscalar channel the contribution from the chiral even part (which is the one we are interested in) is real, while the chiral odd part appears only in the imaginary part. Hence, we can choose to analyze only the part of the signal that does not contain the problematic contributions. Note, however, that the continuum expectation that either the real or the imaginary part (depending on which one corresponds to the chiral odd part) of the signal should be strongly suppressed is not valid for the lattice data. Therefore, the correct identification of the relevant part of the signal is crucial.
We can now concentrate on the correction of the discretization effects in the chiral even part of the propagator (these correspond to the blue points in Fig. 4). First and foremost we have decided to simply discard data points where the free field discretization effect is already larger than 10%, which mainly excludes very small distances (|z| 2a) and directions along the lattice axes. For the remaining data points we use a correction factor c corr (z) determined such that the corrected propagator satisfies the condition where the trace runs over Dirac and color indices. To zeroth order accuracy in α s (where G latt = G free latt is the free propagator) this leads to c corr (z) = tr D / zG free latt (z) Fig. 6 The plots show our results for the scalar-pseudoscalar channel for different fixed values of p · z and different pion momenta. The errorbars include the statistical error only. The solid/dashed lines correspond to predictions taking into account/neglecting higher twist contributions for the different DA models (6). The color-coding is the same as in Fig. 1 which corresponds to multiplying the blue data points in Fig. 4 with a factor such that one obtains the continuum result in the non-interacting case. Looking at Fig. 5 it is clear to the naked eye that this correction leads to a much smoother and less direction-dependent behavior of the data points.

Results
The gauge field ensemble used in this study has been generated (by QCDSF / RQCD) with two mass-degenerate flavors of nonperturbatively improved Wilson fermions and the Wilson gluon action (ensemble IV in Ref. [39]). The dimensions of the lattice are 32 3 × 64 and the hopping parameter is κ = 0.13632. The coupling parameter β = 5.29 translates to the lattice spacing a ≈ 0.071 fm = (2.76 GeV) −1 and the pion mass has been determined in Ref. [40] to the value m π = 0.10675(59)/a ≈ 295 MeV. In order to get a reasonable overlap with the hadron state at large momentum, we have employed the momentum smearing technique (cf. Ref. [2]) with APE smeared links [41]. The operator renormalization is performed as described in Ref. [38]. The local operators are renormalized nonperturbatively in a RI -MOM scheme along with a subtraction of lattice artifacts in one-loop perturbation theory. The final conversion to the MS scheme employs 3 loop continuum perturbation theory. To be consistent, we use the N f = 2 specific running of α s in all perturbative calculations. To this end, we combine the results of Refs. [42] and [43] to obtain a value of α s at 1000/a ≈ 2.76 TeV. From there we evolve it downwards using 5 loop running.
In Figs. 6 and 7 we confront the data points with predictions from continuum perturbation theory corresponding to the pion DAs shown in Fig. 1. For all cases we show a version ignoring higher twist (i.e., twist 4) effects (dashed lines) and one including higher twist corrections (solid lines), where δ π 2 (1 GeV) = 0.2 GeV 2 is set to the QCD sum rule estimate obtained in Ref. [44] (cf. also Ref. [45]). The errorbars only include the statistical error. It is clear that the systematic uncertainty is sizeable: in addition to discretization effects and higher order perturbative corrections there may be excited state contaminations and, possibly, finite volume effects. Therefore, one should refrain from drawing any premature phenomenological conclusions. Nevertheless, the qualitative agreement found in Ref. [46] between our data points at μ = 1.08 GeV and the results obtained using the quasi-DA method is encouraging. In Fig. 6, one immediately notices the higher twist effect, in particular for large distances, while the curves corresponding to the various DAs are hard to distinguish for small values of p · z. In Fig. 7 deviations from the asymptotic form are nicely visible at p · z 4. The latter region, however, can only be reached with larger hadron momenta, since we are limited to perturbatively accessible values of |z| = 2/μ 2 GeV −1 ≈ 5.5a. In Fig. 7 it becomes clear that our data points with |p| = 1.88 GeV can already reach out into this region, but that one still needs higher statistics to be able to differentiate between different DA models.

Summary
In this work we have demonstrated that the coordinate space method for the determination of the pion DA proposed in Ref. [1] is promising, in particular as far as the statistical error is concerned. To this end, we have analyzed lattice data at m π = 295 MeV at a lattice spacing of a = 0.071 fm using dynamic Wilson fermions. We have shown that the large hadron momenta, which are a prerequisite of this method (and also for other related methods), lie just within the scope of the novel momentum smearing technique [2].
We have found particularly large discretization effects for directions along the coordinate axes, which are probably not specific to our calculation but will most likely occur also in other coordinate space calculations. Furthermore, we find that the chiral odd part of the quark propagator leads to large lattice artifacts. This contribution stems from the Wilson term in the propagator and is therefore a peculiarity of using Wilson fermions. We have overcome this problem by analyzing the real part of the scalar-pseudoscalar channel, where only the chiral even part contributes. For the remaining discretization effects stemming from the propagator, we have adopted the correction method described in Sect. 3.3, which has reduced the anisotropy of the data considerably. However, observing large discretization effects on this single intermediate lattice spacing shows that taking the continuum limit will be of vital importance, if one aims at achieving quantitative results in the future.
Unlike the Wilson line approach [5,16], the direction of the separation between the currents z can be chosen arbitrarily using our method. This enables us to realize a large number of different |z| and p · z values and also to study and minimize discretization effects. Intricacies related to the Wilson line renormalization [11][12][13][14][15][16][17][18] are avoided entirely, and the possibility to vary the Dirac structures in the currents offers an additional handle on the higher-order perturbative corrections and higher twist effects.
In the near future we plan to investigate a new algorithm that may reduce the statistical uncertainties. We will also move to a smaller lattice spacing to enable the use of larger momenta |p| π/a (and therefore larger p · z values at a given scale μ = 2/|z| 1 GeV) along with distances |z| a, which will reduce lattice artifacts. Fig. 7 Results from the scalar-pseudoscalar channel for fixed values of the distance z 2 corresponding to perturbative scales around μ ≈ 1GeV. The errorbars only include the statistical error. Note that this comprises only a small subset of the available data. The left-hand-side plots display the real part of the normalized matrix element T (p · z, z 2 )/F π , while those on the right-hand-side show the respective real part of Φ SP π (p · z, z 2 ) defined in Eq. (4). At tree level and up to higher twist corrections, the latter correspond to the position space DA (5). Solid/dashed lines correspond to predictions that include/neglect higher twist contributions. The color-coding is the same as in Fig. 1 deutschen Volkes. The analysis was carried out on the QPACE 2 [47] Xeon Phi installation of the SFB/TRR-55 in Regensburg. We used a modified version of the Chroma [48] software package along with the LibHadronAnalysis library and the multigrid solver implementation of Ref. [49] (see also Ref. [50]). We thank Daniel Richtmann for code development, discussions and software support.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .