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

TL;DR

This paper demonstrates that the momentum smearing technique significantly reduces statistical errors in lattice QCD calculations of the pion distribution amplitude's second moment, improving efficiency over traditional methods.

Contribution

It introduces and tests the momentum smearing technique for lattice calculations, showing its advantages over traditional smearing methods in reducing errors.

Findings

Momentum smearing reduces statistical errors compared to Wuppertal smearing.

The technique is effective across different lattice volumes and pion masses.

It enables more precise lattice calculations of hadronic matrix elements.

Abstract

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.

Figures

Tables

Table 1. Table 1: List of the ensembles used in this work. β = 3.4 𝛽 3.4 \beta=3.4 corresponds to the lattice spacing a ≈ 0.0857   fm 𝑎 times 0.0857 femtometer a\approx$0.0857\text{\,}\mathrm{fm}$ and N c subscript 𝑁 𝑐 N_{c} denotes the number of analysed configurations. A detailed description of these ensembles can be found in Refs.  [ 19 , 20 ] .

$𝒌=\zeta 𝒑.$

(26)

Table 2. Table 2: Equal cost comparison of the errors obtained with both methods on the H105 lattice for R ± superscript 𝑅 plus-or-minus R^{\pm} . We show the squared ratio of the statistical errors of Wuppertal smearing over momentum smearing divided by the number of inversions needed for each method.

$2m_K^2+m_{\pi }^2,$

(29)

Equations

$$\begin{array}{rl}& ⟨0|\mathrm{bar}d(z_{2}n)/n[z_{2}n,z_{1}n]{\gamma }_{5}u(z_{1}n)|{\pi }^{+}(p)⟩=\\ & =i{f}_{\pi }(p\text{mbox}\text{textperiodcentered}n){\int }_0^{1}dx{e}^{-i(z_{1}x+z_2(1-x))p\text{mbox}\text{textperiodcentered}n}\varphi (x,{\mu }^2),\end{array}$$

$$⟨0|\mathrm{bar}d(0){\gamma }_0{\gamma }_{5}u(0)|{\pi }^{+}(p)⟩=i{f}_{\pi }{p}_0,$$

$$\xi =x-(1-x)=2x-1,$$

$$⟨{\xi }^n⟩={\int }_0^{1}dx(2x-1{)}^n\varphi (x,{\mu }^2).$$

$$\varphi (x,{\mu }^2)=6x(1-x)\left[1+\text{smashoperator}[]\sum _{n=1}^{\infty }{a}_n({\mu }^2)C_n^{3/2}(2x-1)\right],$$

$$\varphi (x,\mu \to \infty )={\varphi }^{\mathrm{as}}(x)=6x(1-x).$$

$$\mathcal{P}(x)$$

$${\mathcal{A}}_{\rho }(x)$$

$${\mathcal{O}}_{\rho \mu \nu }^{-}(x)$$

$${\mathcal{O}}_{\rho \mu \nu }^{+}(x)$$

$${\mathcal{O}}_{\rho \mu \nu }^{-}(x)$$

$${\mathcal{O}}_{\rho \mu \nu }^{+}(x)={\partial }_{(\mu }{\partial }_{\nu }{\mathcal{A}}_{\rho )}(x).$$

$${\mathcal{O}}_{4jk}^{\pm },j\ne k.$$

$$C_{\rho }(t,\mathbf{p})$$

$$C_{\rho \mu \nu }^{\pm }(t,\mathbf{p})$$

$$\displaystyle\begin{split}C_{\mathcal{O}}(t,\boldsymbol{p})&=\frac{1}{2E}\langle 0|\mathcal{O}(0)|\pi^{+}(p)\rangle\langle\pi^{+}(p)|J^{\dagger}(0)|0\rangle\\ &\quad\times\bigl{(}e^{-Et}+\tau_{\mathcal{O}}\tau_{J}e^{-E(T-t)}\bigr{)},\end{split}$$

$${\mathcal{R}}_{4ij}^{\pm }$$

$$⟨{\xi }^2{⟩}^{\overline{\text{MS}}}$$

$$a_2^{\overline{\text{MS}}}$$

$$(Fq{)}_{\mathbf{x}}=\sum _{\mathbf{y}}f(\mathbf{x}-\mathbf{y})q_{\mathbf{y}}.$$

$$f_{(\mathbf{k})}(\mathbf{x}-\mathbf{y})=f_{(0)}(0)\exp \left[-\frac{(\mathbf{x}-\mathbf{y}{)}^2}{2{\sigma }^2}+i\mathbf{k}(\mathbf{x}-\mathbf{y})\right],$$

$$\frac{\partial q(\tau )}{\partial \tau }=\alpha (\nabla -i\mathbf{k}{)}^{2}q(\tau ),$$

$$({\Phi }_{(\mathbf{k})}q{)}_{\mathbf{x}}=\frac{1}{1+6\epsilon }\left[q_{\mathbf{x}}+\epsilon \text{smashoperator}[]\sum _{j=\pm 1}^{\pm 3}{U}_{\mathbf{x},j}{e}^{-i\mathbf{k}\hat{\mathbf{ȷ}}}{q}_{\mathbf{x}+\hat{\mathbf{ȷ}}}\right].$$

$$U_{\mathbf{x},i}^{(m+1)}={\mathcal{P}}_{\text{SU}(3)}\left(\delta U_{\mathbf{x},i}^{(m)}+\text{smashoperator}[]\sum _{|j|\ne i}{U}_{\mathbf{x},i}^{(m)}{U}_{\mathbf{x}+\hat{\mathbf{ȷ}},i}^{(m)}{U}_{\mathbf{x}+\hat{\mathbf{ı}},j}^{(m)†}\right),$$

$$\sigma \approx \sqrt{2n{a}^2}\sqrt{\frac{\epsilon }{1+6\epsilon }},$$

$$\displaystyle S(t)\propto p_{i}p_{j}\exp\bigl{(}-\sqrt{m_{\pi}^{2}+\boldsymbol{p}^{2}}t\bigr{)}.$$

$$\displaystyle\boldsymbol{p}^{2}=\frac{2}{t^{2}}\bigl{(}1+\sqrt{1+m_{\pi}^{2}t^{2}}\bigr{)}.$$

$$⟨{\xi }^2⟩$$

$$a_2$$

$$\displaystyle\langle\xi^{2}\rangle^{\mkern 1.5mu\overline{\mkern-1.5mu\text{MS}\mkern-1.5mu}\mkern 1.5mu}($2\text{\,}\mathrm{GeV}$)$$

$$\displaystyle a_{2}^{\mkern 1.5mu\overline{\mkern-1.5mu\text{MS}\mkern-1.5mu}\mkern 1.5mu}($2\text{\,}\mathrm{GeV}$)$$

Full text

1 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\text{\,}\mathrm{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.

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 two-point 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].

2 General formalism

2.1 Continuum definitions

Pseudoscalar mesons like the pion have only one independent leading twist (twist two) DA, $\phi$, which is defined via a meson-to-vacuum matrix element of renormalized non-local quark-antiquark light-ray operators,

[TABLE]

where $z_{1,2}$ are real numbers, $n^{\mu}$ is an auxiliary light-like vector with $n^{2}=0$, and $|\pi^{+}(p)\rangle$ represents the ground state pseudoscalar $\pi^{+}$ meson with on-shell momentum $p^{2}=m_{\pi}^{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 $\phi$ is indicated by the argument $\mu^{2}$.

Neglecting both isospin breaking and electromagnetic effects, the DAs of the charged pseudoscalar $\pi^{\pm}$ and the neutral $\pi^{0}$ are trivially related such that it is sufficient to consider only one of them. The decay constant $f_{\pi}$ appearing in Eq. (1) can be obtained as the matrix element of a local operator,

[TABLE]

and has a value of $f_{\pi}\approx$130\text{,}\mathrm{MeV}$$ [14].

The physical interpretation of Eq. (1) is that the fraction $x$ of the pion momentum is carried by the $u$ quark, while the $\mathrm{bar}{d}$ antiquark carries the remaining fraction $1-x$. Hence the difference of the momentum fractions,

[TABLE]

contains all nontrivial information and its moments are defined as

[TABLE]

Since the Gegenbauer polynomials $C_{n}^{3/2}(2x-1)$, which correspond to irreducible representations of the collinear conformal group $\text{SL}(2,\mathds{R})$, form a complete set of functions, the DAs can be expanded as

[TABLE]

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 $\mu\rightarrow\infty$ only the leading term survives, which gives:

[TABLE]

2.2 Lattice definitions

From now on we will work in Euclidean spacetime and follow the conventions of Ref. [12]. The renormalized light-ray 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

[TABLE]

where $D_{\mu}$ 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 $(\dots)$, for example $\mathcal{O}_{(\mu\nu)}=\frac{1}{2}\left(\mathcal{O}_{\mu\nu}+\mathcal{O}_{\nu\mu}\right)-\frac{1}{4}\delta_{\mu\nu}\mathcal{O}_{\lambda\lambda}$. By using the shorthand notation \mathrlap{\reflectbox{!\vec{,\reflectbox{$D$}}}}\,\,\vec{\!\!D}_{\mu}=\vec{D}_{\mu}-\reflectbox{!\vec{,\reflectbox{$D$}}}_{\mu}, the operator $\mathcal{O}^{-}_{\rho\mu\nu}$ can also be written as

[TABLE]

The operator $\mathcal{O}^{+}_{\rho\mu\nu}$ is, in the continuum, given by the second derivative of the axial-vector current:

[TABLE]

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 $\mathrm{H}(4)$ [10, 9]. For our case, this corresponds to choosing all indices different for the operators $\mathcal{O}^{\pm}$. Identifying one index with the temporal direction, this leaves us with the operators

[TABLE]

In order to extract the desired moments we use two-point correlation functions of the operators $\mathcal{O}_{4jk}^{\pm}$ and $\mathcal{A}_{\rho}$ with an interpolating field,

[TABLE]

where $J=\mathcal{P}$ or $J=\mathcal{A}_{4}$. For sufficiently large $t$, the ground state dominates and the correlation functions give

[TABLE]

where the sign factors $\tau_{\mathcal{O}},\tau_{J}=\pm 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

[TABLE]

where $\langle\xi^{2}\rangle^{\rm{bare}}=R^{-}$ and $a^{\rm{bare}}_{2}=\frac{7}{12}(5R^{-}-R^{+})$. In our calculations we use the interpolator $J=\mathcal{P}$, as this gives a better overlap with the ground state than $\mathcal{A}_{4}$.

The renormalized moments in the $\mkern 1.5mu\overline{\mkern-1.5mu\text{MS}\mkern-1.5mu}\mkern 1.5mu$ scheme read

[TABLE]

where $\zeta_{ij}$ are ratios of renormalization constants which are defined in Ref. [12].

2.3 Momentum smearing

On a lattice of $N_{s}^{3}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\pi/L$. The calculation of the second moment of the DA requires a spatial momentum $\boldsymbol{p}=(2\pi/L)\boldsymbol{n}_{\boldsymbol{p}}$, with at least two non-vanishing components, i.e., $\boldsymbol{n}_{\boldsymbol{p}}^{2}\geq 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$:

[TABLE]

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 $\boldsymbol{k}=\boldsymbol{0}$. If the hadron carries a non-vanishing momentum $\boldsymbol{p}$, it is natural to assume that the quark will carry a momentum fraction $\boldsymbol{k}=\zeta\boldsymbol{p}$. We remark that there is no obvious relation between $\zeta$ and the longitudinal momentum fraction $x$ of the light-cone wave function. A Gaussian wave function with width $\sigma$ that is centred about the momentum $\boldsymbol{k}$ acquires a phase:

[TABLE]

where $f_{(\boldsymbol{0})}=f$. Our periodic lattice appears to imply a quantization of the possible values of $\boldsymbol{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_{(\boldsymbol{k})}q$ can be obtained as the result of evolving the heat equation with a drift term,

[TABLE]

starting from a spatial delta source at $\tau=0$, to the fictitious time $\tau=\sigma^{2}/(2\alpha)$.

One can approximate Eq. (22) by a discrete process, defining $F_{(\boldsymbol{k})}=\Phi_{(\boldsymbol{k})}^{n}$ as the $n$th application of an elementary iteration,

[TABLE]

In practice this smearing is implemented by multiplying the spatial connectors within the timeslice in question by the appropriate phases, $U_{\boldsymbol{x},j}\mapsto e^{-iak_{j}}U_{\boldsymbol{x},j}$. For $\boldsymbol{k}=\boldsymbol{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_{\boldsymbol{x},j}$ and $U_{\boldsymbol{x},-j}\equiv U_{\boldsymbol{x}-{\hat{\boldsymbol{\jmath}\mkern 4.0mu}\mkern-4.0mu}{},j}^{\dagger}$, where ${\hat{\boldsymbol{\jmath}\mkern 4.0mu}\mkern-4.0mu}{}$ denotes a vector of length $a$ and direction $j$, are spatially APE smeared [17]:

[TABLE]

where $i\in\{1,2,3\}$ and $j\in\{\pm 1,\pm 2,\pm 3\}$. The sum is over the four spatial “staples” surrounding $U_{\boldsymbol{x},i}$, and $\mathcal{P}_{\text{SU}(3)}$ is a gauge covariant projector onto the gauge group $\text{SU}(3)$, defined by maximizing $\operatorname{Re}\operatorname{Tr}\{A^{\dagger}\mathcal{P}_{\text{SU}(3)}(A)\}$. If the APE smeared links are close to unit fields then the width parameter of the resulting Gaussian is given as [6]222The root mean squared width of the resulting Gaussian will correspond to $\sqrt{3}\sigma$ as we have three spatial dimensions. This will shrink by a factor $1/\sqrt{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 $\sqrt{2}$ than the individual quark fields.

[TABLE]

where large values of $\varepsilon$ will allow for smaller iteration counts $n$, but the resulting function will be less smooth.

In the meson case the quark creation operator at the source needs to be smeared with $F_{(\boldsymbol{k})}$ and the quark destruction operator with $F_{(-\boldsymbol{k})}$, while for baryons all three quarks should be smeared with $F_{(\boldsymbol{k})}$, see Ref. [6] for details.333The 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^{-i\boldsymbol{p}\boldsymbol{x}}$ and $\boldsymbol{k}=\zeta\boldsymbol{p}$ with $\zeta\geq 0$. The phase used in Ref. [6] corresponds to the non-standard $e^{+i\boldsymbol{p}\boldsymbol{x}}$ convention that is used within the Chroma software suite [18].

3 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 3. These range from $N_{c}\approx 1500$ to $N_{c}\approx 2800$ configurations, separated by four hybrid Monte Carlo molecular dynamics units. The statistical errors were evaluated using the Bootstrap procedure, with $N_{\text{samples}}=500$, combined with the binning method. For the latter we have observed that a binsize of $N_{\text{bin}}=10$ saturates the statistical error.

The gauge links entering the quark smearing were APE smeared according to Eq. (24), employing $25$ iterations with the parameter $\delta=2.5$. We applied both, the standard Wuppertal smearing [15, 16] and the novel momentum smearing, i.e., we implemented Eq. (23) setting $\boldsymbol{k}=\boldsymbol{0}$ and $\boldsymbol{k}\neq\boldsymbol{0}$, respectively, and applied $300$ smearing steps with the smearing parameter $\varepsilon=0.25$. The root mean squared width of the squared pion interpolator wave function can be estimated using Eq. (25). This gives $\sqrt{3}\sigma\approx\sqrt{3}\mathbin{\mbox{\textperiodcentered}}$0.664\text{,}\mathrm{fm}$\approx$1.14\text{,}\mathrm{fm}$$.

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\approx$0.0857\text{,}\mathrm{fm}$$, we cannot as yet perform a continuum limit extrapolation.