Second moment of the pion distribution amplitude with the momentum smearing technique

Using the second moment of the pion distribution amplitude as an example, we investigate whether lattice calculations of matrix elements of local operators involving covariant derivatives may benefit from the recently proposed momentum smearing technique for hadronic interpolators. Comparing the momentum smearing technique to the traditional Wuppertal smearing we find - at equal computational cost - a considerable reduction of the statistical errors. The present investigation was carried out using $N_f=2+1$ dynamical non-perturbatively order $a$ improved Wilson fermions on lattices of different volumes and pion masses down to 220 MeV.


Introduction
Many quantities of interest in high-energy physics involve hadrons carrying large momenta. The prime example is provided by form factors, but also parton distribution functions (PDFs) and their generalizations, in particular transverse momentum dependent parton distribution functions (TMDs) receive their physical interpretation in the large-momentum limit.
Very high accuracy is expected for future experimental data, e.g., on hard exclusive and semi-inclusive reactions at the JLAB 12 GeV upgrade [1] and at the Electron Ion Collider (EIC) [2], as well as on B-meson decay and pion transition form factors at Belle II at KEK [3]. This accuracy has to be matched by an increased theoretical precision. Such processes are usually studied using factorization techniques, where the nonperturbative input is reduced to operator matrix elements which, ideally, should be computed using lattice QCD. Also in the cases where no momentum transfer takes place between the initial and the final state one usually needs to realize hadron sources with nonvanishing momenta in order to have the possibility to employ operators with sufficiently simple renormalization patterns. It has also been argued [4,5] that hadron sources with large momenta offer novel opportunities, enabling a more direct calculation of parton distributions and hadronic light-cone wave functions by performing a collinear factorization of suitably chosen Euclidean correlation functions (e.g., "quasi-PDFs" [5]), thereby circumventing the traditional Wilsonian local operator product expansion. †

RQCD Collaboration
Email address: fabian.hutzler@ur.de (F. Hutzler) Although the problem is known for quite some time, up to very recently [6] no satisfactory techniques for hadrons carrying high momenta on the lattice existed to suppress excited state contributions while maintaining acceptable signal-to-noise ratios. The generic method of reducing excited state overlaps consists of employing carefully tuned extended interpolators. However, for larger momenta the usual smearing techniques become increasingly less effective. The basic idea of Ref. [6] was to modify the usual quark smearing functions by additional phase factors such that the centre of the distribution in momentum space is shifted towards the desired value. By implication, such smearing functions correspond to oscillating wave packets in position space.
It was shown that this technique, which we will refer to as momentum smearing, leads to considerably improved signal-to-noise ratios for the pion and the nucleon twopoint functions [6] as well as for lattice observables that are related to quasi-PDFs [7]. In this letter we address another class of applications, namely computing hadronic matrix elements that contain local operators with covariant derivatives, e.g., moments of parton distributions and distribution amplitudes (DAs) [8][9][10][11][12][13]. In the case that we specifically study, i.e., moments of DAs, the matrix elements of interest are proportional to powers of the hadron momentum and are known, empirically, to be very noisy when using traditional methods. We will demonstrate that momentum smearing results in a major improvement of the quality of the signal for the second moment of the pion DA. In fact, it turns out that this technique is so effective that, at small pion masses and large lattice volumes, statistical fluctuations can be further reduced by deliberately selecting a momentum that is larger than the smallest possible choice.
The scope of the present study is mainly methodological. In addition, we present the first lattice calculation of the 2nd moment of the pion DA using N f = 2+1 dynamical clover Wilson fermions. The results are compatible with the latest N f = 2 study [12], while the second moment is somewhat smaller than what has been reported in a simulation employing N f = 2 + 1 domain wall fermions, which has been carried out at a coarser lattice spacing and at larger quark masses [10].

Continuum definitions
Pseudoscalar mesons like the pion have only one independent leading twist (twist two) DA, φ, which is defined via a meson-to-vacuum matrix element of renormalized non-local quark-antiquark light-ray operators, where z 1,2 are real numbers, n µ is an auxiliary light-like vector with n 2 = 0, and π + (p)⟩ represents the ground state pseudoscalar π + meson with on-shell momentum p 2 = m 2 π . The straight path-ordered Wilson line connecting the quark fields, [z 2 n, z 1 n], is inserted to ensure gauge invariance. The scale dependence of φ is indicated by the argument µ 2 .
Neglecting both isospin breaking and electromagnetic effects, the DAs of the charged pseudoscalar π ± and the neutral π 0 are trivially related such that it is sufficient to consider only one of them. The decay constant f π appearing in Eq. (1) can be obtained as the matrix element of a local operator, and has a value of f π ≈ 130 MeV [14]. The physical interpretation of Eq. (1) is that the fraction x of the pion momentum is carried by the u quark, while thed antiquark carries the remaining fraction 1 − x. Hence the difference of the momentum fractions, contains all nontrivial information and its moments are defined as Since the Gegenbauer polynomials C 3 2 n (2x − 1), which correspond to irreducible representations of the collinear conformal group SL(2, R), form a complete set of functions, the DAs can be expanded as where the Gegenbauer moments a n renormalize multiplicatively in leading logarithmic order. Higher-order contributions in the Gegenbauer expansion are suppressed at large scales, since the anomalous dimensions of a n increase with n. Hence, in the asymptotic limit µ → ∞ only the leading term survives, which gives:

Lattice definitions
From now on we will work in Euclidean spacetime and follow the conventions of Ref. [12]. The renormalized lightray operator on the left-hand side of Eq. (1) generates renormalized local operators. This means that Mellin moments of the DAs, see Eq. (4), can be expressed in terms of matrix elements of local operators and can be evaluated using lattice QCD. In order to calculate the second moment of the pion DA (n = 2), we define the bare operators where D µ is the covariant derivative, which will be replaced by a symmetric discretized version on the lattice. In order to obtain a leading twist projection we symmetrize over all Lorentz indices and subtract all traces. This is indicated by ρµν can also be written as The operator O + ρµν is, in the continuum, given by the second derivative of the axial-vector current: This is not the case on the lattice due to discretization effects of the derivatives which can be numerically sizable. The mixing with operators of lower dimension can be prevented by selecting lattice operators that belong to a suitable irreducible representation of the hypercubic group H(4) [9,10]. For our case, this corresponds to choosing all indices different for the operators O ± . Identifying one index with the temporal direction, this leaves us with the operators In order to extract the desired moments we use twopoint correlation functions of the operators O ± 4jk and A ρ with an interpolating field, where J = P or J = A 4 . For sufficiently large t, the ground state dominates and the correlation functions give where the sign factors τ O , τ J = ±1 depend on the transformation properties of the correlation functions under time reversal. Following Ref. [12], the required matrix elements for the second moments can be extracted from the ratios where ⟨ξ 2 ⟩ bare = R − and a bare 2 = 7 12 (5R − − R + ). In our calculations we use the interpolator J = P, as this gives a better overlap with the ground state than A 4 .
The renormalized moments in the MS scheme read where ζ ij are ratios of renormalization constants which are defined in Ref. [12].

Momentum smearing
On a lattice of N 3 s N t sites, separated by the lattice constant a, the linear spatial extent is given as L = N s a and spatial momentum components are quantized in terms of integer multiples of 2π L. The calculation of the second moment of the DA requires a spatial momentum p = (2π L)n p , with at least two non-vanishing components, i.e., n 2 p ≥ 2. This, in addition to employing two derivatives, considerably deteriorates the signal-to-noise ratio. This problem is ameliorated by using momentum smearing [6]. Here we briefly summarize this method.
It is well known that spatially smearing the quark creation and destruction operators used within the construction of hadronic interpolating fields increases the overlap of the generated superposition of hadronic states with the ground state within a given channel. This is not surprising, as ground state hadrons have smooth, spatially extended wave functions. The smearing operator F should be self-adjoint, gauge covariant and a singlet with respect to all global transformations that act on a timeslice. In the non-interacting case its action on a quark field q can be expressed as a convolution with a scalar kernel function f : In momentum space this convolution becomes a product. If our smearing kernel is a real Gaussian, then in momentum space it will remain a Gaussian centred around k = 0. If the hadron carries a non-vanishing momentum p, it is natural to assume that the quark will carry a momentum fraction k = ζp. We remark that there is no obvious relation between ζ and the longitudinal momentum fraction x of the light-cone wave function. A Gaussian wave function with width σ that is centred about the momentum k acquires a phase: where f (0) = f . Our periodic lattice appears to imply a quantization of the possible values of k. However, Eq. (21) can also be cast into an iterative process, lifting this limitation: It is well known that in infinite volume the above convolution F (k) q can be obtained as the result of evolving the heat equation with a drift term, starting from a spatial delta source at τ = 0, to the fictitious time τ = σ 2 (2α).
One can approximate Eq. (22) by a discrete process, defining F (k) = Φ n (k) as the nth application of an elementary iteration, In practice this smearing is implemented by multiplying the spatial connectors within the timeslice in question by the appropriate phases, U x,j ↦ e −iakj U x,j . For k = 0 Eq. (23) corresponds to the well-known Wuppertal smearing [15,16]. The time coordinate is suppressed as the smearing is local in time.
The gauge connectors within Eq. (23), U x,j and U x,−j ≡ U † x−,j , where denotes a vector of length a and direction j, are spatially APE smeared [17]: where i ∈ {1, 2, 3} and j ∈ {±1, ±2, ±3}. The sum is over the four spatial "staples" surrounding U x,i , and P SU (3) is a gauge covariant projector onto the gauge group SU(3), defined by maximizing Re Tr{A † P SU(3) (A)}. If the APE smeared links are close to unit fields then the width parameter of the resulting Gaussian is given as where large values of ε will allow for smaller iteration counts n, but the resulting function will be less smooth. 1 The root mean squared width of the resulting Gaussian will correspond to √ 3σ as we have three spatial dimensions. This will shrink by a factor 1 √ 2 if we consider the squared wave function and since we will smear both quark and antiquark, the pion interpolator will be wider by a factor √ 2 than the individual quark fields. In the meson case the quark creation operator at the source needs to be smeared with F (k) and the quark destruction operator with F (−k) , while for baryons all three quarks should be smeared with F (k) , see Ref. [6] for details. 2

Results
We illustrate the reduction of statistical errors of the two-point functions that enter the calculation of the second moment of the pion DA, using the momentum smearing technique. For this purpose we consider four Coordinated Lattice Simulations (CLS) ensembles, listed in Table 1. These range from N c ≈ 1500 to N c ≈ 2800 configurations, separated by four hybrid Monte Carlo molecular dynamics units. The statistical errors were evaluated using the Bootstrap procedure, with N samples = 500, combined with the binning method. For the latter we have observed that a binsize of N bin = 10 saturates the statistical error.
After studying the improvement achieved through momentum smearing, we attempt a chiral extrapolation of our results. Since we are working at a fixed lattice spacing a ≈ 0.0857 fm, we cannot as yet perform a continuum limit extrapolation.

Optimizing the smearing and the momentum
In order to obtain the second moment of the pion DA we compute ratios of two-point functions that are smeared at the source and local at the sink (smeared-point), where 2 The sign of the complex phase in Eqs. (21), (22) and (23) is opposite to that of Ref. [6]. Here we assume that the phase of the momentum projection at the sink reads e −ipx and k = ζp with ζ ≥ 0. The phase used in Ref. [6] corresponds to the non-standard e +ipx convention that is used within the Chroma software suite [18]. the physical momenta p and smearing vectors k are parallel: One may naively expect that a value ζ ≲ 1 2 was optimal, evenly distributing the meson momentum between quark and antiquark, however, Ref. [6] indicated that a value ζ ≈ 0.8 was preferable. We confirm this choice: Decreasing ζ from 0.8 to 0.6 we found no improvement of the ground state overlap but slightly increased statistical errors, see Fig. 1 for the example of the ratio R − , defined in Eq. (17). In the following we therefore set ζ = 0.8. Using the momentum smearing for mesons, one needs two inversions per momentum vector n p . In contrast to baryonic two-point functions, where all quarks propagate in the forward direction and therefore are smeared using f (k) , the antiquark in mesonic two-point functions needs to be smeared with f (−k) .
It is instructive to determine which momentum vector n p produces the best signal for a given ensemble. Making the crude approximation of a time-independent noise function, the signal-to-noise ratio of the numerator that dominates the error of the combination Eq. (17) can be estimated as Maximizing this expression with respect to p 2 gives the positive solution Clearly, the optimal choice of momentum for a given correlation function depends on t and lower momenta will always be preferred at large values of t. This means the outcome will depend on the fit window in t and this in turn will depend on the available statistics. To aid in finding the most appropriate momentum, we plot Eq. (28) in Fig. 2 for the typical fit range for our different ensembles, t = 4a-12a. Based on this model, we can read off that for the L = 32a lattices squared momenta in the vicinity of n 2 p = 2 should give reasonable results, whereas for the larger L = 48a lattice values of n 2 p closer to 5 should be investigated. As an example, in Fig. 3 we show the results of the bare observables R ± calculated for different momenta p on the L = 48a C101 ensemble. For small values of t a, larger p 2 exhibit smaller statistical errors, whereas for large values of t a, the error increases with p 2 .
In our further analysis we choose n p = (1, 1, 0), n p = (1, 0, 1) and n p = (0, 1, 1) for the ensembles with a spatial extent of L = 32a and n p = (2, 1, 0), n p = (2, 0, 1) and n p = (0, 2, 1) for the C101 ensemble. For the L = 32a lattices we employ a single source position, while for C101 we realize on average 2 source positions on each configuration. Figure 4 shows the plateau of R + (left) and R − (right) for the H105 lattice with N c = 2830 for both smearing methods. Clearly the momentum smearing generates a much cleaner and longer plateau with very small statistical errors. In contrast to the standard Wuppertal smearing, where errors increase rapidly for high values of t, the signal-to-noise problem is less severe for the momentum smearing. Moreover, in some cases, such as R + (shown in Fig. 4), we notice reduced contaminations from excited states.

Momentum smearing versus Wuppertal smearing
In Fig. 4 we compared the results for the same number of configurations. However, the novel momentum smearing method is computationally more expensive. In general one can average over momenta that are equivalent in terms of the cubic symmetry group O h . Taking into account that the results for p and −p are trivially related, this gives, depending on the momentum, up to 12 possible lattice directions. For the Wuppertal smearing, additional momenta are computationally almost for free, as they only require additional Fourier sums. In contrast, for the momentum smearing each momentum direction requires new, differently smeared sources. For the pion two inversions, with momenta k and −k, are necessary as discussed in Sec. 2.3. For n 2 p = 2 this means that momentum smearing is by a factor of almost 6 more expensive than Wuppertal smearing. Therefore, in Table 2 we provide an equal cost comparison of the ratios of errors obtained using both methods for the H105 lattice. Even at equal cost, we still see a reduction of the squared error by up to a factor 3, in particular for the physically more relevant R − ratio. Note that for mesons containing non-degenerate quarks, the traditional method becomes more expensive as this will also require two inversions, while for baryon interpolators no momentum smearing with −k is required. This means that in terms of a real cost comparison the pion is the least favourable case for momentum smearing.
For a fixed number of measurements the gain of momentum smearing is even larger than at a fixed computational cost. However, the reduction of errors that can be achieved by increasing the number of measurements on each configuration is limited, as additional measurements will become increasingly correlated.

Chiral extrapolation
We use Chiral Perturbation Theory (ChPT) to extrapolate our results to physical quark masses. The CLS lattice ensembles used in this work are chosen such that they lie on the Tr M = const. line [19,20], which means that to next-to-leading order SU(3) ChPT the average quadratic meson mass, is kept fixed at its physical value, up to lattice spacing effects. Thus the extrapolation m π → m phys π also corresponds to m K → m phys K . Up to one-loop order, ⟨ξ 2 ⟩ MS and a MS 2 do not contain chiral logarithms [21], and we assume a linear behaviour in m 2 π , where ⟨ξ 2 ⟩ (n) are LECs of the fit. The chiral extrapolation is depicted in Fig. 5. At the physical point we find  We remark that these numbers were obtained at the fixed lattice spacing a ≈ 0.0857 fm and no continuum limit has been performed yet.

Summary
In this work we have illustrated the effectiveness and advantages of the novel momentum smearing method compared to the standard Wuppertal smearing. For the special case of the pion the momentum smearing requires more inversions, relative to Wuppertal smearing, than for baryons or for mesons consisting of mass non-degenerate quarks. Nevertheless, we have still obtained smaller errors at a similar computational effort. Clearly, using the momentum smearing technique on a fixed number of available configurations, much smaller statistical errors can be achieved. Since for each momentum a new inversion is required in any case, one may suspect that combining the momentum smearing method with the stochastic one-endtrick [22] even bigger gains can be achieved. We have not investigated this possibility as yet.
In future studies the momentum smearing will be applied to mesons and baryons on additional CLS lattices, including ensembles at (nearly) physical quark masses and various lattice spacings down to a ≈ 0.04 fm. This will expand our previous work on meson and baryon DAs [12,13] and enable us to perform systematic continuum limit extrapolations for mesons and octet baryons.