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

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 420MeV down to 214MeV and with 5 different lattice spacings in the range of 0.1fm to 0.05fm. 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.


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 ρ 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., π +N → µ + +µ − +X, which are sensitive to the pion PDF at large x 0.15, and semi-inclusive (tagged) deep inelastic scattering (DIS), e.g., e + N → e + 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 ρ 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 * marius.loeffler@ur.de of the ρ 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 quarkgluon structure of resonances rather than only ground states, the ρ 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 ρ, 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 ρ 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 [18,[26][27][28][29][30][31] used quenched fermion representations, more recent simulations [32][33][34][35][36][37] use, for example, (cloverimproved) dynamical Wilson fermions, where the fermion determinant, and thus the quark sea, is taken into ac-count. 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 (ūu +dd − 2ss) and singlet (ūu +dd +ss) 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 ρ 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 D.
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 cal-culation, 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 where p is the three-momentum and λ 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 with the kinematic factors r µν , s µν , t µν , and u µν , which depend on the momentum transfer q, the target momentum p and the target polarization vector (cf. App. A for our conventions). The quantities r µν , s µν , t µν , and u µν are constructed such that they vanish upon averaging over the target spin, see ref. [17]. For spin-1 2 targets s µ corresponds to the spin four-vector. For spin-1 targets it corresponds to s µ = −i µνρσ e * ν e ρ p σ . Note that, due to current conservation any term proportional to q µ or q ν 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-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 as 1 where z is a lightlike vector with vanishing plus component, z µ = z − n µ − , where n µ − is dimensionless and can be used to project vectors onto their plus component, e.g., n − · 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 ↑ + q ↓ , while the PDF in Eq. (4) corresponds to the difference ∆q = q ↑ − q ↓ of the densities for quarks with opposite helicity. For spin-1 hadrons symmetry implies that distributions for different polarizations, λ = +, 0, −, are related [17,18] 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 In order to see the connection between PDFs and structure functions let us consider the operator product expansion 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 Λ QCD . For any general operator O µ1...µn d,n of dimension d and spin n one can show that the terms in the expansion have 1 In general one finds three quark and three gluon PDFs analogous to Eqs. (3) and (4), see, e.g., [52].
where M is the target hadron mass, ω = (2p · 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 each) a conventional basis for lowest twist t = 2 quark operators can be written in terms of two towers of operators 2 O µ1...µn where 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], 3 e.g., where we define the n-th moment of a function as 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 : In the following we will only write down the quark contribution. The antiquark contribution is obtained by simply substituting q →q. For spin-0 targets one obtains satisfying the Callan-Gross relation [57], F q 2 = 2xF q 1 + O(α 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 Considering the difference between targets with polariza-tion λ = ± and λ = 0 one finds 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 with the so-called reduced matrix element v n . Operators containing γ 5 do not contribute because of symmetry relations. For a spin-1 particle we find three independent structures p, λ|O µ1...µn where we use the convention that ε 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 ∆q.
In the structure functions always the sum of quark and antiquark contributions F q+q = F q + Fq 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 n v q n , n even , while, for spin-1 targets, 2M n (F q+q 1 ) = C (1) n a q n , n even , n a q n , n even , n d q n , n even , The C 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 v q n = x n−1 q+q , n even , v q n = x n−1 q−q , n odd .
For spin-1 hadrons we find (by substituting Eq. (23) into Eq. (11)) that a q n = 1 3 λ=±,0 x n−1 q λ +q λ , n even , i.e., a q n yields the polarization average, while d q n corresponds to the difference between hadrons with polarization λ = ± and λ = 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 ∆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.

A. 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 • TrM = const.: The trace of the quark mass matrix is kept constant near its physical value [59] • m s = 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. I. We use five different lattice spacings from 0.0497 fm up to 0.0984 fm and m π covers a range from ∼ 420 MeV down to ∼ 220 MeV with volumes Lm π between 3.8 and 6.4, see Tab. I. The two-and three-point functions introduced in Sec. III C are computed on the lattice using the gauge configurations in Tab. I. 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-α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 threepoint 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.

B. 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, 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ūu +dd − 2ss) 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 O(4) symmetry is reduced to that of its finite hypercubic subgroup 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 H(4), cf. refs. [18,77]. To be specific, we will use the operators O v2a = O 0i and O v2b = 4 3 O 00 , cf. App. B, where also an explicit definition of the operators is provided.
Our final results will be given in the MS scheme at a scale of 2 GeV. To this end, we adopt a two step procedure: First, we calculate the renormalization factors nonperturbatively in the RI /SMOM scheme. These are then converted to the MS scheme using perturbative QCD. The whole procedure is described in great detail in ref. [78], including subtleties due to the use of open I. CLS gauge ensembles analyzed in this work labeled by their identifier and sorted by the inverse coupling β and mπ. 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π, and the rho meson mass mρ computed in Sec. IV A. 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] fm. The number of configurations analyzed for the ensemble is denoted as n 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].

C. 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: We will consider pion (M = π) and rho mesons (M = ρ), where the Lorentz indices are only necessary in the latter case. The interpolating currentŌ 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 A) to enhance the ground state overlap. In addition to the two interpolating currents the three-point function contains an insertion current O at timeslice τ with 0 < τ < t. The extraction of the ground-state matrix element of O is the key task in the subsequent calculations. In this work we set t = t snk − t src , τ = t ins − t src and hence t src = 0 without loss of generality.

The pion
For the pion case we first define the matrix elements where Z π 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 where we assume that the same smearing setup is used at the source and the sink. For the ground state energies E π p we impose the continuum dispersion relation E π p = m 2 π + p 2 . Similarly, one can show that the spectral decomposition of the three-point function in Eq. (31) for large Euclidean times reads In practice it turns out that especially for the three-point functions the signal-to-noise ratio at large Euclidean time distances t and τ does not only contain the ground state contribution. How to exclude further excited state contributions is explained in Sec. III C 3.

The rho meson
For the rho mesons we define, in analogy to the pseudoscalar case Eq. (33), where the polarization vector µ obeys the general transversality condition Eq. (A3). Therefore, the insertion of a complete set of states into Eq. (30) (including a sum over all possible polarizations), yields 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 D. Inserting two complete sets of states into the threepoint function Eq. (31) we find Excited state contributions will be treated in Sec. III C 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 where ∆E 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 p, while the amplitudes in the three-point function, B 10 , B 01 , and B 11 , also depend on the operator insertion 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) twopion 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 D.

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 which holds for any operator insertion O in the threepoint 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., µ = ν = i in Eqs. (30) and (31)). Defining J p λ λ ≡ p, λ |O|p, λ /(2E ρ p ), one obtains 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 momentap = ±e i one finds the simple formulas for the extraction of the polarization-conserving matrix elements.

D. 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 π threshold. In the non-interacting case, their center of mass energies are given by where k and −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. 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 ⊗ I  [79].
defined as the direct product of the cubic group O (consisting of 24 rotations) and the group of space inversions I. 4 In a moving frame, however, the symmetry is reduced to the so-called little groups (for details see, e.g., refs. [79][80][81]) shown in Tab. III together with the decomposition into irreducible representations. The fourmomenta 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. [79][80][81]. 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 δ 1 is directly related to the quantized two-pion energy levels in finite volume. The latter appear when the condition is satisfied, see Fig. 1, which will be discussed in more detail below. The scattering phase shifts φ d Γ can be taken  from Tab. IV, using For the numeric evaluation of the generalized zeta function Z 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. IV with a phenomenological parametrization, where, for any given parametrization, we define the rho mass and width as [86] cot δ 1 using the Mandelstam variable s = E 2 cm . For instance one can use a relativistic Breit-Wigner (BW) ansatz where k = s/4 − m 2 π , as defined in Eq. (47), and Alternatively, one can use a Gounaris-Sakurai (GS) parametrization [86], where with If one chooses to apply the Kawarabayashi-Suzuki-Riazuddin-Fayyazuddin relation [87,88], m 2 ρ = 2g 2 ρππ F 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 π is well-known and the width of the rho is linked to the rho-pi-pi coupling constant g ρππ 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 × N t = 64 × 128, m π = 201 MeV, m ρ = 746 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 φ d Γ . 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 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 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 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. [81,[92][93][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.

A. 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 π are taken from [96] the values for m ρ 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 ∆E enter the fit as free parameters. However, using the ratio method introduced in Sec. III C 4 the additional amplitudes and also ∆E will not enter the results presented in this work.
The fits are performed using a constant fit window of ∼2 fm for all ensembles with open boundary conditions and we start 1 or 2 timeslices (∼ 0.1 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 cut ≈ 0.65 fm and fit only the single exponential Z ρ p (2E ρ p ) −1 e −E ρ p t for larger times. In case of periodic boundary conditions we choose a symmetric ansatz of the form Z ρ p (2E ρ p ) −1 (e −E ρ p t + e −E ρ p (T −t) ) with t cut ≤ t ≤ T − t 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 ρ p = m 2 ρ + 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.  L 2 ) projected to the rest frame using the continuum dispersion relation.

B. 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.
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 re- liable 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 C 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 q 2 (for the pion) as well as a q 2 and d q 2 for the rho meson. For on-axis momentap = ±e i we obtain (using Eqs. (B2) and (B4))  Fig. 5, i.e., we only plot the data points and individual fits for 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 Ov2a and O v2b using all possible momenta for the corresponding matrix element with n 2 ≤ 1 for the connected, disconnected non-singlet, and disconnected singlet contributions respectively.
with the operator O i v2a as inserted current. For O v2b we find Here we hide the superscripts π, ρ, and the subscript p for the mass and energy. Measurements using O i v2a always require nonzero momentum in direction i. If one uses O 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 n 2 = 1, where n is defined by p = 2π L 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 where the ratio R explicitly depends on p 2 , the operator O ∈ {O v2a , O v2b }, the sink timeslice t, and the insertion timeslice τ . For the operator combination O i v2a , with i = 1, 2, 3, we average the three spatial directions and fit to the data using the parameters B 0 , B 1 O v2a , n 2 = 1 and the excited state energy is given by ∆E n 2 =1 , independent of the operator combination. Note that we require nonzero momentum in direction i for O v2a . Additionally the operator combination O v2b gives rise to the further excited state amplitudes B 1 O v2b , n 2 and TABLE V. Summary of the occurrence of the individual fit parameters in the ansatz Eq. (64) for the extraction of v2 (pion) and a2 (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.
the excited state energy ∆E 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. I. A summary of the individual contributions is given in Tab. V. 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 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 ∆E n 2 =1 by an additional, simultaneous fit to the two-point function Eq. (39) in case of the reduced matrix element d 2 .
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. I. We can use this to our advantage by looking at the non-singlet (ūu +dd − 2ss) and singlet (ūu +dd +ss) 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, con-tent 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 τ . 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 τ for the given final time t.

C. Quark mass dependence and continuum extrapolation
As described in Sec. III A, 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 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 2 π , respectively. As discussed in Sec. IV B, the statistical error of the flavor singlet operator combinations is much larger than in the flavor nonsinglet 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 TrM = const., m s = 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 TrM = 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 β using the weighted where N β is the number of data points A i (ground state matrix elements) per β with corresponding errors σ 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. The final results for the reduced matrix elements are given in Tab. VI. 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 rangē m = (2m 2 K + m 2 π )/3 < 450 MeV and a < 0.09 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.

D. 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 + O(α s ) as a sum over the related reduced matrix element where we assume exact isospin symmetry and fs ≡ u + d + s is the flavor singlet while fns ≡ 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 In Tab. VII 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 only 5 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 B.
For the flavor non-singlet contributions given in Tab. VI 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 nonperturbatively, cf. Sec. III B. As a first step one can compare the connected results in Tab. VII to the connectedonly and the flavor non-singlet results in Tab. VI. 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 and F ρ 1 and the first moment of the structure function b 1 are in very good agreement with the flavor non-singlet results in Tab. VI and of course also reflect the connected only result in Tab. VI. As shown in Sec. IV B 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 SU(3) flavor symmetry for the quark sea, the results in Tab. VI and VII 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) π ∼ F 1 (x) ρ , 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 λ = 0 and λ = +/− rho meson. If the quarks were in a relative S-wave state, cf. the discussion in Sec. III D, 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 ∼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 , F ρ 1 , 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 twopion 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 D, 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. [81,[92][93][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.
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.    8. Extrapolation for the flavor singlet (ūu +dd +ss) and flavor non-singlet operator combinations (ūu +dd − 2ss) using the fits shown in, e.g., Fig. 5. From left to right we show the m 2 π dependence of the reduced matrix elements v2 (pion), a2, and d2 (both rho) for the three different trajectories we use in our analysis. Note that the symmetric trajectory m l = ms with exact SU(3) flavor symmetry approaches the chiral limit in the quark mass plane and not the physical point.  [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 andp the direction of the three-momentum. Then, where such that v T ⊥p. For the momentum we then have p ± = E ± |p| and p T = 0.
For the polarization vectors we use a dimensionless definition. They obey the general transversality condition where m is the hadron mass. For momenta in x direction the polarization vectors in the rest frame are given by (0, e λ ) with where e x 0 corresponds to the longitudinal polarization, while e x ± to the circular polarizations. For momenta in arbitrary direction, we have to rotate these vectors to to obtain the longitudinal polarization vector e 0 and the polarization vectors for the circular polarizations e ± . Last but not least, we have to perform a boost to the laboratory frame. This only affects e 0 (e ± are invariant because they are perpendicular to p), and we obtain In terms of the light-cone coordinates introduced above this yields ± (p, 0) = ±p ± /m , 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 H(4), cf. refs. [18,77]. To be specific, we will use the operators O i v2a = O 0i and O v2b = 4 3 O 00 . For the special case of two indices the action of the symmetrizing and trace-subtracting operator S is defined by Plugging this into the Lorentz decomposition Eq. (22) for the pion one finds Using Eq. (23), we can show for the rho that In the actual computations we use the explicit operators where O µν is defined as The conversion between Minkowski and Euclidean convention is finally given by see ref. [18].
Appendix C: Numerical methods

Stochastic propagator
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 Γ src/snk corresponds to the interpolating currents of the (smeared) meson source or sink and Γ 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,f1 (p , x 4 ) β α α a and the insertion I i,f2,f3 (q, y 4 )αβ β a parts as I i,f2,f3 (q, y 4 )αβ β a = y δ ab δãb γ 5 s i, f2 (y) * α a G f3 (y, r)β β bb e iq · y , 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 fi in Eqs. (C2) and (C3) 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 A. The third propagator connecting the sink timeslice with the insertion current is plotted as a wiggly line and estimated by where the sum runs over N sto realizations of the noise vector η i (x ), with the properties In the implementation presented in this work we use time partitioned Z 2 [106] noise vectors η i (x) which are set to zero unless x 4 = x 4 or x 4 = x 4 . 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].

Disconnected contributions
In addition to connected contributions to the threepoint 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(τ ) [108] which reads where G f denotes a quark propagator of flavor f with y = (y 4 , y). Note that we set r 4 = 0 without loss of generality. To compute the propagator in Eq. (C8) we use stochastic estimators similar to the approach presented in the last section and the corresponding solvers introduced in Sec. III A. More details on our approach of the computation of disconnected loops are given in [109].    11. Extrapolation for the flavor singlet operator (ūu +dd +ss) and the flavor non-singlet operator (ūu +dd − 2ss) using the fits shown in, e.g., Fig. 5. The three different columns show the extrapolation for the lattice spacing dependence along the TrM = const., ms = const., and the m l = ms trajectories for the reduced matrix elements v2, a2, and d2.