RQCD Collaboration

# Mellin moments of spin dependent and independent PDFs of the pion and rho
meson

Marius Löffler<EMAIL_ADDRESS>Philipp Wein  Thomas Wurm  Simon
Weishäupl  Daniel Jenkins  Rudolf Rödl  Andreas Schäfer  Lisa Walter  Institut
für Theoretische Physik, Universität Regensburg, 93040 Regensburg, Germany

###### Abstract

We compute the second moments of pion and rho parton distribution functions
(PDFs) in lattice QCD with $N_{f}=2+1$ flavors of improved Wilson fermions. We
determine both singlet and non-singlet flavor combinations and, for the first
time, take disconnected contributions fully into account. In the case of the
rho, we also calculate the additional contribution arising from the $b_{1}$
structure function. The numerical analysis includes 26 ensembles, mainly
generated by the CLS effort, with pion masses ranging from
$420\text{\,}\mathrm{MeV}$ down to $214\text{\,}\mathrm{MeV}$ and with 5
different lattice spacings in the range of $0.1\text{\,}\mathrm{fm}$ to
$0.05\text{\,}\mathrm{fm}$. This enables us to take the continuum limit, as
well as to resolve the quark mass dependencies reliably. Additionally we
discuss the contaminations of rho correlation functions by two-pion states.

Lattice QCD, Parton distribution functions (PDFs), Deep inelastic scattering,
Photon interactions with hadrons

##  I Introduction

The pion is routinely investigated on the lattice, although, to the best of
our knowledge, disconnected quark contributions, e.g., to quark momentum
fractions, were not included so far. Since the pion is the pseudo-Goldstone
boson of dynamical chiral symmetry breaking (DCSB) its quark structure could
differ substantially from that of other mesons and if so, the flavor singlet
sea quark contribution is a natural place for such a difference to show up. In
contrast to the pion, the quark structure of the $\rho$ is only rarely
analyzed on the lattice, primarily due to the complications caused by its
resonance nature.

For the pion experimental data exists primarily from two classes of
experiments, namely Drell-Yan reactions with (secondary) pion beams, e.g.,
$\pi+N\to\mu^{+}+\mu^{-}+X$, which are sensitive to the pion PDF at large
$x\gtrsim 0.15$, and semi-inclusive (tagged) deep inelastic scattering (DIS),
e.g., $e+N\to e^{\prime}+N+X$, which is sensitive to small $x$ and exploits
the fact that the electron can scatter off the nucleon pion cloud via the
Sullivan process [1]. Experiments of the first type were performed by NA10
[2], E326 [3], E615 [4], and, more recently, by COMPASS [5]. This will be
continued by AMBER at CERN [6, 7]. Experiments of the second type were
performed at HERA [8, 9] (see also [10, 11, 12]) and are currently pursued at
JLab Hall A [13] (cf. the conditionally approved proposal [14]). They are also
under consideration for the physics program at the EIC [15].

In contrast to the pion case, there exists very little relevant experimental
data for the $\rho$ quark PDFs. The rho meson is the lightest strongly
decaying particle with a branching fraction of $>99.9\%$ into 2 pions [16]. It
is spin-1 which implies the existence of novel polarization dependent
structure functions [17]. The unstable nature of the $\rho$ complicates the
analysis of its structure, both, on the lattice – we will discuss some of the
implications in this article – and in experiment. However, as the goal of
hadron physics must be to also determine the quark-gluon structure of
resonances rather than only ground states, the $\rho$ is one of the most
attractive light mesons to explore. To the best of our knowledge, no existing
or planned experiment will investigate the spin structure of the $\rho$, and
so lattice calculations may offer the best, if not only, chance to determine
it. In ref. [18] it was speculated whether one could analyse the spin
structure of the $\rho$ in the meson cloud of a nucleon in a (polarized)
Sullivan process (see also ref. [10]), but the interpretation of such
measurements would be very non-trivial in view of the required analytic
continuation from the $t$ to the $s$ channel [19]. However, the $b_{1}$
structure function of the deuteron was measured by HERMES [20] (using DIS on
tensor-polarized deuteron gas with negligible vector polarization) and turned
out to be surprisingly large for such a loosly bound system. Also, while the
data for a limited $x$ range cannot really test the Close–Kumano sum rule for
the first moment of $b_{1}(x)$[21], an unexpected behavior outside of the
measured $x$ range is needed to fulfill it. Overall the results differ from
the expectation that the deuteron is in an $S$ wave with only a small $D$ wave
admixture (cf., e.g., ref. [22]). Also, there are efforts to measure the
deuteron $b_{1}$ via the proton-deuteron Drell-Yan process (conditionally
approved proposal at JLab Hall C [23] and feasibility studies for Fermilab
[22, 24]), as well as discussions of a measurement via DIS at the EIC [25].

The study of mesonic structure using lattice QCD has, by now, a history of
over three decades. Traditionally, such calculations focused on moments on
PDFs and distribution amplitudes (DAs). While earlier simulations [26, 27, 18,
28, 29, 30, 31] used quenched fermion representations, more recent simulations
[32, 33, 34, 35, 36, 37] use, for example, (clover-improved) dynamical Wilson
fermions, where the fermion determinant, and thus the quark sea, is taken into
account. However, what all these studies have in common is that they neglect
disconnected contributions, because the latter are notoriously difficult to
calculate and usually come with a large statistical error. During our analysis
we found that the noise on the light and strange quark disconnected loops is
highly correlated. We can use this to our advantage by looking at the non-
singlet ($\mathrm{bar}{u}u+\mathrm{bar}{d}d-2\mathrm{bar}{s}s$) and singlet
($\mathrm{bar}{u}u+\mathrm{bar}{d}d+\mathrm{bar}{s}s$) flavor combinations
instead of the light and strange loops themselves. While the large statistical
errors persist for the singlet flavor combination, they are reduced by over an
order of magnitude for the non-singlet operator, which allows us to obtain
quite precise results in this case even though we take the disconnected
contributions fully into account.

The reach of such calculations of moments of PDFs and DAs is limited,
primarily because for higher moments the problems caused by operator mixing
become untraceable. Therefore, in recent years ever more attention has focused
on coordinate space methods [38, 39, 40, 41, 42] which allow to calculate the
full functional forms of DAs and PDFs. As far as we know, no lattice results
for $b_{1}(x,Q^{2})$ of the $\rho$ have been published so far using these
methods, but some work exists for DAs of vector mesons [43], for which there
also exist results for some Mellin moments [44], and other mesons [45]. DAs
and PDFs probe independent aspects of the meson quark structure, and thus
provide valuable complementary information.

In this article we directly calculate the second moments of the pion and rho
PDFs by evaluating operators that contain a covariant derivative. The same
method would not be directly applicable to higher moments, since one would
face the problem of mixing with lower-dimensional operators. Sparked by the
presentation in ref. [39], position space methods have recently fueled a lot
of excitement, since they, in principle, allow for a resolution of the
complete PDF. There are recent studies on the pion PDF exploring possible
methods, such as the current-current method [46, 47], large momentum effective
theory [48, 49] (using quasi PDFs), or Ioffe time distributions [50] (using
pseudo PDFs). Similar to experiment and in contrast to the pion, the rho meson
structure has only been studied once to our knowledge, in the work presented
in ref. [18] (also discussed in refs. [28, 29]) based on a quenched simulation
at large quark masses. The use of large quark masses is probably due to the
additional difficulty raised by the instability of the rho meson at physical
quark masses. Actually, due to the finite volume, there is no continuum of
two-pion final states in a lattice simulation, and therefore the rho meson
cannot decay dynamically into two pions. Nevertheless, the discretized set of
two-pion (or, at higher energies, even multi-pion) finite volume states is
present and might overshadow the contribution from the rho in the correlation
function. We discuss this issue in some detail in Sec. III.4.

This article is structured as follows: We set the stage with a general
discussion of PDF properties and their connection to DIS structure functions
in Sec. II. Next, in Sec. III, we will present the details of the lattice
calculation, including simulation parameters, the analysis of correlation
functions, and possible two-pion contributions. We describe the extrapolation
strategy and our final results in Sec. IV, and summarize in Sec. V.

##  II General properties of PDFs

The cross-section of deep inelastic scattering can be written as a product of
a leptonic and a hadronic part. The hadronic tensor is given by

$\displaystyle W^{\mu\nu}(\bm{p},\lambda)$
$\displaystyle=\int\frac{\mathrm{d}^{4}z}{4\pi}\,\,e^{iq\cdot
z}\langle\bm{p},\lambda|\left[j^{\mu}(z),j^{\nu}(0)\right]|\bm{p},\lambda\rangle\,,$
(1)

where $\bm{p}$ is the three-momentum and $\lambda$ labels the spin of the
target hadron along a quantization direction [17, 51]. Using parity and time-
reversal invariance it is straight forward to show that the most general
hadronic tensor for polarized DIS from targets with spin-1 or less can be
decomposed into eight structure functions

$\displaystyle\begin{split}W^{\mu\nu}=&-F_{1}g_{\mu\nu}+F_{2}\frac{p_{\mu}p_{\nu}}{p\cdot
q}+i\frac{g_{1}}{p\cdot
q}\epsilon_{\mu\nu\lambda\sigma}q^{\lambda}s^{\sigma}\\\
&+i\frac{g_{2}}{(p\cdot
q)^{2}}\epsilon_{\mu\nu\lambda\sigma}q^{\lambda}\left(p\cdot
qs^{\sigma}-s\cdot qp^{\sigma}\right)\\\
&-b_{1}r_{\mu\nu}+\frac{1}{6}b_{2}\left(s_{\mu\nu}+t_{\mu\nu}+u_{\mu\nu}\right)\\\
&+\frac{1}{2}b_{3}\left(s_{\mu\nu}-u_{\mu\nu}\right)+\frac{1}{2}b_{4}\left(s_{\mu\nu}-t_{\mu\nu}\right)\,,\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{split}$ (2)

with the kinematic factors $r_{\mu\nu}$, $s_{\mu\nu}$, $t_{\mu\nu}$, and
$u_{\mu\nu}$, which depend on the momentum transfer $q$, the target momentum
$p$ and the target polarization vector $\epsilon$ (cf. App. A for our
conventions). The quantities $r_{\mu\nu}$, $s_{\mu\nu}$, $t_{\mu\nu}$, and
$u_{\mu\nu}$ are constructed such that they vanish upon averaging over the
target spin, see ref. [17]. For spin-$\frac{1}{2}$ targets $s^{\mu}$
corresponds to the spin four-vector. For spin-$1$ targets it corresponds to
$s^{\mu}=-i\epsilon^{\mu\nu\rho\sigma}e^{*}_{\nu}e_{\rho}p_{\sigma}$. Note
that, due to current conservation any term proportional to $q_{\mu}$ or
$q_{\nu}$ in Eq. (2) would vanish. Which of the structure functions can
contribute depends on the target spin: In case of spin-0 only $F_{1}$ and
$F_{2}$ do. For spin-$\frac{1}{2}$ targets one has $F_{1}$, $F_{2}$, $g_{1}$,
and $g_{2}$, where the measurement of $g_{1}$ and $g_{2}$ requires a polarized
beam. In case of spin-1 targets the full set of eight structure functions can
contribute. Notably, as argued in ref. [17], the additional structure
functions $b_{1-4}$ can be measured using an unpolarized electron beam.

The hadronic tensor can be factorized into a hard scattering kernel, which can
be calculated perturbatively, and in PDFs containing the nonperturbative
information. The PDFs related to the structure functions in Eq. (2) are
defined as111In general one finds three quark and three gluon PDFs analogous
to Eqs. (3) and (4), see, e.g., [52].

$\displaystyle q^{\lambda}_{H}$
$\displaystyle=\int\limits_{-\infty}^{\infty}\frac{dz^{-}}{4\pi}e^{-ixp^{+}z^{-}}\langle\bm{p},\lambda|\mathrm{bar}{q}(z)\gamma^{+}q(0)|\bm{p},\lambda\rangle\,,$
(3) $\displaystyle\Delta q^{\lambda}_{H}$
$\displaystyle=\int\limits_{-\infty}^{\infty}\frac{dz^{-}}{4\pi}e^{-ixp^{+}z^{-}}\langle\bm{p},\lambda|\mathrm{bar}{q}(z)\gamma^{+}\gamma_{5}q(0)|\bm{p},\lambda\rangle\,,$
(4)

where $z$ is a lightlike vector with vanishing plus component,
$z^{\mu}=z^{-}n^{\mu}_{-}$, where $n^{\mu}_{-}$ is dimensionless and can be
used to project vectors onto their plus component, e.g., $n_{-}\cdot p=p^{+}$.
PDF evolution with respect to the factorization scale, which delineates long-
from short-distance physics, is governed by the well-known DGLAP equations
[53, 54, 55]. To assure gauge invariance the fields in the nonlocal operators
are connected by Wilson lines which we do not write out explicitly. The PDF in
Eq. (3) corresponds to the sum $q=q_{\uparrow}+q_{\downarrow}$, while the PDF
in Eq. (4) corresponds to the difference $\Delta
q=q_{\uparrow}-q_{\downarrow}$ of the densities for quarks with opposite
helicity. For spin-$1$ hadrons symmetry implies that distributions for
different polarizations, $\lambda=+,0,-$, are related [17, 18]

$\displaystyle q^{+}$ $\displaystyle=q^{-}\,,$ $\displaystyle\Delta q^{+}$
$\displaystyle=-\Delta q^{-}\,,$ $\displaystyle\Delta q^{0}$
$\displaystyle=0\,,$ (5)

such that only three independent quark PDFs remain. The quark PDFs defined
above support $-1<x<1$, where the values at negative $x$ have to be
interpreted as momentum fractions of anti-quarks

$\displaystyle q^{\lambda}(x)$ $\displaystyle=-\mathrm{bar}q^{\lambda}(-x)\,,$
$\displaystyle\Delta q^{\lambda}(x)$
$\displaystyle=\Delta\mathrm{bar}q^{\lambda}(-x)\,,$ $\displaystyle\text{ for
$x<0$}\,.$ (6)

In order to see the connection between PDFs and structure functions let us
consider the operator product expansion

$\displaystyle\mathcal{O}_{a}(z)\,\mathcal{O}_{b}(0)$
$\displaystyle=\sum_{k}c_{abk}(z)\,\mathcal{O}_{k}(0)\,,$ (7)

which allows us to rewrite the product of two operators as a sum over local
operators assuming that the momentum components of the external states under
consideration are small compared to the inverse separation $1/z$. Using this
concept allows us to expand the product of the electromagnetic currents in Eq.
(1) into a series of local operators multiplied by coefficient functions
depending solely on the momentum transfer $q$. However, this is only valid for
target matrix elements provided that the momentum transfer $q$ is much larger
than the typical hadronic mass scale $\Lambda_{\rm{QCD}}$.

For any general operator $\mathcal{O}_{d,n}^{\mu_{1}\dots\mu_{n}}$ of
dimension $d$ and spin $n$ one can show that the terms in the expansion have
the structure

$\displaystyle c_{\mu_{1}\dots\mu_{n}}\mathcal{O}_{d,n}^{\mu_{1}\dots\mu_{n}}$
$\displaystyle\rightarrow\omega^{n}\left(\frac{Q}{M}\right)^{2-t}\,,$ (8)

where $M$ is the target hadron mass, $\omega=(2p\cdot q)/(-q^{2})$, and the
twist $t=d-n$. Taking into account that QCD operators contain at least two
quark fields ($t=1$ each) and an arbitrary number of covariant, symmetrized
derivatives
$\overleftrightarrow{D}^{\mu}=\overrightarrow{D}^{\mu}-\overleftarrow{D}^{\mu}$
($t=0$ each) a conventional basis for lowest twist $t=2$ quark operators can
be written in terms of two towers of operators222In general one finds six
towers of twist $t=2$ operators, see, e.g., [52].

$\displaystyle\mathcal{O}^{\mu_{1}\dots\mu_{n}}$
$\displaystyle=\frac{1}{2^{n-1}}\,\mathcal{S}\,\mathrm{bar}{q}\,\gamma^{\mu_{1}}i\overleftrightarrow{D}^{\mu_{2}}\dots
i\overleftrightarrow{D}^{\mu_{n}}\,q\,,$ (9)
$\displaystyle\mathcal{O}_{5}^{\mu_{1}\dots\mu_{n}}$
$\displaystyle=\frac{1}{2^{n-1}}\,\mathcal{S}\,\mathrm{bar}{q}\,\gamma^{\mu_{1}}\,\gamma_{5}i\overleftrightarrow{D}^{\mu_{2}}\dots
i\overleftrightarrow{D}^{\mu_{n}}\,q\,,$ (10)

where $\mathcal{S}$ projects out the completely symmetrized and traceless
components of the r.h.s. tensor. It is straight forward to confirm that the
matrix elements of these operators correspond to Mellin moments of the PDFs
[56],333See Eq. (33) in ref. [56] and note that the forward matrix elements in
Eqs. (3) and (4) are invariant under translation. e.g.,

$\displaystyle\begin{split}&\int_{-1}^{1}dx\,x^{n-1}q(x)=\langle
x^{n-1}\rangle_{q}+(-1)^{n}\langle x^{n-1}\rangle_{\mathrm{bar}q}\\\
&=\frac{1}{2p_{+}^{n}}n^{-}_{\mu_{1}}\cdots
n^{-}_{\mu_{n}}\langle\bm{p}|\mathcal{O}^{\mu_{1}\dots\mu_{n}}|\bm{p}\rangle\,,\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{split}$ (11)
$\displaystyle\begin{split}&\int_{-1}^{1}dx\,x^{n-1}\Delta q(x)=\langle
x^{n-1}\rangle_{\Delta q}-(-1)^{n}\langle
x^{n-1}\rangle_{\Delta\mathrm{bar}q}\\\
&=\frac{1}{2p_{+}^{n}}n^{-}_{\mu_{1}}\cdots
n^{-}_{\mu_{n}}\langle\bm{p}|\mathcal{O}^{\mu_{1}\dots\mu_{n}}_{5}|\bm{p}\rangle\,,\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{split}$ (12)

where we define the $n$-th moment of a function as

$\displaystyle M_{n}(f)$
$\displaystyle=\int_{0}^{1}\mathrm{d}x\,x^{n-1}\,f(x)=\langle
x^{n-1}\rangle_{f}\,.$ (13)

In perturbation theory and to leading twist accuracy the structure functions
are directly related to the PDFs, see, e.g., refs. [17, 18], where a generic
structure function $F$ is always obtained as the sum over the contributions
from quarks and antiquarks for the individual quark flavors weighted by the
square of their electric charge $e_{q}$:

$\displaystyle F$
$\displaystyle=\sum_{q}e_{q}^{2}\bigl{(}F^{q}+F^{\mathrm{bar}q}\bigr{)}\,.$
(14)

In the following we will only write down the quark contribution. The antiquark
contribution is obtained by simply substituting $q\rightarrow\mathrm{bar}q$.
For spin-$0$ targets one obtains

$\displaystyle F_{1}^{q}(x)$
$\displaystyle=\frac{1}{2}q(x)+\mathcal{O}(\alpha_{s})\,,$ (15) $\displaystyle
F_{2}^{q}(x)$ $\displaystyle=xq(x)+\mathcal{O}(\alpha_{s})\,,$ (16)

satisfying the Callan-Gross relation [57],
$F_{2}^{q}=2xF_{1}^{q}+\mathcal{O}(\alpha_{s})$. The gluon PDF does not appear
at leading order, since the gluons do not carry electric charge, and thus can
only couple through a quark loop. For spin-$1$ targets the hadronic tensor
depends on the hadron spin. Taking an average over the target spins one finds

$\displaystyle F_{1}^{q}(x)$
$\displaystyle=\frac{1}{6}\Bigl{(}q^{+}(x)+q^{0}(x)+q^{-}(x)\Bigr{)}+\mathcal{O}(\alpha_{s})\,,$
(17) $\displaystyle F_{2}^{q}(x)$
$\displaystyle=\frac{x}{3}\Bigl{(}q^{+}(x)+q^{0}(x)+q^{-}(x)\Bigr{)}+\mathcal{O}(\alpha_{s})\,.$
(18)

Considering the difference between targets with polarization $\lambda=\pm$ and
$\lambda=0$ one finds

$\displaystyle g_{1}^{q}(x)$ $\displaystyle=\frac{1}{2}\Delta
q^{+}(x)+\mathcal{O}(\alpha_{s})\,,$ (19) $\displaystyle b_{1}^{q}(x)$
$\displaystyle=\frac{1}{2}\Bigl{(}q^{0}(x)-q^{+}(x)\Bigr{)}+\mathcal{O}(\alpha_{s})\,,$
(20) $\displaystyle b_{2}^{q}(x)$
$\displaystyle=x\Bigl{(}q^{0}(x)-q^{+}(x)\Bigr{)}+\mathcal{O}(\alpha_{s})\,,$
(21)

which means that $b_{1}$ and $b_{2}$ are sensitive to a possible dependence of
the quark densities on the hadron polarization. The structure functions
$g_{2}$, $b_{3}$, and $b_{4}$ do not contribute at leading twist.

Next, we perform a Lorentz decomposition for the forward matrix elements of
the operators Eqs. (9) and (10). For a spin-0 particle this yields

$\displaystyle\langle\bm{p}|\mathcal{O}^{\mu_{1}\dots\mu_{n}}|\bm{p}\rangle=2\,\mathcal{S}\left[v^{q}_{n}p^{\mu_{1}}\cdots
p^{\mu_{n}}\right],$ (22)

with the so-called reduced matrix element $v_{n}$. Operators containing
$\gamma_{5}$ do not contribute because of symmetry relations. For a spin-$1$
particle we find three independent structures

$\displaystyle\langle\bm{p},\lambda|\mathcal{O}^{\mu_{1}\dots\mu_{n}}|\bm{p},\lambda\rangle$
$\displaystyle=2\,\mathcal{S}\left[a^{q}_{n}p^{\mu_{1}}\cdots
p^{\mu_{n}}+d^{q}_{n}\left(m^{2}\epsilon^{*\mu_{1}}(\bm{p},\lambda)\,\epsilon^{\mu_{2}}(\bm{p},\lambda)-\frac{1}{3}p^{\mu_{1}}p^{\mu_{2}}\right)p^{\mu_{3}}\cdots
p^{\mu_{n}}\right]\,,$ (23)
$\displaystyle\langle\bm{p},\lambda|\mathcal{O}^{\mu_{1}\dots\mu_{n}}_{5}|\bm{p},\lambda\rangle$
$\displaystyle=2i\,\mathcal{S}\left[r^{q}_{n}\epsilon^{\rho\sigma\tau\mu_{1}}\epsilon^{*}_{\rho}(\bm{p},\lambda)\epsilon_{\sigma}(\bm{p},\lambda)p_{\tau}\,p^{\mu_{2}}\cdots
p^{\mu_{n}}\right]\,,$ (24)

where we use the convention that $\varepsilon^{0123}=-1$. Here, $a_{n}$ is
related to the polarization averaged contribution and $d_{n}$ to the polarized
contribution of the quark PDF $q$. The reduced matrix element $r_{n}$ in Eq.
(24) is related to the (quark-)spin dependent PDF $\Delta q$.

In the structure functions always the sum of quark and antiquark contributions
$F^{q+\mathrm{bar}q}=F^{q}+F^{\mathrm{bar}q}$ is relevant, see Eq. (14).
Comparing this to Eqs. (11) and (12) one notices that the matrix elements
given above yield information about either the even or (exclusive) the odd
moments of a given structure function. For spin-$0$ targets one finds

$\displaystyle\begin{split}2M_{n}(F_{1}^{q+\mathrm{bar}q})&=C_{n}^{(1)}v^{q}_{n}\,,\qquad
n\text{ even}\,,\\\
M_{n-1}(F_{2}^{q+\mathrm{bar}q})&=C_{n}^{(2)}v^{q}_{n}\,,\qquad n\text{
even}\,,\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{split}$ (25)

while, for spin-$1$ targets,

$\displaystyle\begin{split}2M_{n}(F_{1}^{q+\mathrm{bar}q})&=C_{n}^{(1)}a^{q}_{n}\,,\qquad
n\text{ even}\,,\\\
M_{n-1}(F_{2}^{q+\mathrm{bar}q})&=C_{n}^{(2)}a^{q}_{n}\,,\qquad n\text{
even}\,,\\\ 2M_{n}(b_{1}^{q+\mathrm{bar}q})&=C_{n}^{(1)}d^{q}_{n}\,,\qquad
n\text{ even}\,,\\\
M_{n-1}(b_{2}^{q+\mathrm{bar}q})&=C_{n}^{(2)}d^{q}_{n}\,,\qquad n\text{
even}\,,\\\ 2M_{n}(g_{1}^{q+\mathrm{bar}q})&=C_{n}^{(3)}r^{q}_{n}\,,\qquad
n\text{ odd}\,.\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{split}$ (26)

The $C^{(k)}_{n}=1+\mathcal{O}(\alpha_{s})$ are the Wilson coefficients of the
OPE.

We can also relate the moments of the PDFs to the reduced matrix elements. By
substituting Eq. (22) into Eq. (11), we find for spin-$0$

$\displaystyle\begin{split}v^{q}_{n}&=\langle
x^{n-1}\rangle_{q+\mathrm{bar}q}\,,\qquad n\text{ even}\,,\\\
v^{q}_{n}&=\langle x^{n-1}\rangle_{q-\mathrm{bar}q}\,,\qquad n\text{
odd}\,.\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{split}$ (27)

For spin-$1$ hadrons we find (by substituting Eq. (23) into Eq. (11)) that

$\displaystyle\begin{split}a^{q}_{n}&=\frac{1}{3}\sum_{\lambda=\pm,0}\langle
x^{n-1}\rangle_{q^{\lambda}+\mathrm{bar}q^{\lambda}}\,,\qquad n\text{
even}\,,\\\ a^{q}_{n}&=\frac{1}{3}\sum_{\lambda=\pm,0}\langle
x^{n-1}\rangle_{q^{\lambda}-\mathrm{bar}q^{\lambda}}\,,\qquad n\text{
odd}\,,\\\ d^{q}_{n}&=\langle x^{n-1}\rangle_{q^{0}+\mathrm{bar}q^{0}}-\langle
x^{n-1}\rangle_{q^{+}+\mathrm{bar}q^{+}}\,,\qquad n\text{ even}\,,\\\
d^{q}_{n}&=\langle x^{n-1}\rangle_{q^{0}-\mathrm{bar}q^{0}}-\langle
x^{n-1}\rangle_{q^{+}-\mathrm{bar}q^{+}}\,,\qquad n\text{
odd}\,,\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{split}$ (28)

i.e., $a^{q}_{n}$ yields the polarization average, while $d^{q}_{n}$
corresponds to the difference between hadrons with polarization $\lambda=\pm$
and $\lambda=0$. In the following we will be particularly interested in the
second moments, since the corresponding operator (cf. Eq. (9) with $n=2$) is
equivalent to the quark part of the energy-momentum tensor [58], and describes
the distribution of the momentum within the hadron. For instance, in the
spin-$1$ case a non-zero value of $d^{q}_{2}$ would indicate that the portion
of the momentum carried by quarks of flavor $q$ depends on the polarization
direction of the hadron.

In analogy to the relations in Eqs. (27) and (28) the reduced matrix element
$r_{n}$ in Eq. (24) can be related to the moments of $\Delta q$ defined in Eq.
(12). We will, however, restrict ourselves to the computation of the second
moments of the vector PDF Eq. (3) for the rest of this work.

##  III Computation on the lattice

###  III.1 Lattice setup and numerical methods

To calculate the second moment of the structure functions introduced in the
last section we analyzed a subset of the lattice gauge ensembles generated
within the Coordinated Lattice Simulations (CLS) effort [59]. The ensembles
have been generated using a tree-level Symanzik improved gauge action with
$N_{f}=2+1$ flavors of nonperturbatively order $a$ improved Wilson (clover)
fermions. Stable Monte-Carlo sampling is achieved by applying twisted-mass
determinant reweighting [60] to avoid near zero modes of the Wilson-Dirac
operator.

To avoid freezing of the topological charge and large autocorrelation times
for the very fine lattices we use open boundary conditions for most of our
simulations [60, 61]. Only some of the coarser lattices are simulated using
periodic boundaries. The CLS gauge ensembles are generated along three
different trajectories in the renormalized quark mass plane

  * •

${\rm TrM}=\text{const.}$: The trace of the quark mass matrix is kept constant
near its physical value [59]

  * •

$m_{s}=\text{const.}$: The strange quark mass is kept constant close to its
physical value [62]

  * •

$m_{l}=m_{s}$: The symmetric line.

This strategy is explained in [62] while an additional graphical illustration
can be found in [63]. A complete list of the gauge ensembles used in this work
is shown in Tab. 1. We use five different lattice spacings from
$0.0497\text{\,}\mathrm{fm}$ up to $0.0984\text{\,}\mathrm{fm}$ and $m_{\pi}$
covers a range from $\sim$420\text{\,}\mathrm{MeV}$$ down to
$\sim$220\text{\,}\mathrm{MeV}$$ with volumes $Lm_{\pi}$ between $3.8$ and
$6.4$, see Tab. 1.

Table 1: CLS gauge ensembles analyzed in this work labeled by their identifier and sorted by the inverse coupling $\beta$ and $m_{\pi}$. We also label lattice volume in spatial and temporal direction, the boundary condition in time, open (o) or periodic (p), the lattice spacing, the pion mass, the volume in terms of $Lm_{\pi}$, and the rho meson mass $m_{\rho}$ computed in Sec. IV.1. The list of source-sink distances analyzed for the connected three-point function of the ensemble is labeled by $t$ and if more than one measurement for a set of source-sink separations was performed we denote the corresponding number of measurements in parenthesis. In physical units these distances roughly correspond to $[0.7,0.9,1.0,1.2]\text{\,}\mathrm{fm}$. The number of configurations analyzed for the ensemble is denoted as $n_{\rm{cnfgs}}$ and traj. specifies the trajectory in the quark mass plane of the ensemble [62]. An in-depth description of the ensemble generation can be found in [59]. Ens. | $\beta$ | $N_{s}\times N_{t}$ | bc |  $a$[fm] |  $m_{\pi}$ [MeV] | $Lm_{\pi}$ |  $m_{\rho}$ [MeV] | $t/a$ | $n_{\rm{cnfgs}}$ | traj.
---|---|---|---|---|---|---|---|---|---|---
A653 | 3.34 | $24\times 48$ | p | 0.0984 | 426 | 5.1 | 870 | [7, 9, 11, 13] | 2525 | sym
A650 | 3.34 | $24\times 48$ | p | 0.0984 | 368 | 4.4 | 813 | [7, 9, 11, 13] | 2216 | sym
H101 | 3.4 | $32\times 96$ | o | 0.0859 | 420 | 5.9 | 860 | [8, 10, 12, 14] | 2016 | trm/sym
H102r001 | 3.4 | $32\times 96$ | o | 0.0859 | 352 | 4.9 | 828 | [8, 10, 12, 14] (2) | 1992 | trm
H102r002 | 3.4 | $32\times 96$ | o | 0.0859 | 356 | 5.0 | 820 | [8, 10, 12, 14] (2) | 2016 | trm
H105 | 3.4 | $32\times 96$ | o | 0.0859 | 279 | 3.9 | 793 | [8, 10, 12, 14] | 2052 | trm
H106 | 3.4 | $32\times 96$ | o | 0.0859 | 272 | 3.8 | 805 | [8, 10, 12, 14] | 1543 | ms
H107 | 3.4 | $32\times 96$ | o | 0.0859 | 366 | 5.1 | 860 | [8, 10, 12, 14] | 1564 | ms
C101 | 3.4 | $48\times 96$ | o | 0.0859 | 220 | 4.6 | 753 | [8, 10, 12, 14] | 1997 | trm
C102 | 3.4 | $48\times 96$ | o | 0.0859 | 222 | 4.6 | 756 | [8, 10, 12, 14] | 1465 | ms
B450 | 3.46 | $32\times 64$ | p | 0.0760 | 418 | 5.2 | 869 | [9, 11, 14, 16] | 1612 | trm/sym
B452 | 3.46 | $32\times 64$ | p | 0.0760 | 350 | 4.3 | 856 | [9, 11, 14, 16] | 1944 | ms
D450 | 3.46 | $64\times 128$ | p | 0.0760 | 214 | 5.3 | 741 | [9, 11, 14, 16] | 617 | trm
N450 | 3.46 | $48\times 128$ | p | 0.0760 | 285 | 5.3 | 812 | [9, 11, 14, 16] | 1131 | ms
S400 | 3.46 | $32\times 128$ | o | 0.0760 | 352 | 4.3 | 839 | [9, 11, 14, 16] | 2001 | ms
X450 | 3.46 | $48\times 64$ | p | 0.0760 | 263 | 4.9 | 737 | [9, 11, 14, 16] | 400 | sym
rqcd030 | 3.46 | $32\times 64$ | p | 0.0760 | 317 | 3.9 | 795 | [9, 11, 14, 16] | 1222 | sym
N201 | 3.55 | $48\times 128$ | o | 0.0643 | 285 | 4.5 | 822 | [11, 14, 16, 19] | 1522 | ms
N202 | 3.55 | $48\times 128$ | o | 0.0643 | 411 | 6.4 | 860 | [11, 14, 16, 19] | 899 | trm/sym
N203 | 3.55 | $48\times 128$ | o | 0.0643 | 345 | 5.4 | 833 | [11, 14, 16, 19] (2) | 3086 | trm
N204 | 3.55 | $48\times 128$ | o | 0.0643 | 351 | 5.5 | 859 | [11, 14, 16, 19] | 1500 | ms
N200 | 3.55 | $48\times 128$ | o | 0.0643 | 284 | 4.4 | 804 | [11, 14, 16, 19] (2) | 3424 | trm
X250 | 3.55 | $48\times 64$ | p | 0.0643 | 348 | 5.4 | 816 | [11, 14, 16, 19] | 345 | sym
X251 | 3.55 | $48\times 64$ | p | 0.0643 | 267 | 4.2 | 757 | [11, 14, 16, 19] | 434 | sym
J303 | 3.7 | $64\times 192$ | o | 0.0497 | 257 | 4.2 | 802 | [14, 17, 21, 24] | 1068 | trm
N300 | 3.7 | $48\times 128$ | o | 0.0497 | 422 | 5.1 | 891 | [14, 17, 21, 24] | 1539 | trm/sym
N304 | 3.7 | $48\times 128$ | o | 0.0497 | 351 | 4.3 | 884 | [14, 17, 21, 24] | 1726 | ms

The two- and three-point functions introduced in Sec. III.3 are computed on
the lattice using the gauge configurations in Tab. 1. While we get the two-
point functions by an inversion of the lattice Dirac operator using common
numerical solvers (in particular we use a modified version of the Wuppertal
adaptive algebraic multigrid code DD-$\alpha$AMG [64, 65] on SIMD
architectures [66, 67, 68, 69] and the IDFLS solver [70, 71] on other
architectures) the computation of the three-point functions is more involved.
The three-point function connected parts of all ensembles are computed using
stochastic estimators as described in App. C.1. The computation of the three-
point function disconnected contributions is described in App. C.2. To improve
the overlap of the interpolating currents at the source and the sink timeslice
we use Wuppertal smeared [72] quarks in the source and sink interpolators
employing APE-smoothed gauge links [73]. All the computations are performed
using the Chroma software package [74] and additional libraries implemented by
our group.

###  III.2 Renormalization

In order to obtain physically meaningful results, the bare operators
introduced in Eqs. (9) and (10) have to be renormalized. In this context, one
faces the additional difficulty that the isosinglet quark operators will mix
under renormalization with the gluonic operators, schematically,

$\displaystyle\mathcal{O}^{\rm ren}$
$\displaystyle=Z^{qq}\mathcal{O}+Z^{qg}\mathcal{O}_{g}\,,$ (29)

where we have suppressed the Lorentz indices for better readability. For the
opposite direction (i.e., admixture of quark operators into gluon operators)
it has been shown in ref. [75], using one-loop perturbative renormalization,
that this admixture is a few percent effect (see also the discussion in ref.
[76]). We will assume that the same is true for the admixture of gluonic
operators into quark operators, and that, as a consequence, its effect is
negligible within the statistical accuracy of this work. Still, this caveat
has to be kept in mind and needs to be addressed in future work. Note,
however, that operators without an isosinglet part (e.g., with flavor
structure $\mathrm{bar}uu+\mathrm{bar}dd-2\mathrm{bar}ss$) are not affected.
Furthermore, we will approximate the isosinglet renormalization factor by the
(non-perturbatively calculated) renormalization factor for isovector currents.
This is exact to NLO accuracy (within a perturbative renormalization
procedure).

On the lattice, the continuous Euclidean $\mathrm{O}(4)$ symmetry is reduced
to that of its finite hypercubic subgroup $\mathrm{H}(4)$. Therefore, symmetry
imposes much weaker constraints on the mixing of operators under
renormalization. In order to avoid mixing as far as possible, in particular
mixing with lower-dimensional operators, we use operators from suitably chosen
multiplets that possess a definite C-parity and transform according to
irreducible representations of $\mathrm{H}(4)$, cf. refs. [77, 18]. To be
specific, we will use the operators $\mathcal{O}_{\rm v2a}=\mathcal{O}^{0i}$
and $\mathcal{O}_{\rm v2b}=\frac{4}{3}\mathcal{O}^{00}$, cf. App. B, where
also an explicit definition of the operators is provided.

Our final results will be given in the $\mkern
1.5mu\overline{\mkern-1.5mu\text{MS}\mkern-1.5mu}\mkern 1.5mu$ scheme at a
scale of $2\text{\,}\mathrm{GeV}$. To this end, we adopt a two step procedure:
First, we calculate the renormalization factors nonperturbatively in the
$\text{RI}^{\prime}\mkern-3.0mu\raisebox{1.0pt}{$/$}\text{SMOM}$ scheme. These
are then converted to the $\mkern
1.5mu\overline{\mkern-1.5mu\text{MS}\mkern-1.5mu}\mkern 1.5mu$ scheme using
perturbative QCD. The whole procedure is described in great detail in ref.
[78], including subtleties due to the use of open boundaries in time direction
and details of the perturbative subtraction of lattice artifacts. To be
specific, we use the values for $Z_{\rm v2a}$ and $Z_{\rm v2b}$ based on
$\text{RI}^{\prime}\mkern-3.0mu\raisebox{1.0pt}{$/$}\text{SMOM}$ as
intermediate scheme, with the perturbative subtraction of lattice artifacts
(without the use of the so-called fixed-scale method). The explicit values of
$Z^{qq}$ for the operator combinations $\mathcal{O}_{\rm v2a}$ and
$\mathcal{O}_{\rm v2b}$ used in this work are given in Tab. 2.

###  III.3 Correlation functions

In order to calculate the DIS structure functions on the lattice one has to
compute two- and three-point correlation functions in the forward limit:

$\displaystyle C_{{\rm 2},\bm{p},t}^{(\mu\nu)}$
$\displaystyle=a^{3}\sum_{\bm{x}}e^{-i\bm{p}\cdot\bm{x}}\langle
O_{M}^{(\mu)}(\bm{x},t)\bar{O}_{M}^{(\nu)}(\bm{0},0)\rangle\,,$ (30)
$\displaystyle C_{{\rm 3},\bm{p},t,\tau}^{(\mu\nu)}$
$\displaystyle=a^{6}\mathop{\smash{\sum_{\bm{x},\bm{y}}}}e^{-i\bm{p}\cdot\bm{x}}\langle
O_{M}^{(\mu)}(\bm{x},t)\mathcal{O}(\bm{y},\tau)\bar{O}_{M}^{(\nu)}(\bm{0},0)\rangle\,.$
(31)

We will consider pion ($M=\pi$) and rho mesons ($M=\rho$), where the Lorentz
indices are only necessary in the latter case. The interpolating current
$\smash{\bar{O}_{M}^{(\mu)}}$ creates a meson state with matching quantum
numbers at the source timeslice $t_{\rm{src}}$ while $\smash{O_{M}^{(\nu)}}$
annihilates the meson at the sink time slice $t_{\rm{snk}}$. They read

$\displaystyle O_{\pi}$ $\displaystyle=\mathrm{bar}q_{f}\gamma_{5}q_{g}\,,$
$\displaystyle O_{\rho}^{i}$
$\displaystyle=\mathrm{bar}q_{f}\gamma^{i}q_{g}\,,$ (32)

with appropriately chosen quark flavors $f$ and $g$ and $i=1,2,3$. The quark
fields in the interpolating currents are spatially smeared (see Sec. III.1) to
enhance the ground state overlap. In addition to the two interpolating
currents the three-point function contains an insertion current $\mathcal{O}$
at timeslice $\tau$ with $0<\tau<t$. The extraction of the ground-state matrix
element of $\mathcal{O}$ is the key task in the subsequent calculations. In
this work we set $t=t_{\rm{snk}}-t_{\rm{src}}$,
$\tau=t_{\rm{ins}}-t_{\rm{src}}$ and hence $t_{\rm{src}}=0$ without loss of
generality.

####  III.3.1 The pion

For the pion case we first define the matrix elements

$\displaystyle\langle 0|O_{\pi}|\bm{p}\rangle$
$\displaystyle=\sqrt{\vphantom{Z_{\bm{p}}}\smash{Z^{\pi}_{\bm{p}}}}\,,$ (33)

where $Z^{\pi}_{\bm{p}}$ is smearing dependent and encodes the overlap of the
ground state with the interpolating currents at the source and the sink.
Inserting a complete set of states into Eq. (30) allows us to expand the two-
point correlation function in terms of hadronic matrix elements. At large
Euclidian times the correlation function can be approximated by the ground
state contribution

$\displaystyle C_{{\rm 2},\bm{p},t}$
$\displaystyle=Z^{\pi}_{\bm{p}}\frac{e^{-E^{\pi}_{\bm{p}}t}}{2E^{\pi}_{\bm{p}}}+\dots\,,$
(34)

where we assume that the same smearing setup is used at the source and the
sink. For the ground state energies $E^{\pi}_{\bm{p}}$ we impose the continuum
dispersion relation ${E^{\pi}_{\bm{p}}=\sqrt{m_{\pi}^{2}+\bm{p}^{2}}}$.

Similarly, one can show that the spectral decomposition of the three-point
function in Eq. (31) for large Euclidean times reads

$\displaystyle\begin{split}C_{{\rm
3},\bm{p},t,\tau}&=Z^{\pi}_{\bm{p}}\frac{e^{-E^{\pi}_{\bm{p}}t}}{(2E^{\pi}_{\bm{p}})^{2}}\langle\bm{p}|\mathcal{O}|\bm{p}\rangle+\dots\,.\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{split}$ (35)

In practice it turns out that especially for the three-point functions the
signal-to-noise ratio at large Euclidean time distances $t$ and $\tau$ does
not only contain the ground state contribution. How to exclude further excited
state contributions is explained in Sec. III.3.3.

####  III.3.2 The rho meson

For the rho mesons we define, in analogy to the pseudoscalar case Eq. (33),

$\displaystyle\langle 0|O_{\rho}^{\mu}|\bm{p},\lambda\rangle$
$\displaystyle=\sqrt{\vphantom{Z_{\bm{p}}}\smash{Z^{\rho}_{\bm{p}}}}\epsilon^{\mu}(\bm{p},\lambda)\,,$
(36)

where the polarization vector $\epsilon^{\mu}$ obeys the general
transversality condition Eq. (75). Therefore, the insertion of a complete set
of states into Eq. (30) (including a sum over all possible polarizations),
yields

$\displaystyle C_{{\rm 2},\bm{p},t}^{\mu\nu}$
$\displaystyle=-Z^{\rho}_{\bm{p}}\frac{e^{-E^{\rho}_{\bm{p}}t}}{2E^{\rho}_{\bm{p}}}\left(g^{\mu\nu}-\frac{p^{\mu}p^{\nu}}{m_{\rho}^{2}}\right)+\dots\,,$
(37)

where we have only written out the contribution from the rho meson, which is
the leading one-particle state at large Euclidean times. Note, however, that
in this case also two-pion states can occur, which, depending on the
simulation parameters, can have smaller energies than the rho meson. This
problem will be discussed in Sec. III.4.

Inserting two complete sets of states into the three-point function Eq. (31)
we find

$\displaystyle C_{{\rm 3},\bm{p},t,\tau}^{\mu\nu}$
$\displaystyle=Z^{\rho}_{\bm{p}}\frac{e^{-E^{\rho}_{\bm{p}}t}}{(2E^{\rho}_{\bm{p}})^{2}}\mathop{\smash{\sum_{\lambda^{\prime},\lambda}}}\epsilon^{\mu}(\bm{p},\lambda^{\prime})\epsilon^{\nu*}(\bm{p},\lambda)\,\langle\bm{p},\lambda^{\prime}|\mathcal{O}|\bm{p},\lambda\rangle\,.$
(38)

Excited state contributions will be treated in Sec. III.3.3.

Table 2: Renormalization factors $Z^{qq}$ for the operator combinations $\mathcal{O}_{\rm v2a}$ and $\mathcal{O}_{\rm v2b}$ used in this work. All values taken from Tab. XIII in [78]. $\beta$ | 3.34 | 3.4 | 3.46 | 3.55 | 3.7
---|---|---|---|---|---
$Z^{qq}_{\rm v2a}$ | 1.0731 | 1.1010 | 1.1251 | 1.1578 | 1.2053
$Z^{qq}_{\rm v2b}$ | 1.0672 | 1.0938 | 1.1170 | 1.1485 | 1.1949

####  III.3.3 Excited states analysis

In the three-point functions Eqs. (34) and (37) the signal-to-noise ratio
decreases exponentially with the source-sink separation in time. At small time
distances between the operators, however, there are still noticeable excited
state effects. We take these into account by allowing for a generic excited
state contribution in the spectral decomposition of the correlation functions.
For the pseudoscalar correlation functions Eqs. (34) and (35) our ansatz reads

$\displaystyle C_{{\rm 2},\bm{p},t}$
$\displaystyle=Z^{\pi}_{\bm{p}}\frac{e^{-E^{\pi}_{\bm{p}}t}}{2E^{\pi}_{\bm{p}}}\Bigl{(}1+Ae^{-\Delta
E_{\bm{p}}t}\Bigr{)}\,,$ (39) $\displaystyle\begin{split}C_{{\rm
3},\bm{p},t,\tau}&=Z^{\pi}_{\bm{p}}\frac{e^{-E^{\pi}_{\bm{p}}t}}{(2E^{\pi}_{\bm{p}})^{2}}\langle\bm{p}|\mathcal{O}|\bm{p}\rangle\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\\\
&\times\Bigl{(}1+B_{10}e^{-\Delta E_{\bm{p}}(t-\tau)}+B_{01}e^{-\Delta
E_{\bm{p}}\tau}+B_{11}e^{-\Delta E_{\bm{p}}t}\Bigr{)}\,,\end{split}$ (40)

where $\Delta E_{\bm{p}}$ denotes the energy difference to the first excited
state. The excited state amplitude in the two-point function, $A$, depends on
the interpolating currents at the source and the sink, their smearing, and the
momentum $\bm{p}$, while the amplitudes in the three-point function, $B_{10}$,
$B_{01}$, and $B_{11}$, also depend on the operator insertion $\mathcal{O}$.

For the rho meson case we perform the analysis analogously. However, in
particular for ensembles with small quark masses and large volumes, one would
in this situation expect a contribution from (possibly multiple) two-pion
states, which can have even smaller energy than the “ground state” rho meson
itself. Despite the fact that we do not find any trace of these two-pion
states in our numerical analysis, we cannot claim to have this problem fully
under control; cf. the discussion of this delicate issue in Sec. III.4.

####  III.3.4 Ratios

Instead of performing a fit to three-point functions, one can equivalently fit
to ratios of two- and three-point functions. As discussed in ref. [63], this
can be advantageous due to a cancellation of unwanted correlations between
two- and three-point functions. Furthermore, the ratio can be chosen in such a
way that contributions from the ground state directly corresponds to the
matrix element we are interested in. For the pseudoscalar correlation
functions we define

$\displaystyle R_{\bm{p}}$
$\displaystyle=\frac{C_{3,\bm{p},t,\tau}}{C_{2,\bm{p},t}}\xrightarrow{t\gg\tau\gg
0}\frac{\langle\bm{p}|\mathcal{O}|\bm{p}\rangle}{2E^{\pi}_{\bm{p}}}\,,$ (41)

which holds for any operator insertion $\mathcal{O}$ in the three-point
function. For the vector meson case we will consider the diagonal case with
the same Lorentz indices at the sink and at the source (i.e., $\mu=\nu=i$ in
Eqs. (30) and (31)). Defining
$J^{\bm{p}}_{\lambda^{\prime}\lambda}\equiv\langle\bm{p},\lambda^{\prime}|\mathcal{O}|\bm{p},\lambda\rangle/(2E^{\rho}_{\bm{p}})$,
one obtains

$\displaystyle\begin{split}&R^{i}_{\bm{p}}=\frac{C^{ii}_{3,\bm{p},t,\tau}}{C^{ii}_{2,\bm{p},t}}\\\
&\xrightarrow{t\gg\tau\gg
0}\frac{m_{\rho}^{2}}{\left(E^{\rho}_{\bm{p}}\right)^{2}}\sum_{\lambda,\lambda^{\prime}}\epsilon^{i}(\bm{p},\lambda^{\prime})\epsilon^{i*}(\bm{p},\lambda)J^{\bm{p}}_{\lambda^{\prime}\lambda}\,,\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{split}$ (42)

where $i$ is fixed (no summation). On the right-hand side a sum over multiple
matrix elements occurs, which can be evaluated explicitly for the chosen
three-momentum. For on-axis momenta $\hat{\bm{p}}=\pm\bm{e}_{i}$ one finds the
simple formulas

$\displaystyle J^{\bm{p}}_{00}$ $\displaystyle=R^{i}_{\bm{p}}\,,$
$\displaystyle J^{\bm{p}}_{++}+J^{\bm{p}}_{--}$ $\displaystyle=\sum_{j\neq
i}R^{j}_{\bm{p}}\,.$ (43)

for the extraction of the polarization-conserving matrix elements.

###  III.4 Two-pion state contribution in the vector meson case

In an infinite volume, above the particle creation threshold, a continuum of
states would contribute to the spectral decomposition of the rho meson. In
particular in Eq. (37), a continuum of two-pion states would contribute above
the $2m_{\pi}$ threshold. In the non-interacting case, their center of mass
energies are given by

$\displaystyle E_{\rm cm}$ $\displaystyle=2\sqrt{m_{\pi}^{2}+\bm{k}^{2}}\,,$
(44)

where $\bm{k}$ and $-\bm{k}$ are the momenta of the two pions in the center of
mass frame. In a finite volume momenta are quantized such that one gets a sum
over a discrete set of states that contribute.

Table 3: Little groups and decomposition of angular momentum $1$ in irreducible representations for all momentum sectors $\bm{d}^{2}\leq 4$, where $\bm{p}=\frac{2\pi}{L}\bm{d}$. The groups are isomorphic for each representative of a sector. Table taken from ref. [79]. $\bm{d}^{2}$ | $\rm{LG}(\bm{p})$ | $\Gamma$
---|---|---
$0$ | $O_{h}$ | $\mathrm{T{1u}}$
$1$ | $C_{4v}$ | $\mathrm{A1}\oplus\mathrm{E}$
$2$ | $C_{2v}$ | $\mathrm{A1}\oplus\mathrm{B1}\oplus\mathrm{B2}$
$3$ | $C_{3v}$ | $\mathrm{A1}\oplus\mathrm{E}$
$4$ | $C_{4v}$ | $\mathrm{A1}\oplus\mathrm{E}$
Table 4: Scattering phase shifts (assuming that only the $P$-wave contributes) for momentum sectors $\bm{d}^{2}$ and irreducible representations $\Gamma$. See refs. [79, 80] for more details. $\bm{d}^{2}$ | $\Gamma$ | $\phi^{\bm{d}}_{\Gamma}$
---|---|---
$0$ | $\mathrm{T1u}$ | $w_{0,0}-w_{2,0}-\sqrt{6}w_{2,2}$
$1$ | $\mathrm{A1}$ | $w_{0,0}+2w_{2,0}$
$1$ | $\mathrm{E}$ | $w_{0,0}-w_{2,0}$
… | … | …
Figure 1: Scattering phase shifts $\phi^{\bm{d}}_{\Gamma}$ for different total
momenta $\bm{p}=\frac{2\pi}{L}\bm{d}$ and irreducible representations $\Gamma$
calculated based on the relations given in Tab. 4. As input we have used the
pion mass, the rho mass and the volume of the ensemble D200. The dashed
vertical line indicates the rho mass we have measured on this ensemble. The
poles occur at the positions of the noninteracting center of mass energies.
The red and blue curves correspond to the Breit-Wigner ansatz Eq. (49) and to
the Gounaris-Sakurai ansatz Eq. (52), respectively, see also Fig. 2. The
energy levels are at the intersections of the curves. Figure 2: Comparison of
the Breit-Wigner (red) and Gounaris-Sakurai (blue) parametrizations of the
phase shift. Here we use the pion mass and (naively measured) rho mass from
D200 as input. The rho-pi-pi coupling constant is set to the phenomenological
value $g_{\rho\pi\pi}=5.96$.

For particles of integer spin and at zero momentum (i.e., in the center of
mass frame), the full symmetry group on the lattice is the octahedral group
$O_{h}=O\otimes I$ defined as the direct product of the cubic group $O$
(consisting of $24$ rotations) and the group of space inversions $I$.444For
half-integer spin one would have to consider the corresponding double covers
of $O$ and $O_{h}$. In a moving frame, however, the symmetry is reduced to the
so-called little groups (for details see, e.g., refs. [80, 81, 79]) shown in
Tab. 3 together with the decomposition into irreducible representations. The
four-momenta in the laboratory frame are related to those in the center-of-
mass frame by a Lorentz boost, such that the states that obey the quantization
condition on the lattice in the different moving frames will in general
correspond to different center-of-mass energies.

The connection between the finite volume energy spectrum of two-pion states
and infinite volume scattering phase shifts has been established by M. Lüscher
in his seminal articles [82, 83]. Recent discussions of this topic are also
found in refs. [80, 81, 79]. Being interested in the rho resonance in the
vector channel, we may restrict ourselves to the $P$-wave ($l=1$)
contribution, since it is usually found that nonzero phase shifts in higher
odd partial waves are not required to describe the two-pion spectrum [84, 85].
In this simplified situation the $P$-wave phase shift $\delta_{1}$ is directly
related to the quantized two-pion energy levels in finite volume. The latter
appear when the condition

$\displaystyle\cot{\delta_{1}}\overset{!}{=}\cot{\phi^{\bm{d}}_{\Gamma}}$ (45)

is satisfied, see Fig. 1, which will be discussed in more detail below. The
scattering phase shifts $\phi^{\bm{d}}_{\Gamma}$ can be taken from Tab. 4,
using

$\displaystyle w_{lm}$
$\displaystyle=\frac{Z^{\bm{d}}_{lm}(1,q^{2})}{\pi^{3/2}\,\sqrt{2l+1}\,\gamma\,q^{l+1}}\,,$
$\displaystyle q$ $\displaystyle=\frac{Lk}{2\pi}\,,$ (46) $\displaystyle k$
$\displaystyle=\sqrt{\frac{E_{\rm cm}^{2}}{4}-m_{\pi}^{2}}\,.$ (47)

For the numeric evaluation of the generalized zeta function
$Z^{\bm{d}}_{lm}(1,q^{2})$ we use the representation derived in ref. [80].

Using Eq. (45) we obtain the energy levels in the interactive case via
equating the phase shifts given in Tab. 4 with a phenomenological
parametrization, where, for any given parametrization, we define the rho mass
and width as [86]

$\displaystyle\cot{\delta_{1}}\biggr{|}_{s=m_{\rho}^{2}}$ $\displaystyle=0\,,$
$\displaystyle m_{\rho}\Gamma_{\rho}$
$\displaystyle=\biggl{(}\frac{d\delta_{1}}{ds}\biggr{)}^{-1}_{s=m_{\rho}^{2}}\,,$
(48)

using the Mandelstam variable $s=E_{\rm cm}^{2}$. For instance one can use a
relativistic Breit-Wigner (BW) ansatz

$\displaystyle\frac{k^{3}}{\sqrt{s}}\cot{\delta_{1}^{\rm BW}}$
$\displaystyle=\frac{6\pi\bigl{(}m_{\rho}^{2}-s\bigr{)}}{g_{\rho\pi\pi}^{2}}=\frac{k_{\rho}^{3}\bigl{(}m_{\rho}^{2}-s\bigr{)}}{\Gamma_{\rho}m_{\rho}^{2}}\,,$
(49) with $\displaystyle\Gamma_{\rho}$
$\displaystyle=\frac{g_{\rho\pi\pi}^{2}k^{3}_{\rho}}{6\pi m_{\rho}^{2}}\,,$
(50)

where $k=\sqrt{s/4-m_{\pi}^{2}}$, as defined in Eq. (47), and

$\displaystyle k_{\rho}$
$\displaystyle=k\biggr{|}_{s=m_{\rho}^{2}}=\sqrt{m_{\rho}^{2}/4-m_{\pi}^{2}}\,.$
(51)

Alternatively, one can use a Gounaris-Sakurai (GS) parametrization [86], where

$\displaystyle\frac{k^{3}}{\sqrt{s}}\cot{\delta_{1}^{\rm GS}}$
$\displaystyle=k^{2}(h(s)-h(m_{\rho}^{2}))+(k_{\rho}^{2}-k^{2})c\,,$ (52) with
$\displaystyle h(s)$
$\displaystyle=\frac{2}{\pi}\frac{k}{\sqrt{s}}\ln{\biggl{(}\frac{\sqrt{s}+2k}{2m_{\pi}}\biggr{)}}\,,$
(53) $\displaystyle c$
$\displaystyle=\frac{4k_{\rho}^{3}}{m_{\rho}^{2}\Gamma_{\rho}}+4k_{\rho}^{2}h^{\prime}(m_{\rho}^{2})\,.$
(54)

If one chooses to apply the Kawarabayashi-Suzuki-Riazuddin-Fayyazuddin
relation [87, 88], $m_{\rho}^{2}=2g_{\rho\pi\pi}^{2}F_{\pi}^{2}$, (which, as
argued in ref. [89], is a consequence of chiral symmetry and the requirement
of consistency of the effective field theory with respect to
renormalizability), both the BW and the GS parametrizations are determined
solely by the rho mass (given that the pion decay constant $F_{\pi}$ is well-
known and the width of the rho is linked to the rho-pi-pi coupling constant
$g_{\rho\pi\pi}$ via Eq. (50)).

In Fig. 2 we plot the $P$-wave phase shifts for both, the BW (red) and the GS
(blue) parametrization, using the properties of our ensemble D200
(${N_{s}\times N_{t}=64\times
128,m_{\pi}=$201\text{\,}\mathrm{MeV}$,m_{\rho}=$746\text{\,}\mathrm{MeV}$}$)
as input. As one can see, the two parametrizations yield quite similar
results. In Fig. 1, we illustrate the quantization condition Eq. (45): the
energy levels are situated at the intersections between the phase shift
parametrizations and the curves for $\cot{\phi^{\bm{d}}_{\Gamma}}$. The pole
positions correspond to the center of mass energies of the noninteracting
system.

As pointed out in ref. [86], the phase shifts are linked to the pion form
factor via

$\displaystyle F_{\pi}(s)$ $\displaystyle=\frac{f(0)}{f(s)}\,,$ with
$\displaystyle f(s)$
$\displaystyle=\frac{k^{3}}{\sqrt{s}}\bigl{(}\cot{\delta_{1}}-i\bigr{)}\,.$
(55)

Note that this formula only works for the GS ansatz, which we will use in the
following, and not for the BW ansatz, because in the latter case $f(s)$
diverges at $s=0$. Using the GS parametrization one finds

Figure 3: The pion form factor obtained using Eq. (56) and the corresponding
overlap factors $Z^{\pi\pi}_{\Gamma}$ (again using D200 input parameters) for
the local currents from Eq. (57). The red line on the lower panel is plotted
for comparison and corresponds to the estimated
$Z^{\rho}=2f_{\rho}^{2}m_{\rho}^{2}$ for the local currents using the
phenomenological value $f_{\rho}=$222\text{\,}\mathrm{MeV}$$ as input.

$\displaystyle F_{\pi}(s)$
$\displaystyle=\frac{m_{\pi}^{2}h(m_{\rho}^{2})+\frac{c}{4}m_{\rho}^{2}-\frac{1}{\pi}m_{\pi}^{2}}{\frac{k^{3}}{\sqrt{s}}\bigl{(}\cot{\delta_{1}^{\rm
GS}}-i\bigr{)}}\,.$ (56)

As shown in ref. [90] (which is a generalization of the original derivation
given in ref. [91] for moving frames), the form factor can be determined from
the overlap factor of the two-pion states with a local (unsmeared) vector
current. By inverting this relation (and adapting it to our conventions) we
obtain the overlap factors

$\displaystyle Z^{\pi\pi}_{\bm{p}}(s)$
$\displaystyle=\biggl{(}q\frac{\partial\phi^{\bm{d}}_{\Gamma}}{\partial
q}+k\frac{\partial\delta_{1}(k)}{\partial
k}\biggr{)}^{-1}\frac{k^{5}}{6\pi\sqrt{s}}|F_{\pi}(s)|^{2}\,,$ (57)

from a given form factor, which, in turn, can be determined from the phase
shift.

In Fig. 3 we show the form factor and the corresponding estimate for the
overlap factor of local (i.e., unsmeared) vector currents at the source and
the sink with the two-pion states at a given center of mass energy. The first
thing to notice is, that two-pion states whose center of mass energy is much
smaller (or larger) than the rho mass are strongly suppressed. In the example
shown here the overlap of these states is (roughly) smaller by a factor of 100
compared to our estimate for overlap of the rho meson itself (horizontal red
line). This means that these states will not yield large contributions to the
correlation functions at the intermediate time distances available in our
simulation, despite being energetically favored, which may explain why we do
not see these states in our numerical analysis. More problematic is the
possible contribution of states that have a center of mass energy close to the
rho mass. As one can see in Fig. 3, the overlap of these states is strongly
enhanced, and can be of the same size or even larger than the overlap of the
rho meson. This is particularly concerning, because we would not be able to
distinguish such a state in the spectral decomposition within our numerical
analysis. It is important to keep this caveat in mind when interpreting our
results.

That being said, we want to stress that the analysis provided above is
actually only valid for unsmeared currents. Obviously, the situation might be
less critical for the smeared currents that we use in our simulation. A
posteriori, the trustworthiness of the numerical results presented in the
following could be enhanced significantly, if future studies (e.g., by using
the generalized eigenvalue method with two-pion interpolating currents, cf.
refs. [92, 93, 81, 94, 95]) can show that the overlap of smeared vector
interpolating currents with the two-pion states is much smaller than for the
local currents.

##  IV Analysis and results

###  IV.1 Pion and rho mass

To compute the reduced matrix elements introduced, e.g., in Eq. (23), we need
the mass (energy) of the meson in the rest (boosted) frame. While the values
for $m_{\pi}$ are taken from [96] the values for $m_{\rho}$ are obtained by a
direct fit to the correlation function using the spectral decomposition
presented in Eq. (39). Beside the mass (energy) itself two additional
amplitudes ($Z$ and $A$) and also the energy gap to the first excited state
$\Delta E$ enter the fit as free parameters. However, using the ratio method
introduced in Sec. III.3.4 the additional amplitudes and also $\Delta E$ will
not enter the results presented in this work.

The fits are performed using a constant fit window of
${\sim}$2\text{\,}\mathrm{fm}$$ for all ensembles with open boundary
conditions and we start 1 or 2 timeslices
(${\sim}\,$0.1\text{\,}\mathrm{fm}$$) away from the source for the coarser or
finer lattices. Due to the structure of Eq. (39) the values obtained by the
fit for the ground state and excited state energy can be interchanged. To
overcome this technical issue we have introduced a cutoff for the double
exponential fit at $t_{\mathrm{cut}}\approx$0.65\text{\,}\mathrm{fm}$$ and fit
only the single exponential
$Z^{\rho}_{\bm{p}}(2E^{\rho}_{\bm{p}})^{-1}e^{-E^{\rho}_{\bm{p}}t}$ for larger
times. In case of periodic boundary conditions we choose a symmetric ansatz of
the form
$Z^{\rho}_{\bm{p}}(2E^{\rho}_{\bm{p}})^{-1}(e^{-E^{\rho}_{\bm{p}}t}+e^{-E^{\rho}_{\bm{p}}(T-t)})$
with $t_{\mathrm{cut}}\leq t\leq T-t_{\mathrm{cut}}$ where we only fit the
amplitude and the ground state energy. The final results of these fits are
shown in Fig. 4. We depict the rest frame results by triangles and the boosted
frame results by circles. Note that the continuum dispertion relation
$E_{\bm{p}}^{\rho}=\sqrt{m_{\rho}^{2}+\bm{p}^{2}}$ is used to project the
energies onto their corresponding mass values and to check that the dispersion
relation is well satisfied for the momenta in use.

Figure 4: Rho masses for all ensembles analyzed from a double exponential fit
(open boundaries) or a single exponential fit (periodic boundaries) to the
two-point function correlator. The triangles depict the results in the rest
frame ($\bm{p}^{2}=0$) while the circles correspond to fits in the boosted
frame ($\bm{p}^{2}=\frac{4\pi^{2}}{L^{2}}$) projected to the rest frame using
the continuum dispersion relation.

###  IV.2 Extraction of ground state matrix elements

The observables studied in this article are affected by disconnected quark
loops. Often one can circumvent this problem by considering isovector current
insertions, where the up and down quark disconnected loops cancel each other
identically in the limit of exact isospin symmetry. This is not a viable
solution in this case, since also the connected part vanishes for the
isovector currents. Unfortunately, the disconnected contributions are
notorious for having a large statistical error. However, as will be discussed
later in this section, this is not true in general.

Figure 5: Extraction of renormalized values for $v_{2}$, $a_{2}$, and $d_{2}$
from the ratios obtained in Eqs. (58)–(60) using the operator combination
$\mathcal{O}_{\rm v2a}$. For illustrative purposes we only show the data
points and individual fits (solid blue lines in the first column) for all
momentum combinations with $\bm{n}^{2}=1$, however, the ground state results
(orange lines) are obtained by a simultaneous fit to the operator combinations
$\mathcal{O}_{\rm v2a}$ and $\mathcal{O}_{\rm v2b}$ using all possible momenta
for the corresponding matrix element with $\bm{n}^{2}\leq 1$ for the
connected, disconnected non-singlet, and disconnected singlet contributions
respectively. The analysis shown in this plot has been performed on the
ensemble N204. The solid blue lines in the first column correspond to a
simultaneous fit to the four source-sink separations of the ensemble, cf. Tab.
1, for the insertion current $(\mathrm{bar}{u}u+\mathrm{bar}{d}d)$ needed to
construct the flavor (non-)singlet operator contributions. In case of $v_{2}$
and $a_{2}$, the fits allow for a generic excited state on top of the ground
state, while in case of $d_{2}$ the excited state energy $\Delta E$ is fixed
by the two-point function, see Eq. (39). The orange line depicts the extracted
ground state contribution and directly correspond to the values of the reduced
matrix elements. In the second and third column we show the disconnected
contributions for the non-singlet
($\mathrm{bar}{u}u+\mathrm{bar}{d}d-2\mathrm{bar}{s}s$) and singlet
($\mathrm{bar}{u}u+\mathrm{bar}{d}d+\mathrm{bar}{s}s$) operators as function
of the final timeslice $t$. In addition to the original data points (grayed
out) we also show an average over all insertion times for every final
timeslice $t$ (black triangle markers). Figure 6: Extraction of renormalized
values for $v_{2}$, $a_{2}$, and $d_{2}$ from the ratios obtained in Eqs.
(61)–(63) using the operator combination $\mathcal{O}_{\rm v2b}$. The data is
visualized in the same way as in Fig. 5, i.e., we only plot the data points
and individual fits for $\bm{n}^{2}=1$. Also here we want to stress, that the
ground state results (orange line) are obtained by a simultaneous fit to the
operator combinations $\mathcal{O}_{\rm v2a}$ and $\mathcal{O}_{\rm v2b}$
using all possible momenta for the corresponding matrix element with
$\bm{n}^{2}\leq 1$ for the connected, disconnected non-singlet, and
disconnected singlet contributions respectively.

In cases where the disconnected contributions are zero within the error one
might be tempted to simply drop them. However, in situations where the
statistical error is large (for instance the flavor singlet operators), their
inclusion can shift the mean and, even more important, can increase the error
for the final result substantially. I.e., they have to be included, if one
wants to provide reliable error estimates for phenomenological applications.
Nevertheless, we perform a second analysis in these cases, where we solely use
the connected part, which allows us to compare to other lattice results for
connected contributions.

In principle, one would like to add the connected and the disconnected
contributions already at the correlation function level. However, this is not
feasible since the connected and the disconnected parts are calculated in
different ways. The connected part is calculated with the stochastic
propagator estimation presented in App. C.1 using a setup with a fixed sink
timeslice to obtain the result for all possible insertion times. The
disconnected loops are calculated on fixed timeslices also using stochastic
propagators (see App. C.2 for details). For any computed two-point function we
can therefore obtain the disconnected three-point function at a fixed
insertion time but for arbitrary final times. In order to sum up the
contributions at the correlation function level, we would have to throw away a
large part of our data, since we could only use the insertion and final times
where we have data for both. This would be prohibitively wasteful. Therefore,
we will perform the extraction of the ground state matrix elements separately
for the connected and the disconnected contribution.

Next, we extract the ground state contribution from the ratios defined in Sec.
III.3.4. We use kinematic prefactors and, in case of the rho, take appropriate
linear combinations such that the ground state contributions directly
correspond to the reduced matrix elements $v_{2}^{q}$ (for the pion) as well
as $a_{2}^{q}$ and $d_{2}^{q}$ for the rho meson. For on-axis momenta
$\hat{\bm{p}}=\pm\bm{e}_{i}$ we obtain (using Eqs. (83) and (85))

$\displaystyle v_{2}$
$\displaystyle=\frac{1}{p^{i}}R_{\bm{p}}\bigl{(}\mathcal{O}_{\rm
v2a}^{i}\bigr{)}\,,$ (58) $\displaystyle a_{2}$
$\displaystyle=\frac{1}{3p^{i}}\sum_{j}R^{j}_{\bm{p}}\bigl{(}\mathcal{O}_{\rm
v2a}^{i}\bigr{)}\,,$ (59) $\displaystyle d_{2}$
$\displaystyle=\frac{3}{4p^{i}}\bigl{(}2R^{i}_{\bm{p}}\bigl{(}\mathcal{O}_{\rm
v2a}^{i}\bigr{)}-\sum_{j\neq i}R^{j}_{\bm{p}}\bigl{(}\mathcal{O}_{\rm
v2a}^{i}\bigr{)}\bigr{)}\,,$ (60)

with the operator $\mathcal{O}_{\rm v2a}^{i}$ as inserted current. For
$\mathcal{O}_{\rm v2b}$ we find

$\displaystyle v_{2}$
$\displaystyle=\frac{3E}{4E^{2}-m^{2}}R_{\bm{p}}\bigl{(}\mathcal{O}_{\rm
v2b}\bigr{)}\,,$ (61) $\displaystyle a_{2}$
$\displaystyle=\frac{E}{4E^{2}-m^{2}}\sum_{j}R^{j}_{\bm{p}}\bigl{(}\mathcal{O}_{\rm
v2b}\bigr{)}\,,$ (62) $\displaystyle d_{2}$
$\displaystyle=\frac{3E}{8(E^{2}-m^{2})}\bigl{(}2R^{i}_{\bm{p}}\bigl{(}\mathcal{O}_{\rm
v2b}\bigr{)}-\sum_{j\neq i}R^{j}_{\bm{p}}\bigl{(}\mathcal{O}_{\rm
v2b}\bigr{)}\bigr{)}\,.$ (63)

Here we hide the superscripts $\pi$, $\rho$, and the subscript $\bm{p}$ for
the mass and energy. Measurements using $\mathcal{O}_{\rm v2a}^{i}$ always
require nonzero momentum in direction $i$. If one uses $\mathcal{O}_{\rm v2b}$
the reduced matrix elements $v_{2}$ and $a_{2}$ can be measured for vanishing
three-momentum. The extraction of $d_{2}$, however, always requires nonzero
momentum. The reason for this is that $d_{2}$ corresponds to the difference of
the PDF moment between longitudinally and transversally polarized rho mesons,
which is no useful concept for mesons in their rest frame. Explicit operator
definitions can be found in App. B.

In Figs. 5 and 6, we show examples of the ratio fits for the ensemble N204.
For the statistical analysis we generate 500 bootstrap samples per ensemble
using a bin size of 40 molecular dynamics units to eliminate autocorrelations.
To visualize the fits and corresponding data points we only present plots for
$\bm{n}^{2}=1$, where $\bm{n}$ is defined by $\bm{p}=\frac{2\pi}{L}\,\bm{n}$,
using the operator combinations Eqs. (58)–(63). However, the results for the
reduced matrix elements $v_{2}$, $a_{2}$, and $d_{2}$ in the left column are
obtained by combined fits to all ratios using an ansatz similar to Eq. (40)
with the definitions Eqs. (41) and (42). In case of pseudoscalar mesons this
reads

$\displaystyle R(\mathcal{O},\bm{p}^{2},t,\tau)=B_{0}$
$\displaystyle+B_{1}\left(\mathcal{O},\bm{p}^{2}\right)\,\,e^{-\Delta
E_{\bm{p}^{2}}\,\,\left(t-\tau\right)}$
$\displaystyle+B_{1}\left(\mathcal{O},\bm{p}^{2}\right)\,\,e^{-\Delta
E_{\bm{p}^{2}}\,\,\tau},$ (64)

where the ratio $R$ explicitly depends on $\bm{p}^{2}$, the operator
$\mathcal{O}\in\left\\{\mathcal{O}_{\rm v2a},\mathcal{O}_{\rm v2b}\right\\}$,
the sink timeslice $t$, and the insertion timeslice $\tau$. For the operator
combination $\mathcal{O}_{\rm v2a}^{i}$, with $i=1,2,3$, we average the three
spatial directions and fit to the data using the parameters $B_{0}$,
$B_{1}\left(\mathcal{O}_{\rm v2a},\bm{n}^{2}=1\right)$ and the excited state
energy is given by $\Delta E_{\bm{n}^{2}=1}$, independent of the operator
combination. Note that we require nonzero momentum in direction $i$ for
$\mathcal{O}_{\rm v2a}$. Additionally the operator combination
$\mathcal{O}_{\rm v2b}$ gives rise to the further excited state amplitudes
$B_{1}\left(\mathcal{O}_{\rm v2b},\bm{n}^{2}\right)$ and the excited state
energy $\Delta E_{\bm{n}^{2}=0}$. All in all this yields a simultaneous fit to
three operator combinations for each source sink separation, to resolve the
individual parameters. The actual fit is performed simultaneously to all
source sink separations, cf. Tab. 1. A summary of the individual contributions
is given in Tab. 5. One can easily deduce that the ground state contribution
$B_{0}$ is present in all operator combinations while the excited state
amplitudes and energies depend on the operator combination and $\bm{n}^{2}$
respectively. A similar approach holds for the extraction of $a_{2}$ and
$d_{2}$. However, for $d_{2}$ the statistical error is much larger and we are
not able to resolve reliable excited state energies from the ratios. Therefore
we fix the excited state energies $\Delta E_{\bm{n}^{2}=1}$ by an additional,
simultaneous fit to the two-point function Eq. (39) in case of the reduced
matrix element $d_{2}$.

Table 5: Summary of the occurrence of the individual fit parameters in the ansatz Eq. (IV.2) for the extraction of $v_{2}$ (pion) and $a_{2}$ (rho). A green check mark indicates that the fit parameter is present in the corresponding operator combination, whereas the red crosses indicate that the fit parameter is not present. Fit Parameter | $\mathcal{O}_{\rm v2a}(\bm{n}^{2}=1)$ | $\mathcal{O}_{\rm v2b}(\bm{n}^{2}=1)$ | $\mathcal{O}_{\rm v2b}(\bm{n}^{2}=0)$
---|---|---|---
$B_{0}$ | ✓ | ✓ | ✓
$B_{1}\left(\mathcal{O}_{\rm v2a},\bm{n}^{2}=1\right)$ | ✓ | ✗ | ✗
$B_{1}\left(\mathcal{O}_{\rm v2b},\bm{n}^{2}=1\right)$ | ✗ | ✓ | ✗
$B_{1}\left(\mathcal{O}_{\rm v2b},\bm{n}^{2}=0\right)$ | ✗ | ✗ | ✓
$\Delta E_{\bm{n}^{2}=\,1}$ | ✓ | ✓ | ✗
$\Delta E_{\bm{n}^{2}=\,0}$ | ✗ | ✗ | ✓

In summary this implies that the fits do not only take into account the data
points shown in Figs. 5 and 6, but are actually based on a larger data set
stemming from various operator and momentum configurations. However, the
prefactors in Eqs. (58) – (63) are solely determined by the meson masses and
corresponding momentum contributions.

As discussed above, we perform the extraction of the ground state matrix
elements separately for the connected (left column) and the disconnected
(middle and right column) contributions. While analyzing the disconnected
contribution we found that the considerable noise on the light and strange
quark loops is highly correlated for all included ensembles, cf. Tab. 1. We
can use this to our advantage by looking at the non-singlet
($\mathrm{bar}{u}u+\mathrm{bar}{d}d-2\mathrm{bar}{s}s$) and singlet
($\mathrm{bar}{u}u+\mathrm{bar}{d}d+\mathrm{bar}{s}s$) flavor combinations
instead of the light and strange loops themselves. As depicted impressively in
the middle and the right column of Figs. 5 and 6 (note the difference in
scale), the statistical error is smaller by more than one order of magnitude
for the non-singlet operator. For the disconnected contributions we do not see
an indication for a significant excited state contribution and, consequently,
content ourselves with a constant fit to extract the ground state signal
employing the fit strategy described above.

For the disconnected contribution we have data points for a large number of
combinations of final times $t$ and insertion times $\tau$. If one plots the
data for various insertions times in one plot (cf. the grayed out points in
Figs. 5 and 6), the statistical scattering of these data points alone can lead
to the misconception that the statistical error of the extracted ground state
(yellow band) is underestimated. In order to convince the viewer of the plots
that this is not the case we also plot the black points, which are obtained by
taking the average over data at all insertion times $\tau$ for the given final
time $t$.

###  IV.3 Quark mass dependence and continuum extrapolation

As described in Sec. III.1, the ensembles we analyze have been generated along
multiple trajectories in the quark mass plane and at different lattice
spacings $a$. We obtain the final results by extrapolating to the continuum
limit (at $a=0$) and to physical masses. To this end we employ the
parametrization

$\displaystyle f(a,m_{\pi}^{2},m_{K}^{2})$
$\displaystyle=c_{0}+c_{1}a+c_{2}m_{\pi}^{2}+c_{3}m_{K}^{2}\,.$ (65)

Note that we have to use a linear term in the lattice spacing as leading
contribution despite the fact that our lattice action is order $a$ improved,
because we lack the order $a$ improvement for the inserted currents.

In Figs. 7 and 8, the yellow bands show the extrapolations for the flavor non-
singlet and flavor singlet operators as a function of $a$ and $m_{\pi}^{2}$,
respectively. As discussed in Sec. IV.2, the statistical error of the flavor
singlet operator combinations is much larger than in the flavor non-singlet
case, which makes it hard to draw any convincing conclusions. Nevertheless,
the results for the flavor singlet combinations will allow us to give at least
an upper bound for the reduced matrix elements $v_{2}$, $a_{2}$, and $d_{2}$.
However, for the flavor non-singlet combinations, the situation is much better
and we get meaningful results and errors. Fig. 8 shows the extrapolations for
the quark mass dependence along the ${\rm TrM}=\text{const.}$,
$m_{s}=\text{const.}$, and $m_{l}=m_{s}$ trajectories, from left to right. The
data points are corrected for lattice spacing effects and are shifted to the
corresponding trajectories. For, e.g., the ${\rm TrM}=\text{const.}$
trajectory this keeps the average quark mass fixed (again using the fitted
model shown above) and thus allows us to see the effect of flavor symmetry
breaking. The lattice spacing dependence is depicted in Fig. 7. To visualize
solely the discretization effects, the data points in Fig. 7 are corrected for
mass effects (using the fitted model), i.e., they are translated to physical
masses along the fitted curve, and finally averaged for all ensembles with the
same values of $\beta$ using the weighted average

$\displaystyle\mathrm{bar}{A}_{i}$
$\displaystyle=\sum_{i}^{N_{\beta}}w_{i}\,A_{i},$ with $\displaystyle w_{i}$
$\displaystyle=\frac{1/\sigma_{i}^{2}}{\sum_{i}^{N_{\beta}}1/\sigma_{i}^{2}},$
(66)

where $N_{\beta}$ is the number of data points $A_{i}$ (ground state matrix
elements) per $\beta$ with corresponding errors $\sigma_{i}$. Note that this
procedure is only applied to the points in these plots for illustrative
purposes, while the bands are obtained from the actual fit performed using the
original data points, cf. App. D.

Figure 7:  Lattice spacing dependence of the extrapolations for the flavor
singlet ($\mathrm{bar}{u}u+\mathrm{bar}{d}d+\mathrm{bar}{s}s$) and flavor non-
singlet ($\mathrm{bar}{u}u+\mathrm{bar}{d}d-2\mathrm{bar}{s}s$) operator
combinations using the fits shown in, e.g., Figs. 5 and 6. The reduced matrix
elements $v_{2}$ (pion), $a_{2}$, and $d_{2}$ (both rho) have been obtained by
a translation to physical quark masses and averaging measurements with the
same values of $\beta$ using the weighted average given in Eq. (66). From
coarsest to finest lattice spacing, this corresponds to averaging the data of
2, 8, 7, 7, and 3 independent ensembles, cf. Tab. 1. The non-averaged plots
for the individual trajectories can be found in appendix D. Figure 8:
Extrapolation for the flavor singlet
($\mathrm{bar}{u}u+\mathrm{bar}{d}d+\mathrm{bar}{s}s$) and flavor non-singlet
operator combinations ($\mathrm{bar}{u}u+\mathrm{bar}{d}d-2\mathrm{bar}{s}s$)
using the fits shown in, e.g., Fig. 5. From left to right we show the
$m_{\pi}^{2}$ dependence of the reduced matrix elements $v_{2}$ (pion),
$a_{2}$, and $d_{2}$ (both rho) for the three different trajectories we use in
our analysis. Note that the symmetric trajectory $m_{l}=m_{s}$ with exact
$\rm{SU}(3)$ flavor symmetry approaches the chiral limit in the quark mass
plane and not the physical point.

The final results for the reduced matrix elements are given in Tab. 6. In
addition to the statistical error $()_{s}$ we provide estimates for the
systematic uncertainties due to the quark mass extrapolation $()_{m}$ and the
continuum extrapolation $()_{a}$. To this end, we have performed additional
fits with cuts in the mass range
$\mathrm{bar}{m}=\sqrt{(2m_{K}^{2}+m_{\pi}^{2})/3}<$450\text{\,}\mathrm{MeV}$$
and $a<$0.09\text{\,}\mathrm{fm}$$, respectively. We then take the difference
between the results from these fits and our main result as an estimate of the
corresponding systematic uncertainties.

###  IV.4 Discussion

Using Eq. (14) in combination with Eqs. (25) and (26), one can write the ratio
of the moments of the structure functions (at leading twist accuracy) with the
corresponding Wilson coefficient $C^{(k)}_{n}=1+\mathcal{O}(\alpha_{s})$ as a
sum over the related reduced matrix element

$\displaystyle\frac{2M_{2}(F_{1})_{\pi}}{C_{2}^{(1)}}=\frac{M_{1}(F_{2})_{\pi}}{C_{2}^{(2)}}$
$\displaystyle=\sum_{q}e_{q}^{2}v_{2}^{q}=\frac{2}{9}\Bigl{(}v_{2}^{\rm
fs}+\frac{1}{4}v_{2}^{\rm fns}\Bigr{)}\,,$ (67)
$\displaystyle\frac{2M_{2}(F_{1})_{\rho}}{C_{2}^{(1)}}=\frac{M_{1}(F_{2})_{\rho}}{C_{2}^{(2)}}$
$\displaystyle=\sum_{q}e_{q}^{2}a_{2}^{q}=\frac{2}{9}\Bigl{(}a_{2}^{\rm
fs}+\frac{1}{4}a_{2}^{\rm fns}\Bigr{)}\,,$ (68)
$\displaystyle\frac{2M_{2}(b_{1})_{\rho}}{C_{2}^{(1)}}=\frac{M_{1}(b_{2})_{\rho}}{C_{2}^{(2)}}$
$\displaystyle=\sum_{q}e_{q}^{2}d_{2}^{q}=\frac{2}{9}\Bigl{(}d_{2}^{\rm
fs}+\frac{1}{4}d_{2}^{\rm fns}\Bigr{)}\,,$ (69)

where we assume exact isospin symmetry and ${\rm fs}\equiv u+d+s$ is the
flavor singlet while ${\rm fns}\equiv u+d-2s$ is the flavor non-singlet
combination. If we only consider the connected part, the strange quark
contribution drops out entirely. The result can be written in terms of the
light quark connected contribution as

$\displaystyle\frac{2M_{2}(F_{1})_{\pi}^{\rm
conn.}}{C_{2}^{(1)}}=\frac{M_{1}(F_{2})_{\pi}^{\rm conn.}}{C_{2}^{(2)}}$
$\displaystyle=\frac{5}{9}v_{2}^{\ell,{\rm conn.}}\,,$ (70)
$\displaystyle\frac{2M_{2}(F_{1})_{\rho}^{\rm
conn.}}{C_{2}^{(1)}}=\frac{M_{1}(F_{2})_{\rho}^{\rm conn.}}{C_{2}^{(2)}}$
$\displaystyle=\frac{5}{9}a_{2}^{\ell,{\rm conn.}}\,,$ (71)
$\displaystyle\frac{2M_{2}(b_{1})_{\rho}^{\rm
conn.}}{C_{2}^{(1)}}=\frac{M_{1}(b_{2})_{\rho}^{\rm conn.}}{C_{2}^{(2)}}$
$\displaystyle=\frac{5}{9}d_{2}^{\ell,{\rm conn.}}\,.$ (72)

Table 6: Results obtained from the extrapolations in Figs. 7 and 8 and the coresponding connected-only contributions for the flavor combination $(u+d)$, all at $\mu=$2\text{\,}\mathrm{GeV}$$. The final statistical error is given by $()_{s}$ and estimates of the systematic uncertainties due to the quark mass by $()_{m}$, and due to the continuum extrapolation by $()_{a}$. The values of $\chi^{2}$ per degrees of freedom are obtained from the corresponding extrapolations. Matrix element | Final result | $\chi^{2}/\rm{d.o.f.}$
---|---|---
$v_{2}^{(u+d+s)}$ | 0.220 $(95)_{s}(98)_{m}(155)_{a}$ | 1.71
$a_{2}^{(u+d+s)}$ | 0.285 $(90)_{s}(76)_{m}(271)_{a}$ | 1.78
$d_{2}^{(u+d+s)}$ | 0.226 $(112)_{s}(6)_{m}(54)_{a}$ | 0.72
$v_{2}^{(u+d-2s)}$ | 0.344 $(20)_{s}(3)_{m}(19)_{a}$ | 1.72
$a_{2}^{(u+d-2s)}$ | 0.384 $(29)_{s}(11)_{m}(42)_{a}$ | 1.33
$d_{2}^{(u+d-2s)}$ | 0.163 $(38)_{s}(5)_{m}(7)_{a}$ | 0.58
$v_{2}^{(u+d),\rm{conn.}}$ | 0.357 $(16)_{s}(2)_{m}(15)_{a}$ | 1.74
$a_{2}^{(u+d),\rm{conn.}}$ | 0.393 $(29)_{s}(10)_{m}(35)_{a}$ | 1.47
$d_{2}^{(u+d),\rm{conn.}}$ | 0.180 $(38)_{s}(5)_{m}(7)_{a}$ | 0.59

In Tab. 7 we give our final results for these linear combinations. As shown in
the last section the flavor singlet contributions contain relatively large
errors which affect Eqs. (67)–(69). However, treating the connected part
only555The connected only results are presented as a comparison option for
other studies neglecting disconnected contributions. reduces the errors
significantly. The reader should be aware of the fact that we use (in both
cases) the flavor non-singlet renormalization constants, which is only an
approximation, cf. the discussion in Sec. III.2.

For the flavor non-singlet contributions given in Tab. 6 we obtain very
precise results, despite the fact that all disconnected quark loops are fully
taken into account. The non-singlet operators do not mix with gluonic
operators under renormalization, and the necessary renormalization factors
have been calculated non-perturbatively, cf. Sec. III.2. As a first step one
can compare the connected results in Tab. 7 to the connected-only and the
flavor non-singlet results in Tab. 6. Multiplying the latter two by the
prefactors given in Eqs. (70)–(72) one finds that both the first moments of
the unpolarized structure functions $F_{1}^{\pi}$ and $F_{1}^{\rho}$ and the
first moment of the structure function $b_{1}$ are in very good agreement with
the flavor non-singlet results in Tab. 6 and of course also reflect the
connected only result in Tab. 6. As shown in Sec. IV.2 the flavor non-singlet
contributions of the quark line disconnected diagrams in the extraction of the
ground state matrix elements are small compared to the connected
contributions. In leading order the structure function $F^{q}_{1}(x)$
corresponds to one half of the probability to find a quark of flavor $q$ with
momentum fraction $x$. If we assume exact $\rm{SU}(3)$ flavor symmetry for the
quark sea, the results in Tab. 6 and 7 imply that in the pion the valence
quarks carry about 35% of the total momentum, while in the rho they carry
about 40% of the total momentum. It is remarkable that these values justify
the assumption $F_{1}(x)^{\pi}\sim F_{1}(x)^{\rho}$, which is often used in
phenomenological estimates. The structure functions $b_{1}(x)$ and $b_{2}(x)$
are sensitive to a possible dependence of the quark densities on the hadron
polarization, i.e., they measure the difference in quark distributions of a
spin projected $\lambda=0$ and $\lambda=+/-$ rho meson. If the quarks were in
a relative S-wave state, cf. the discussion in Sec. III.4, one would expect
$b_{1}=b_{2}=0$. However, our results show a large contribution (compared to
the scale of $a_{2}$) to the approximated valence quark contribution $d_{2}$
with a relative error of only $\sim$10%. This confirms the conclusion in [18]
that the quarks carry substantial angular momentum and also reflects the
results of the various phenomenological studies cited in Sec. I.

##  V Summary and outlook

In this article we have presented the computation of the first moments for the
structure functions $F_{1}^{\pi}$, $F_{1}^{\rho}$, and $b_{1}$ including quark
line disconnected contributions. Despite the fact that our final results are
tainted with large statistical errors due to the flavor singlet disconnected
contributions we are able to provide very accurate results for the flavor non-
singlet combination $u+d-2s$. As an additional subtlety we had a closer look
at possible two-pion contributions which might occur in our analysis. We do
not find any evidence for the contribution of two-pion states. However, as
discussed at the end of Sec. III.4, using our analysis technique we cannot
fully exclude them either. This is particularly true for two-pion states close
to the resonance energy. To this end, the trustworthiness of our numerical
results could be enhanced a posteriori, if future studies (e.g., by using the
generalized eigenvalue method with two-pion interpolating currents, cf. refs.
[92, 93, 81, 94, 95]) can show that the overlap of smeared vector
interpolating currents with the two-pion states is much smaller than for the
local currents.

Table 7: Estimated results for the first moments of the structure functions $F_{1}$ and $b_{1}$ exploiting Eqs. (67)–(72) and performing the extrapolations as discussed in Sec. IV.3. Further we use the abbreviation $C\equiv 2/C_{2}^{(1)}$. Structure function | Final result | $\chi^{2}/\rm{d.o.f.}$
---|---|---
$CM_{2}(F_{1})_{\pi}$ | 0.132 $(33)_{s}(32)_{m}(57)_{a}$ | 1.75
$CM_{2}(F_{1})_{\rho}$ | 0.156 $(33)_{s}(23)_{m}(102)_{a}$ | 1.84
$CM_{2}(b_{1})_{\rho}$ | 0.108 $(41)_{s}(1)_{m}(13)_{a}$ | 0.72
$CM_{2}(F_{1})_{\pi}^{\text{conn.}}$ | 0.099 $(5)_{s}(0)_{m}(4)_{a}$ | 1.74
$CM_{2}(F_{1})_{\rho}^{\text{conn.}}$ | 0.109 $(8)_{s}(2)_{m}(9)_{a}$ | 1.47
$CM_{2}(b_{1})_{\rho}^{\text{conn.}}$ | 0.050 $(10)_{s}(1)_{m}(2)_{a}$ | 0.59

Despite the fact that we for the first time present comprehensive results
including disconnected contributions we have reduced the statistical error
considerably. This can be seen comparing the error of the connected
contribution alone with earlier studies. However, to determine the
phenomenologically important moments of the structure functions (at leading
twist), one needs the flavor singlet combination, where the statistical error
is still large. Future studies will have to aim at a further reduction of
these statistical errors. Once this is achieved, also a nonperturbative
calculation of the singlet renormalization factors and the inclusion of mixing
with gluonic operators might be worthwhile.

###### Acknowledgements.

Support of this project was granted by the German DFG (SFB/TRR 55). In
addition this project has received funding from the European Union’s Horizon
2020 research and innovation program under the Marie Skłodowska-Curie grant
agreement No. 813942. We are grateful to Gunnar S. Bali, Lorenzo Barca, Sara
Collins, Vladimir Braun, Meinulf Göckeler, and Christoph Lehner for the
various fruitful discussions and to Wolfgang Söldner for providing parts of
the intermediate results to be published in [96]. In addition we thank
Benjamin Gläßle and Simon Heybrock for co-developing some of the software used
here and Daniel Richtmann, Peter Georg and Jakob Simeth for software
development and software support and Enno E. Scholz for the support in data
management. We gratefully acknowledge computing time granted by the John von
Neumann Institute for Computing (NIC), provided on the Booster partition of
the supercomputer JURECA [97] and use of computing time and services on the
HDF Cloud [98], funded as part of the Helmholtz Data Federation (HDF)
strategic initiative, at Jülich Supercomputing Centre (JSC). Additional
simulations were carried out at the QPACE2 and QPACE 3 Xeon Phi cluster of
SFB/TRR 55 und the Regensburg computing cluster QPACE B and the Regensburg HPC
cluster Athene2. We owe special thanks to Randy Rückner for the software and
runtime support concerning the Athene2 compute cluster. The authors also
gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for providing
computing time for GCS Large-Scale Projects on SuperMUC and SuperMUC NG at
Leibniz Supercomputing Centre (LRZ). GCS is the alliance of the three national
supercomputing centres HLRS (Universität Stuttgart), JSC (Forschungszentrum
Jülich), and LRZ (BayerischeAkademie der Wissenschaften), funded by the German
Federal Ministry of Education and Research (BMBF) and the German State
Ministries for Research of Baden-Württemberg(MWK), Bayern (StMWFK) and
Nordrhein-Westfalen (MIWF). We thank all our CLS colleagues for the joint
generation of the gauge ensembles. The ensembles were generated using the
OpenQCD [60] software package.

##  Appendix A Light-cone coordinates and polarization vectors

We define the light-cone coordinates used in Sec. II in such a way that the
perpendicular (or transverse) part of the momentum always vanishes, i.e.,
$p_{T}=0$ and $p_{+}$ is its large component. To achieve this let $v$ be any
four-vector and $\hat{\bm{p}}$ the direction of the three-momentum. Then,

$\displaystyle v^{\mu}$
$\displaystyle=v^{+}n_{+}^{\mu}+v^{-}n_{-}^{\mu}+v_{T}^{\mu}\,,$
$\displaystyle\text{ with }n_{\pm}^{\mu}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\\
\pm\hat{\bm{p}}\end{pmatrix}^{\mu},$ (73)

where $v^{\pm}=n_{\mp}\cdot v$ and

$\displaystyle v_{T}^{\mu}$ $\displaystyle=\begin{pmatrix}0\\\
\bm{v_{T}}\end{pmatrix}^{\mu}\,,$ $\displaystyle\text{ with }\bm{v}_{T}$
$\displaystyle=\bm{v}-(\hat{\bm{p}}\cdot\bm{v})\hat{\bm{p}}\,,$ (74)

such that $\bm{v}_{T}\perp\hat{\bm{p}}$. For the momentum we then have
$p^{\pm}=E\pm|\bm{p}|$ and $p_{T}=0$.

For the polarization vectors we use a dimensionless definition. They obey the
general transversality condition

$\displaystyle\sum_{\lambda}\,\epsilon_{\mu}(\bm{p},\lambda)\epsilon^{*}_{\nu}(\bm{p},\lambda)=-\left(g_{\mu\nu}-\frac{p_{\mu}p_{\nu}}{m^{2}}\right)\,,$
(75)

where $m$ is the hadron mass. For momenta in $x$ direction the polarization
vectors in the rest frame are given by $(0,\bm{e}_{\lambda})$ with

$\displaystyle\bm{e}_{0}^{x}=\bm{e}_{0}\biggr{|}_{\hat{\bm{p}}=\hat{\bm{e}}_{1}}$
$\displaystyle=\begin{pmatrix}1\\\ 0\\\ 0\end{pmatrix}\,,$
$\displaystyle\bm{e}_{\pm}^{x}=\bm{e}_{\pm}\biggr{|}_{\hat{\bm{p}}=\hat{\bm{e}}_{1}}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{pmatrix}0\\\ \mp 1\\\
-i\end{pmatrix}\,,$ (76)

where $\bm{e}_{0}^{x}$ corresponds to the longitudinal polarization, while
$\bm{e}_{\pm}^{x}$ to the circular polarizations. For momenta in arbitrary
direction, we have to rotate these vectors to

$\displaystyle\bm{e}_{0}$ $\displaystyle=\hat{\bm{p}}\,,$
$\displaystyle\bm{e}_{\pm}$
$\displaystyle=\bm{e}_{\pm}^{x}-\frac{\hat{\bm{p}}\cdot\bm{e}_{\pm}^{x}}{1+\hat{\bm{p}}\cdot\bm{e}_{0}^{x}}(\bm{e}_{0}^{x}+\hat{\bm{p}})\,,$
(77)

to obtain the longitudinal polarization vector $\bm{e}_{0}$ and the
polarization vectors for the circular polarizations $\bm{e}_{\pm}$. Last but
not least, we have to perform a boost to the laboratory frame. This only
affects $\bm{e}_{0}$ ($\bm{e}_{\pm}$ are invariant because they are
perpendicular to $\bm{p}$), and we obtain

$\displaystyle\epsilon^{\mu}(\bm{p},0)$
$\displaystyle=\begin{pmatrix}\frac{|\bm{p}|}{m}\\\
\frac{E}{m}\hat{\bm{p}}\end{pmatrix}^{\mu}\,,$
$\displaystyle\epsilon^{\mu}(\bm{p},\pm)$ $\displaystyle=\begin{pmatrix}0\\\
\bm{e}_{\pm}\end{pmatrix}^{\mu}\,.$ (78)

In terms of the light-cone coordinates introduced above this yields

$\displaystyle\epsilon^{\pm}(\bm{p},0)$ $\displaystyle=\pm p^{\pm}/m\,,$
$\displaystyle\epsilon^{\mu}_{T}(\bm{p},0)$ $\displaystyle=0\,,$ (79)
$\displaystyle\epsilon^{\pm}(\bm{p},+)$ $\displaystyle=0\,,$
$\displaystyle\epsilon^{\mu}_{T}(\bm{p},+)$
$\displaystyle=\epsilon^{\mu}(\bm{p},+)\,,$ (80)
$\displaystyle\epsilon^{\pm}(\bm{p},-)$ $\displaystyle=0\,,$
$\displaystyle\epsilon^{\mu}_{T}(\bm{p},-)$
$\displaystyle=\epsilon^{\mu}(\bm{p},-)\,.$ (81)

##  Appendix B Operators and matrix element decomposition

To avoid mixing as far as possible we use operators from suitably chosen
multiplets that possess a definite C-parity and transform according to
irreducible representations of $\mathrm{H}(4)$, cf. refs. [77, 18]. To be
specific, we will use the operators $\mathcal{O}_{\rm
v2a}^{i}=\mathcal{O}^{0i}$ and $\mathcal{O}_{\rm
v2b}=\frac{4}{3}\mathcal{O}^{00}$. For the special case of two indices the
action of the symmetrizing and trace-subtracting operator $\mathcal{S}$ is
defined by

$\displaystyle\mathcal{S}\mathcal{O}^{\mu\nu}$
$\displaystyle=\mathcal{S}^{\mu\nu}_{\rho\sigma}\mathcal{O}^{\rho\sigma}=\frac{1}{2}\Bigl{(}g^{\mu}_{\rho}g^{\nu}_{\sigma}+g^{\mu}_{\sigma}g^{\nu}_{\rho}-\frac{2}{d}g^{\mu\nu}g_{\rho\sigma}\Bigr{)}\mathcal{O}^{\rho\sigma}\,.$
(82)

Plugging this into the Lorentz decomposition Eq. (22) for the pion one finds

$\displaystyle\langle\bm{p}|\mathcal{O}_{\rm v2a}^{i}|\bm{p}\rangle$
$\displaystyle=2v^{q}_{2}\,Ep^{i}\,,$ (83)
$\displaystyle\langle\bm{p}|\mathcal{O}_{\rm v2b}|\bm{p}\rangle$
$\displaystyle=2v^{q}_{2}\,\frac{4E^{2}-m^{2}}{3}\,.$ (84)

Using Eq. (23), we can show for the rho that

$\displaystyle\langle\bm{p},\lambda|\mathcal{O}_{\rm
v2a}^{i}|\bm{p},\lambda\rangle$
$\displaystyle=2Ep^{i}\begin{cases}a^{q}_{2}+\frac{2}{3}d^{q}_{2}&\text{ for
}\lambda=0\\\ a^{q}_{2}-\frac{1}{3}d^{q}_{2}&\text{ for
}\lambda=\pm\end{cases}\,,$ (85)
$\displaystyle\begin{split}\langle\bm{p},\lambda|\mathcal{O}_{\rm
v2b}|\bm{p},\lambda\rangle&=2\Bigl{(}a^{q}_{2}-\frac{d^{q}_{2}}{3}\Bigr{)}\frac{4E^{2}-m^{2}}{3}\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\\\
&\quad+2\frac{d^{q}_{2}}{3}\begin{cases}4E^{2}-3m^{2}&\text{ for }\lambda=0\\\
m^{2}&\text{ for }\lambda=\pm\end{cases}\,.\end{split}$ (86)

In the actual computations we use the explicit operators

$\displaystyle\mathcal{O}_{\rm v2a}^{i}$
$\displaystyle=O^{\\{0i\\}},\quad\quad\,\,\,\,\text{with }i=1,2,3,$ (87)
$\displaystyle\mathcal{O}_{\rm v2b}$
$\displaystyle=O^{00}+\frac{1}{3}\bigl{(}O^{11}+O^{22}+O^{33}\bigr{)},$ (88)

where $O^{\mu\nu}$ is defined as

$\displaystyle
O_{\mu\nu}=\frac{i}{2}\,\mathrm{bar}{q}\,\gamma^{\mu}\,\overleftrightarrow{D}^{\nu}\,q.$
(89)

The conversion between Minkowski and Euclidean convention is finally given by

$\displaystyle\mathcal{O}_{\rm v2a}^{(E),i}=-i\,\mathcal{O}_{\rm
v2a}^{i}\,\,,$ $\displaystyle\mathcal{O}_{\rm v2b}^{(E)}=-\mathcal{O}_{\rm
v2b}\,\,,$ (90)

see ref. [18].

##  Appendix C Numerical methods

###  C.1 Stochastic propagator

Figure 9: Sketch of a generic meson three-point function in forward and
backward direction. The source timeslice is $r_{4}$, the backward/forward sink
timeslice $x_{4}^{\prime}$/$x_{4}^{\prime\prime}$ and the current is located
at timeslice $y_{4}$. While the solid lines represent point-to-all propagators
the wiggly line illustrates the stochastic timeslice-to-all propagator
connecting the sink and the operator insertion.

Using the common sequential source method [99] the computational cost of
evaluating three-point functions is high because a new inversion is necessary
for each sink setup (timeslice, momentum and interpolating current). To reduce
the computational cost and maximize synergies in calculating matrix elements
we implemented a stochastic algorithm [100, 101] which circumvents this
limitations. The implementation we propose parallelizes the computations in
such a way that those for multiple source positions and multiple insertion
positions can be done simultaneously. A similar approach was already used in
[102, 103, 104, 105], however, by storing the uncontracted data, with all spin
indices open, on disk, our implementation enables the user to analyze any
channel of interest at a later stage.

First of all we factorize the three-point correlation function into two
largely independent parts denoted as spectator $S$ and insertion $I$ part
which can be computed separately. The generic expression of our factorized
meson three-point function with open spin indices (Greek letters) as well as
color, stochastic and flavor indices reads

$\displaystyle\begin{split}&C_{\rm{3pt}}(\bm{p}^{\prime\\!},\bm{q},x_{4}^{\prime},y_{4})_{f_{1},f_{2},f_{3}}^{\alpha^{\prime\\!}\beta^{\prime\\!}\tilde{\alpha}\tilde{\beta}\beta\alpha}=\Gamma_{\rm{snk}}^{\alpha^{\prime\\!}\beta^{\prime\\!}}\,\Gamma_{\rm{ins}}^{\tilde{\alpha}\tilde{\beta}}\,\Gamma_{\rm{src}}^{\beta\alpha}\\\
&\times\frac{1}{N_{\rm{sto}}}\sum_{i=1}^{N_{\rm{sto}}}S_{i,f_{1}}(\bm{p}^{\prime\\!},x^{\prime}_{4})^{\beta^{\prime\\!}\alpha^{\prime\\!}\alpha}_{a}\,I_{i,f_{2},f_{3}}(\bm{q},y_{4})^{\tilde{\alpha}\tilde{\beta}\beta}_{a}\,.\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{split}$ (91)

$\Gamma_{\rm{src/snk}}$ corresponds to the interpolating currents of the
(smeared) meson source or sink and $\Gamma_{\rm{ins}}$ corresponds to the
local operator insertion which can contain additional derivatives (at the
moment only first derivatives are implemented). We define the spectator
$S_{i,f_{1}}(\bm{p}^{\prime},x^{\prime}_{4})^{\beta^{\prime}\alpha^{\prime}\alpha}_{a}$
and the insertion
$I_{i,f_{2},f_{3}}(\bm{q},y_{4})^{\tilde{\alpha}\tilde{\beta}\beta}_{a}$ parts
as

$\displaystyle\begin{split}&S_{i,f_{1}}(\bm{p}^{\prime\\!},x^{\prime}_{4})^{\beta^{\prime\\!}\alpha^{\prime\\!}\alpha}_{a}=\\\
&\sum_{\bm{x}^{\prime}}\delta_{a^{\prime\\!}b^{\prime\\!}}\Bigl{[}\eta_{i}(x^{\prime})\gamma_{5}\Bigr{]}^{\beta^{\prime\\!}}_{b^{\prime\\!}}\,\Bigl{[}\gamma_{5}G_{f_{1}}^{\dagger}(x^{\prime},r)\gamma_{5}\Bigr{]}^{\alpha^{\prime\\!}\alpha}_{a^{\prime\\!}a}\,e^{-i\bm{p}^{\prime\\!}\cdot\bm{x}^{\prime\\!}}\,,\\\
\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip
0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{split}$ (92)
$\displaystyle\begin{split}&I_{i,f_{2},f_{3}}(\bm{q},y_{4})^{\tilde{\alpha}\tilde{\beta}\beta}_{a}=\\\
&\sum_{\bm{y}}\delta_{ab}\delta_{\tilde{a}\tilde{b}}\Bigl{[}\gamma_{5}s_{i,\,f_{2}}(y)\Bigr{]}^{*\tilde{\alpha}}_{\phantom{*}\tilde{a}}\,G_{f_{3}}(y,r)^{\tilde{\beta}\beta}_{\tilde{b}b}\,e^{i\bm{q}\cdot\bm{y}}\,,\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{split}$ (93)

where we assume that the spatial source is located at the origin without loss
of generality. In Fig. 9 we show a sketch of a generic meson three-point
function to further illustrate the factorization. We have two quark
propagators $G_{f_{i}}$ in Eqs. (92) and (93) depicted as solid lines
connecting the source position $r$ with all other points of the lattice. These
point-to-all propagators are computed using the solver methods introduced in
Sec. III.1. The third propagator connecting the sink timeslice with the
insertion current is plotted as a wiggly line and estimated by

$\displaystyle
G_{f_{2}}(y,x^{\prime})^{\tilde{\alpha}\beta^{\prime\\!}}_{\tilde{a}b^{\prime\\!}}$
$\displaystyle\approx\frac{1}{N}\sum_{i=1}^{N_{\rm{sto}}}s_{i,\,f_{2}}(y)^{\tilde{\alpha}}_{\tilde{a}}\,\,\eta_{i}^{*}(x^{\prime})^{\beta^{\prime\\!}}_{b^{\prime\\!}},$
(94)

where the sum runs over $N_{\rm{sto}}$ realizations of the noise vector
$\eta_{i}(x^{\prime})$, with the properties

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\eta_{i}(x)^{\alpha}_{a}$
$\displaystyle=0+\mathcal{O}\left(\frac{1}{\sqrt{N}}\right)\,,$ (95)
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\eta_{i}(x)^{\alpha}_{a}\eta_{i}^{*}(x^{\prime})^{\alpha^{\prime}}_{a^{\prime}}$
$\displaystyle=\delta_{xx^{\prime}}\delta_{\alpha\alpha^{\prime}}\delta_{aa^{\prime}}+\mathcal{O}\left(\frac{1}{\sqrt{N}}\right)\,.$
(96)

In the implementation presented in this work we use time partitioned
$\mathbb{Z}_{2}$ [106] noise vectors $\eta_{i}(x)$ which are set to zero
unless $x_{4}=x_{4}^{\prime\prime}$ or $x_{4}=x_{4}^{\prime}$. Seeding the
noise vectors in forward and backward temporal direction enables us to
increase statistics by a factor of two with only little computational
overhead. Moreover, the insertion part of our factorization is constructed
such that it can be reused in the calculation of baryon three-point functions
by contracting an appropriate spectator part. The results in [107] for the
baryon computation look very promising. In this work we showed that this holds
for the meson computations as well and that the stochastic approach is a
serious alternative to the sequential source method. For more details about
the implementation the interested reader is referred to [100].

###  C.2 Disconnected contributions

Figure 10: Sketch of a generic meson disconnected three-point function. The
source timeslice is $r_{4}$, the sink timeslice $x_{4}^{\prime}$ and the
current (loop) is located at timeslice $y_{4}$. While the solid lines
represent point-to-all propagators the solid circle represents a quark loop.
Figure 11: Extrapolation for the flavor singlet operator
($\mathrm{bar}{u}u+\mathrm{bar}{d}d+\mathrm{bar}{s}s$) and the flavor non-
singlet operator ($\mathrm{bar}{u}u+\mathrm{bar}{d}d-2\mathrm{bar}{s}s$) using
the fits shown in, e.g., Fig. 5. The three different columns show the
extrapolation for the lattice spacing dependence along the ${\rm
TrM}=\text{const.}$, $m_{s}=\text{const.}$, and the $m_{l}=m_{s}$ trajectories
for the reduced matrix elements $v_{2}$, $a_{2}$, and $d_{2}$.

In addition to connected contributions to the three-point function, which we
treated in Sec. C.1, we also compute the disconnected contributions
illustrated in Fig. 10. To get the disconnected contribution we multiply the
two-point function by a disonnected loop $L(\tau)$ [108] which reads

$\displaystyle
C^{\rm{discon}}_{\rm{3pt}}(\bm{p}^{\prime\\!},\bm{q},x_{4}^{\prime},y_{4})=$
$\displaystyle\,\langle
C^{\rm{c}}_{\rm{2pt}}(x_{4}^{\prime},r_{4})L^{\rm{c}}(y_{4})\rangle_{\rm{c}}-\langle
C^{\rm{c}}_{\rm{2pt}}(x_{4}^{\prime},r_{4})\rangle_{\rm{c}}\langle
L^{\rm{c}}(y_{4})\rangle_{\rm{c}},$ (97)

where $\langle\rangle_{\rm{c}}$ makes the configuration average explicit. The
loop is given by

$\displaystyle
L^{\rm{c}}(y_{4})=\sum_{\bm{y}}\,\operatorname{Tr}[G_{f}(y,y)\Gamma_{\rm{ins}}],$
(98)

where $G_{f}$ denotes a quark propagator of flavor $f$ with
$y=(y_{4},\bm{y})$. Note that we set $r_{4}=0$ without loss of generality. To
compute the propagator in Eq. (98) we use stochastic estimators similar to the
approach presented in the last section and the corresponding solvers
introduced in Sec. III.1. More details on our approach of the computation of
disconnected loops are given in [109].

##  Appendix D Additional plots

In Fig. 11 we explicitly show the $a$ dependence of the global extrapolation
plots for the reduced matrix elements $v_{2}$, $a_{2}$, and $d_{2}$. We still
use physical masses but in contrast to the main text we do not average over
equal lattice spacings at this point.

## References

  * Sullivan [1972] J. D. Sullivan, One-Pion Exchange and Deep-Inelastic Electron-Nucleon Scattering, Phys. Rev. D 5, 1732 (1972).
  * Betev  _et al._ [1985] B. Betev _et al._ (NA10), Differential Cross-Section of High Mass Muon Pairs Produced by a $194\text{\,}\mathrm{G}\mathrm{eV}\mathrm{/}\mathrm{c}$ $\pi^{-}$ Beam on a Tungsten Target, Z. Phys. C 28, 9 (1985).
  * Greenlee  _et al._ [1985] H. B. Greenlee _et al._ , The Production of Massive Muon Pairs in $\pi^{-}$-Nucleus Collisions, Phys. Rev. Lett. 55, 1555 (1985).
  * Conway  _et al._ [1989] J. S. Conway _et al._ , Experimental study of muon pairs produced by 252-GeV pions on tungsten, Phys. Rev. D 39, 92 (1989).
  * Aghasyan  _et al._ [2017] M. Aghasyan _et al._ (COMPASS), First Measurement of Transverse-Spin-Dependent Azimuthal Asymmetries in the Drell-Yan Process, Phys. Rev. Lett. 119, 112002 (2017), arXiv:1704.00488 [hep-ex] .
  * Adams  _et al._ [2018] B. Adams _et al._ , Letter of Intent: A New QCD facility at the M2 beam line of the CERN SPS (COMPASS++/AMBER), (2018), arXiv:1808.00848 [hep-ex] .
  * CERN [2021] CERN, Amber – a new qcd facility at the m2 beam line of the cern sps, https://nqf-m2.web.cern.ch/ (2021), [Online; accessed 29-July-2021].
  * Chekanov  _et al._ [2002] S. Chekanov _et al._ (ZEUS), Leading neutron production in $e^{+}$ $p$ collisions at HERA, Nucl. Phys. B 637, 3 (2002), arXiv:hep-ex/0205076 .
  * Aaron  _et al._ [2010] F. D. Aaron _et al._ (H1), Measurement of leading neutron production in deep-inelastic scattering at HERA, Eur. Phys. J. C 68, 381 (2010), arXiv:1001.0532 [hep-ex] .
  * Holtmann  _et al._ [1994] H. Holtmann, G. Levman, N. N. Nikolaev, A. Szczurek, and J. Speth, How to measure the pion structure function at HERA, Phys. Lett. B 338, 363 (1994).
  * McKenney  _et al._ [2016] J. R. McKenney, N. Sato, W. Melnitchouk, and C.-R. Ji, Pion structure function from leading neutron electroproduction and SU(2) flavor asymmetry, Phys. Rev. D 93, 054011 (2016), arXiv:1512.04459 [hep-ph] .
  * Barry  _et al._ [2018] P. C. Barry, N. Sato, W. Melnitchouk, and C.-R. Ji, First Monte Carlo Global QCD Analysis of Pion Parton Distributions, Phys. Rev. Lett. 121, 152001 (2018), arXiv:1804.01965 [hep-ph] .
  * Montgomery  _et al._ [2017] R. A. Montgomery, J. R. M. Annand, D. Dutta, C. E. Keppel, P. King, B. Wojtsekhowski, and J. Zhang (TDIS, SBS), Proposed Measurement of Tagged Deep Inelastic Scattering in Hall A of Jefferson lab, AIP Conf. Proc. 1819, 030004 (2017).
  * Keppel  _et al._ [2015] C. E. Keppel _et al._ , Measurement of Tagged Deep Inelastic Scattering (TDIS): Hall A and SBS Collaboration Proposal (2015), PR12-15-006.
  * Aguilar  _et al._ [2019] A. C. Aguilar _et al._ , Pion and kaon structure at the electron-ion collider, Eur. Phys. J. A 55, 190 (2019), arXiv:1907.08218 [nucl-ex] .
  * Zyla  _et al._ [2020] P. A. Zyla _et al._ (PDG), Review of Particle Physics, Prog. Theor. Exp. Phys. 2020, 083C01 (2020).
  * Hoodbhoy  _et al._ [1989] P. Hoodbhoy, R. L. Jaffe, and A. Manohar, Novel effects in deep inelastic scattering from spin-one hadrons, Nucl. Phys. B 312, 571 (1989).
  * Best  _et al._ [1997a] C. Best, M. Göckeler, R. Horsley, E.-M. Ilgenfritz, H. Perlt, P. Rakow, A. Schäfer, G. Schierholz, A. Schiller, and S. Schramm, $\pi$ and $\rho$ structure functions from lattice QCD, Phys. Rev. D 56, 2743 (1997a), arXiv:hep-lat/9703014 .
  * Ji [2021] X.-D. Ji, private communication (2021), [July-2021].
  * Airapetian  _et al._ [2005] A. Airapetian _et al._ (HERMES), Measurement of the Tensor Structure Function b(1) of the Deuteron, Phys. Rev. Lett. 95, 242001 (2005), arXiv:hep-ex/0506018 .
  * Close and Kumano [1990] F. E. Close and S. Kumano, Sum rule for the spin-dependent structure function $b_{1}(x)$ for spin-one hadrons, Phys. Rev. D 42, 2377 (1990).
  * Kumano and Song [2017] S. Kumano and Q.-T. Song, Estimate on Spin Asymmetry for Drell-Yan Process at Fermilab with Tensor-Polarized Deuteron, JPS Conf. Proc. 13, 020048 (2017), arXiv:1611.00474 [hep-ph] .
  * Alleda  _et al._ [2013] K. Alleda _et al._ , The Deuteron Tensor Structure Function $b_{1}$: A Proposal to Jefferson Lab PAC-40 (2013), PR12-13-011.
  * Song [2019] Q.-T. Song, Structure of Deuteron by Polarized Proton-Deuteron Drell-Yan Process, JPS Conf. Proc. 26, 031001 (2019).
  * Boer [2019] D. Boer, Overview of Spin Physics at EIC, PoS SPIN2018, 167 (2019), arXiv:1903.01119 [hep-ph] .
  * Martinelli and Sachrajda [1987] G. Martinelli and C. T. Sachrajda, Pion structure functions from lattice QCD, Phys. Lett. B 196, 184 (1987).
  * Martinelli and Sachrajda [1988] G. Martinelli and C. T. Sachrajda, A lattice calculation of the pion’s form-factor and structure function, Nucl. Phys. B 306, 865 (1988).
  * Best  _et al._ [1997b] C. Best _et al._ , Hadron structure functions from lattice QCD - 1997, AIP Conf. Proc. 407, 997 (1997b), arXiv:hep-ph/9706502 .
  * Best  _et al._ [1998] C. Best, M. Göckeler, R. Horsley, H. Perlt, P. Rakow, A. Schäfer, G. Schierholz, A. Schiller, and S. Schramm, The Deep-Inelastic Structure Functions of $\pi$ and $\rho$ Mesons, Nucl. Phys. B (Proc. Suppl.) 63A-C, 236 (1998), arXiv:hep-lat/9710037 .
  * Guagnelli  _et al._ [2005] M. Guagnelli, K. Jansen, F. Palombi, R. Petronzio, A. Shindler, and I. Wetzorke (ZeRo), Non-perturbative pion matrix element of a twist-2 operator from the lattice, Eur. Phys. J. C 40, 69 (2005), arXiv:hep-lat/0405027 .
  * Capitani  _et al._ [2006] S. Capitani, K. Jansen, M. Papinutto, A. Shindler, C. Urbach, and I. Wetzorke, Parton distribution functions with twisted mass fermions, Phys. Lett. B 639, 520 (2006), arXiv:hep-lat/0511013 .
  * Brömmel  _et al._ [2008] D. Brömmel _et al._ (QCDSF/UKQCD), Quark distributions in the pion, PoS LATTICE2007, 140 (2008).
  * Bali  _et al._ [2014a] G. S. Bali, S. Collins, B. Glässle, M. Göckeler, N. Javadi-Motaghi, J. Najjar, W. Söldner, and A. Sternbeck, Pion structure from lattice QCD, PoS LATTICE2013, 447 (2014a), arXiv:1311.7639 [hep-lat] .
  * Abdel-Rehim  _et al._ [2015] A. Abdel-Rehim _et al._ , Nucleon and pion structure with lattice QCD simulations at physical value of the pion mass, Phys. Rev. D 92, 114513 (2015), [Erratum: Phys. Rev. D 93, 039904(E) (2016)], arXiv:1507.04936 [hep-lat] .
  * Oehm  _et al._ [2019] M. Oehm, C. Alexandrou, M. Constantinou, K. Jansen, G. Koutsou, B. Kostrzewa, F. Steffens, C. Urbach, and S. Zafeiropoulos, $\langle x\rangle$ and $\langle x^{2}\rangle$ of the pion PDF from lattice QCD with $N_{f}=2+1+1$ dynamical quark flavors, Phys. Rev. D 99, 014508 (2019), arXiv:1810.09743 [hep-lat] .
  * Alexandrou  _et al._ [2020] C. Alexandrou, S. Bacchio, I. Clöet, M. Constantinou, K. Hadjiyiannakou, G. Koutsou, and C. Lauer, The Mellin moments $\langle x\rangle$ and $\langle x^{2}\rangle$ for the pion and kaon from lattice QCD, (2020), arXiv:2010.03495 [hep-lat] .
  * Alexandrou  _et al._ [2021] C. Alexandrou, S. Bacchio, I. Cloët, M. Constantinou, K. Hadjiyiannakou, G. Koutsou, and C. Lauer (ETM), The pion and kaon $\langle x^{3}\rangle$ from lattice QCD and PDF reconstruction from Mellin moments, (2021), arXiv:2104.02247 [hep-lat] .
  * Braun and Müller [2008] V. Braun and D. Müller, Exclusive processes in position space and the pion distribution amplitude, Eur. Phys. J. C 55, 349 (2008), arXiv:0709.1348 [hep-ph] .
  * Ji [2013] X. Ji, Parton Physics on a Euclidean Lattice, Phys. Rev. Lett. 110, 262002 (2013), arXiv:1305.1539 [hep-ph] .
  * Ma and Qiu [2018] Y.-Q. Ma and J.-W. Qiu, Extracting Parton Distribution Functions from Lattice QCD Calculations, Phys. Rev. D 98, 074021 (2018), arXiv:1404.6860 [hep-ph] .
  * Radyushkin [2017] A. V. Radyushkin, Quasi-parton distribution functions, momentum distributions, and pseudo-parton distribution functions, Phys. Rev. D 96, 034025 (2017), arXiv:1705.01488 [hep-ph] .
  * Ji  _et al._ [2021] X. Ji, Y.-S. Liu, Y. Liu, J.-H. Zhang, and Y. Zhao, Large-momentum effective theory, Rev. Mod. Phys. 93, 035005 (2021), arXiv:2004.03543 [hep-ph] .
  * Hua  _et al._ [2021] J. Hua, M.-H. Chu, P. Sun, W. Wang, J. Xu, Y.-B. Yang, J.-H. Zhang, and Q.-A. Zhang (Lattice Parton), Distribution Amplitudes of K* and $\phi$ at the Physical Pion Mass from Lattice QCD, Phys. Rev. Lett. 127, 062002 (2021), arXiv:2011.09788 [hep-lat] .
  * Braun  _et al._ [2017] V. M. Braun _et al._ , The $\rho$-meson light-cone distribution amplitudes from lattice QCD, JHEP 04, 082, arXiv:1612.02955 [hep-lat] .
  * Cichy [2021] K. Cichy, Progress in $x$-dependent partonic distributions from lattice QCD, in  _38th International Symposium on Lattice Field Theory_ (2021) arXiv:2110.07440 [hep-lat] .
  * Sufian  _et al._ [2019] R. S. Sufian, J. Karpie, C. Egerer, K. Orginos, J.-W. Qiu, and D. G. Richards, Pion valence quark distribution from matrix element calculated in lattice QCD, Phys. Rev. D 99, 074507 (2019), arXiv:1901.03921 [hep-lat] .
  * Sufian  _et al._ [2020] R. S. Sufian, C. Egerer, J. Karpie, R. G. Edwards, B. Joó, Y.-Q. Ma, K. Orginos, J.-W. Qiu, and D. G. Richards, Pion valence quark distribution from current-current correlation in lattice QCD, Phys. Rev. D 102, 054508 (2020), arXiv:2001.04960 [hep-lat] .
  * Zhang  _et al._ [2019] J.-H. Zhang, J.-W. Chen, L. Jin, H.-W. Lin, A. Schäfer, and Y. Zhao, First direct lattice-QCD calculation of the $x$-dependence of the pion parton distribution function, Phys. Rev. D 100, 034505 (2019), arXiv:1804.01483 [hep-lat] .
  * Izubuchi  _et al._ [2019] T. Izubuchi, L. Jin, C. Kallidonis, N. Karthik, S. Mukherjee, P. Petreczky, C. Shugert, and S. Syritsyn, Valence parton distribution function of pion from fine lattice, Phys. Rev. D 100, 034516 (2019), arXiv:1905.06349 [hep-lat] .
  * Joó  _et al._ [2019] B. Joó, J. Karpie, K. Orginos, A. V. Radyushkin, D. G. Richards, R. S. Sufian, and S. Zafeiropoulos, Pion valence structure from Ioffe-time parton pseudodistribution functions, Phys. Rev. D 100, 114512 (2019), arXiv:1909.08517 [hep-lat] .
  * Manohar [1992] A. V. Manohar, An introduction to spin dependent deep inelastic scattering, (1992), arXiv:hep-ph/9204208 .
  * Ji [1998] X.-D. Ji, Off forward parton distributions, J. Phys. G 24, 1181 (1998), arXiv:hep-ph/9807358 .
  * Dokshitzer [1977] Y. L. Dokshitzer, Calculation of the structure functions for deep inelastic scattering and $e^{+}e^{-}$ annihilation by perturbation theory in quantum chromodynamics, Sov. Phys. JETP 46, 641 (1977).
  * Gribov and Lipatov [1972] V. N. Gribov and L. N. Lipatov, Deep inelastic $ep$ scattering in perturbation theory, Sov. J. Nucl. Phys. 15, 438 (1972).
  * Altarelli and Parisi [1977] G. Altarelli and G. Parisi, Asymptotic freedom in parton language, Nucl. Phys. B 126, 298 (1977).
  * Diehl [2003] M. Diehl, Generalized parton distributions, Phys. Rept. 388, 41 (2003), arXiv:hep-ph/0307382 .
  * Callan and Gross [1969] J. Callan, C. G. and D. J. Gross, High-energy electroproduction and the constitution of the electric current, Phys. Rev. Lett. 22, 156 (1969).
  * Ji [1997] X. Ji, Gauge-Invariant Decomposition of Nucleon Spin, Phys. Rev. Lett. 78, 610 (1997), arXiv:hep-ph/9603249 .
  * Bruno  _et al._ [2015] M. Bruno _et al._ , Simulation of QCD with $N_{f}=2+1$ flavors of non-perturbatively improved Wilson fermions, JHEP 02, 043, arXiv:1411.3982 [hep-lat] .
  * Lüscher and Schaefer [2013] M. Lüscher and S. Schaefer, Lattice QCD with open boundary conditions and twisted-mass reweighting, Comput. Phys. Commun. 184, 519 (2013), arXiv:1206.2809 [hep-lat] .
  * Lüscher and Schaefer [2011] M. Lüscher and S. Schaefer, Lattice QCD without topology barriers, JHEP 07, 036, arXiv:1105.4749 [hep-lat] .
  * Bali  _et al._ [2016a] G. S. Bali, E. E. Scholz, J. Simeth, and W. Söldner (RQCD), Lattice simulations with $N_{f}=2+1$ improved Wilson fermions at a fixed strange quark mass, Phys. Rev. D 94, 074501 (2016a), arXiv:1606.09039 [hep-lat] .
  * Bali  _et al._ [2020a] G. S. Bali, L. Barca, S. Collins, M. Gruber, M. Löffler, A. Schäfer, W. Söldner, P. Wein, S. Weishäupl, and T. Wurm (RQCD), Nucleon axial structure from lattice QCD, JHEP 05, 126, arXiv:1911.13150 [hep-lat] .
  * Babich  _et al._ [2010] R. Babich, J. Brannick, R. C. Brower, M. A. Clark, T. A. Manteuffel, S. F. McCormick, J. C. Osborn, and C. Rebbi, Adaptive multigrid algorithm for the lattice Wilson-Dirac operator, Phys. Rev. Lett. 105, 201602 (2010), arXiv:1005.3043 [hep-lat] .
  * Frommer  _et al._ [2014] A. Frommer, K. Kahl, S. Krieg, B. Leder, and M. Rottmann, Adaptive aggregation based domain decomposition multigrid for the lattice Wilson Dirac operator, SIAM J. Sci. Comput. 36, A1581 (2014), arXiv:1303.1377 [hep-lat] .
  * Heybrock  _et al._ [2016] S. Heybrock, M. Rottmann, P. Georg, and T. Wettig, Adaptive algebraic multigrid on SIMD architectures, PoS LATTICE2015, 036 (2016), arXiv:1512.04506 [physics.comp-ph] .
  * Richtmann  _et al._ [2016] D. Richtmann, S. Heybrock, and T. Wettig, Multiple right-hand-side setup for the DD-$\alpha$AMG, PoS LATTICE2015, 035 (2016), arXiv:1601.03184 [hep-lat] .
  * Georg  _et al._ [2017] P. Georg, D. Richtmann, and T. Wettig, pMR: A high-performance communication library, PoS LATTICE2016, 361 (2017), arXiv:1701.08521 [hep-lat] .
  * Georg  _et al._ [2018] P. Georg, D. Richtmann, and T. Wettig, DD-$\alpha$AMG on QPACE 3, EPJ Web Conf. 175, 02007 (2018), arXiv:1710.07041 [hep-lat] .
  * Lüscher [2007a] M. Lüscher, Local coherence and deflation of the low quark modes in lattice QCD, JHEP 07, 081, arXiv:0706.2298 [hep-lat] .
  * Lüscher [2007b] M. Lüscher, Deflation acceleration of lattice QCD simulations, JHEP 12, 011, arXiv:0710.5417 [hep-lat] .
  * Güsken [1990] S. Güsken, A study of smearing techniques for hadron correlation functions, Nucl. Phys. B (Proc. Suppl.) 17, 361 (1990).
  * Falcioni  _et al._ [1985] M. Falcioni, M. L. Paciello, G. Parisi, and B. Taglienti, Again on SU(3) glueball mass, Nucl. Phys. B 251, 624 (1985).
  * Edwards and Joó [2005] R. G. Edwards and B. Joó (SciDAC, LHPC, UKQCD), The Chroma Software System for Lattice QCD, Nucl. Phys. B (Proc. Suppl.) 140, 832 (2005), arXiv:hep-lat/0409003 .
  * Alexandrou  _et al._ [2017] C. Alexandrou, M. Constantinou, K. Hadjiyiannakou, K. Jansen, H. Panagopoulos, and C. Wiese, Gluon momentum fraction of the nucleon from lattice QCD, Phys. Rev. D 96, 054503 (2017), arXiv:1611.06901 [hep-lat] .
  * Shanahan and Detmold [2019] P. E. Shanahan and W. Detmold, Gluon gravitational form factors of the nucleon and the pion from lattice QCD, Phys. Rev. D 99, 014511 (2019), arXiv:1810.04626 [hep-lat] .
  * Göckeler  _et al._ [1996] M. Göckeler, R. Horsley, E.-M. Ilgenfritz, H. Perlt, P. E. L. Rakow, G. Schierholz, and A. Schiller, Lattice operators for moments of the structure functions and their transformation under the hypercubic group, Phys. Rev. D 54, 5705 (1996), arXiv:hep-lat/9602029 .
  * Bali  _et al._ [2021] G. S. Bali, S. Bürger, S. Collins, M. Göckeler, M. Gruber, S. Piemonte, A. Schäfer, A. Sternbeck, and P. Wein, Nonperturbative Renormalization in Lattice QCD with three Flavors of Clover Fermions: Using Periodic and Open Boundary Conditions, Phys. Rev. D 103, 094511 (2021), arXiv:2012.06284 [hep-lat] .
  * Werner  _et al._ [2020] M. Werner _et al._ (ETMC), Hadron–Hadron interactions from $N_{f}=2+1+1$ lattice QCD: the $\rho$-resonance, Eur. Phys. J. A56, 61 (2020), arXiv:1907.01237 [hep-lat] .
  * Göckeler  _et al._ [2012] M. Göckeler, R. Horsley, M. Lage, U.-G. Meißner, P. E. L. Rakow, A. Rusetsky, G. Schierholz, and J. M. Zanotti, Scattering phases for meson and baryon resonances on general moving-frame lattices, Phys. Rev. D86, 094513 (2012), arXiv:1206.4141 [hep-lat] .
  * Erben  _et al._ [2020] F. Erben, J. R. Green, D. Mohler, and H. Wittig, Rho resonance, timelike pion form factor, and implications for lattice studies of the hadronic vacuum polarization, Phys. Rev. D101, 054504 (2020), arXiv:1910.01083 [hep-lat] .
  * Lüscher [1986] M. Lüscher, Volume Dependence of the Energy Spectrum in Massive Quantum Field Theories. 2. Scattering States, Commun. Math. Phys. 105, 153 (1986).
  * Lüscher [1991] M. Lüscher, Two-particle states on a torus and their relation to the scattering matrix, Nucl. Phys. B354, 531 (1991).
  * Dudek  _et al._ [2013] J. J. Dudek, R. G. Edwards, and C. E. Thomas (HSC), Energy dependence of the $\rho$ resonance in $\pi\pi$ elastic scattering from lattice QCD, Phys. Rev. D87, 034505 (2013), [Erratum: Phys. Rev. D90 (2014) 099902(E)], arXiv:1212.0830 [hep-ph] .
  * Wilson  _et al._ [2015] D. J. Wilson, R. A. Briceño, J. J. Dudek, R. G. Edwards, and C. E. Thomas, Coupled $\pi\pi,K\mathrm{bar}{K}$ scattering in $P$-wave and the $\rho$ resonance from lattice QCD, Phys. Rev. D 92, 094502 (2015), arXiv:1507.02599 [hep-ph] .
  * Gounaris and Sakurai [1968] G. J. Gounaris and J. J. Sakurai, Finite-width corrections to the vector-meson-dominance prediction for $\rho\to e^{+}e^{-}$, Phys. Rev. Lett. 21, 244 (1968).
  * Kawarabayashi and Suzuki [1966] K. Kawarabayashi and M. Suzuki, Partially conserved axial-vector current and the decays of vector mesons, Phys. Rev. Lett. 16, 255 (1966).
  * Riazuddin and Fayyazuddin [1966] Riazuddin and Fayyazuddin, Algebra of current components and decay widths of $\rho$ and $K^{*}$ mesons, Phys. Rev. 147, 1071 (1966).
  * Djukanovic  _et al._ [2004] D. Djukanovic, M. R. Schindler, J. Gegelia, G. Japaridze, and S. Scherer, Universality of the $\rho$ meson coupling in effective field theory, Phys. Rev. Lett. 93, 122002 (2004), arXiv:hep-ph/0407239 .
  * Feng  _et al._ [2015] X. Feng, S. Aoki, S. Hashimoto, and T. Kaneko, Timelike pion form factor in lattice QCD, Phys. Rev. D91, 054504 (2015), arXiv:1412.6319 [hep-lat] .
  * Meyer [2011] H. B. Meyer, Lattice QCD and the Timelike Pion Form Factor, Phys. Rev. Lett. 107, 072002 (2011), arXiv:1105.1892 [hep-lat] .
  * Lüscher and Wolff [1990] M. Lüscher and U. Wolff, How to calculate the elastic scattering matrix in two-dimensional quantum field theories by numerical simulation, Nucl. Phys. B 339, 222 (1990).
  * Bali  _et al._ [2016b] G. S. Bali, S. Collins, A. Cox, G. Donald, M. Göckeler, C. B. Lang, and A. Schäfer (RQCD), $\rho$ and $K^{*}$ resonances on the lattice at nearly physical quark masses and $N_{f}=2$, Phys. Rev. D 93, 054509 (2016b), arXiv:1512.08678 [hep-lat] .
  * Fischer  _et al._ [2020a] M. Fischer, B. Kostrzewa, M. Mai, M. Petschlies, F. Pittler, M. Ueding, C. Urbach, and M. Werner (ETM), The $\rho$-resonance with physical pion mass from $N_{f}=2$ lattice QCD, (2020a), arXiv:2006.13805 [hep-lat] .
  * Fischer  _et al._ [2020b] M. Fischer, B. Kostrzewa, J. Ostmeyer, K. Ottnad, M. Ueding, and C. Urbach, On the generalised eigenvalue method and its relation to Prony and generalised pencil of function methods, Eur. Phys. J. A 56, 206 (2020b), arXiv:2004.10472 [hep-lat] .
  * [96] G. S. Bali _et al._ , Scale setting and the light hadron spectrum in $N_{f}=2+1$ QCD with Wilson fermions, to be published in preparation.
  * Jülich Supercomputing Centre [2018] Jülich Supercomputing Centre, JURECA: Modular supercomputer at Jülich Supercomputing Centre, Journal of large-scale research facilities 4, A132 (2018).
  * Jülich Supercomputing Centre [2019] Jülich Supercomputing Centre, HDF Cloud – Helmholtz Data Federation Cloud Resources at the Jülich Supercomputing Centre, Journal of Large-Scale Research Facilities 5, A137 (2019).
  * Martinelli and Sachrajda [1989] G. Martinelli and C. T. Sachrajda, A lattice study of nucleon structure, Nucl. Phys. B 316, 355 (1989).
  * Bali  _et al._ [2018] G. S. Bali, S. Collins, B. Gläßle, S. Heybrock, P. Korcyl, M. Löffler, R. Rödl, and A. Schäfer, Baryonic and mesonic 3-point functions with open spin indices, EPJ Web Conf. 175, 06014 (2018), arXiv:1711.02384 [hep-lat] .
  * Rödl [2020] R. H. Rödl,  _Nucleon generalized form factors from two-flavor lattice QCD and stochastic three-point functions with open indices_ , Ph.D. thesis, Regensburg U. (2020).
  * Evans  _et al._ [2010] R. Evans, G. S. Bali, and S. Collins, Improved semileptonic form factor calculations in lattice QCD, Phys. Rev. D 82, 094501 (2010), arXiv:1008.3293 [hep-lat] .
  * Alexandrou  _et al._ [2014] C. Alexandrou, M. Constantinou, S. Dinter, V. Drach, K. Jansen, K. Hadjiyiannakou, and D. B. Renner (ETM), A stochastic method for computing hadronic matrix elements, Eur. Phys. J. C 74, 2692 (2014), arXiv:1302.2608 [hep-lat] .
  * Bali  _et al._ [2014b] G. S. Bali _et al._ , Nucleon generalized form factors and sigma term from lattice QCD near the physical quark mass, PoS Lattice2013, 291 (2014b), arXiv:1312.0828 [hep-lat] .
  * Yang  _et al._ [2016] Y.-B. Yang, A. Alexandru, T. Draper, M. Gong, and K.-F. Liu, Stochastic method with low mode substitution for nucleon isovector matrix elements, Phys. Rev. D 93, 034503 (2016), arXiv:1509.04616 [hep-lat] .