Quark and gluon momentum fractions in the pion from $N_f=2+1+1$ lattice QCD

We perform the first full decomposition of the pion momentum into its gluon and quark contributions. We employ an ensemble generated by the Extended Twisted Mass Collaboration with $N_f=2 + 1 +1$ Wilson twisted mass clover fermions at maximal twist tuned to reproduce the physical pion mass. We present our results in the $\overline{\mathrm{MS}}$ scheme at $2\gev$. We find $\avgx_{u+d}=0.601(28)$, $\avgx_s=0.059(13)$, $\avgx_c=0.019(05)$, and $\avgx_g=0.52(11)$ for the separate contributions, respectively, whose sum saturates the momentum sum rule.

(Dated: September 23, 2021) We perform the first full decomposition of the pion momentum into its gluon and quark contributions. We employ an ensemble generated by the Extended Twisted Mass Collaboration with N f = 2 + 1 + 1 Wilson twisted mass clover fermions at maximal twist tuned to reproduce the physical pion mass. We present our results in the MS scheme at 2 GeV. We find x u+d = 0.601(28), x s = 0.059 (13), x c = 0.019(05), and x g = 0.52 (11) for the separate contributions, respectively, whose sum saturates the momentum sum rule.
Introduction.-Quantum Chromodynamics (QCD) manifests itself in the form of a plethora of states -so called hadrons, formed by quarks and gluons. Pions are particularly interesting hadrons: they are the lightest and simplest of the QCD bound states composed out of quark and antiquark. At the same time pions are also pseudo-Goldstone bosons, with the spontaneous breaking of chiral symmetry playing a fundamental role in the emergence of their mass. Yet, in contrast to the nucleon (proton and neutron), a first principles computation of the pion structure, and in particular how quarks and gluons contribute to its mass and momentum decomposition is still lacking. The importance of this topic is well represented in the Electron Ion Collider (EIC) yellow report [1]: eight main science questions concerning pions (and kaons) are prominently put forward. Let us highlight two of these questions here: "What are the quark and gluon energy contributions to the pion mass?", and "Is the pion full or empty of gluons as viewed at large Q 2 ?" The results presented in this letter on the momentum decomposition of the pion using lattice QCD simulations address both questions. As mentioned before, there is a wealth of studies on the nucleon momentum decomposition available in the literature using phenomenological analyses of experimental data [2][3][4][5], and, more recently, from precise simulations of lattice QCD at the physical point [6,7]. The reason for the pion being much less well investigated is that proton and neutron structure is experimentally well accessible, while the pion is significantly more challenging because there is no pion target available. For that reason, only recently the first Monte Carlo global QCD analysis for pion PDFs has been presented in Ref. [8], which includes leading neutron electroproduction (LNE) data from HERA and Drell-Yan data from CERN and Fermilab. One of their interesting findings is that the decomposition of the pion momentum, x π , into its valence, x v , sea, x s , and gluon, x g components depends strongly on which data set is included in the analysis. In particular, the inclusion of LNE data, which induces a model dependence in the extraction of the pion PDFs, has a significant effect on the average momentum carried by gluons and sea quarks in the pion. Precise lattice QCD data for both quark and gluon momentum fractions has, thus, the potential to add new model independent constraints on the extraction of pion PDFs from experimental data. Finally, new data coming from planned EIC experiments, as well as from COMPASS++/Amber [9] will help to clarify the quark and gluon dynamics within the pion. On the theory side one has to resort to models [10,11] or to nonperturbative methods as provided by lattice QCD. Also from the lattice side, there is surprisingly little known for the pion. Most of the computations available so far [12][13][14][15][16][17][18][19] neglect potentially important contributions, the so-called quark disconnected contributions. While we were finalising the present work a first com-putation including disconnected contributions was put forward [20]. Also for the gluon contributions there exists only a single computation and only in the quenched approximation [21]. Thus, systematics are certainly not sufficiently controlled. More recently, there are studies using quasi- [22][23][24][25] and pseudo-distributions [26,27] as well as so-called good lattice cross sections [28][29][30] approaches to compute the x dependence of the pion PDFs directly on the lattice. These studies, however, are restricted to connected contributions only. In this letter we present the first calculation of the quark and gluon momentum fractions in the pion based on lattice QCD simulations with N f = 2 + 1 + 1 dynamical quark flavours including all required contributions. The computation is performed using one ensemble with physical values of all four quark mass parameters. This allows us to check the momentum sum rule, i.e. whether all four quark and the gluon fractions sum up to one. This result can pave the way towards a global QCD analysis including experimental data as well as lattice QCD results of the pion, which will help to sort out the discrepancy found between different experimental data sets. Lattice Computation.-Our computation is based on an ensemble [31] generated by the Extended Twisted Mass Collaboration (ETMC) using N f = 2 + 1 + 1 dynamical Wilson twisted mass clover fermions at maximal twist [32,33] and Iwasaki gauge action [34]. With this discretisation, lattice artefacts are of O(a 2 ) only [35]. The lattice volume is L 3 × T = 64 3 × 128 and the lattice spacing a = 0.08029(41) fm. For strange and charm quarks we use a mixed action approach following Ref. [36], and all quark mass parameters are tuned to assume approximately physical values [31,37]. We give further details on quark mass tuning in the appendix. For all estimates we used 745 well-separated gauge configurations. The relevant elements of the traceless Euclidean energymomentum tensor (EMT) for quark flavour q with the symmetrised covariant derivative ↔ Dµ read with κ µν = δ µ,4 δ ν,4 . Analogously for the gluon (2) For X = u, d, s, c, g, one then obtains x X from π(p)|T X µν |π(p) = 2 x X p µ p ν − δ µν with on-shell momentum p = E π = m 2 π + p 2 , p . We extract these matrix elements from ratios of Euclidean three-and two-point functions which are related to the matrix element should be maintained, otherwise finite size effects become sizable via excited states contaminations. R depends on t f − t i and t − t i only, and in the following we set t i = 0. According to Eq. (3), x can be extracted with zero pion momentum from tensor elements with µ = ν, whereas for µ = ν nonzero momentum is required. In general, one might expect the signal to be noisier with nonzero momentum, and this is indeed the case for the connectedonly contribution. However, due to the fact that for µ = ν the signal requires the subtraction of the trace of the EMT, the quark disconnected and gluon contributions are better determined from the off-diagonal elements of the EMT, see also Ref. [7]. Therefore, we determine the connected-only light contribution to x from T 44 at zero pion momentum p = 0, and all the other contributions fromT 4k with smallest nonzero momentum |p| = 2π/L, averaged over all six spatial directions. Further justification for using (off-) diagonal tensor elements for (dis-)connected diagrams is given in the appendix. For the light-quark connected part both two-and threepoint functions are constructed using stochastic timeslice sources with spin-color-site components independently and identically distributed, according to (Z 2 + i Z 2 ) / √ 2, with random Z 2 noise and eight stochastic samples per gauge field configuration. For the quark-disconnected part of any quark flavour, as well as the gluon operator part, we employ point-to-all propagators with 200 randomly distributed source points per configuration and full spin-color dilution to estimate the pion two-point function. The light-quark loop diagrams with covariant derivative insertion are determined based on the combination of low-mode deflation [38] of the Dirac operator with 200 eigenmodes, and hierarchical probing [39] with one stochastic volume source decomposed to coloring distance of eight lattice sites in each spatial direction. Spin-color dilution is also employed. The strange quark is treated with the same hierarchical probing setup, but without deflation. Analogously, for the charm-quark loops we use 12 spin-color diluted volume sources with coloring distance 4 [7]. The gluon field strength tensor in the gluon operator matrix element is computed with the clover field definition, see e.g. [40]. We apply ten levels of stout smearing [41] to the gauge links in order to sufficiently reduce ultraviolet fluctuations, see Ref. [6,7]. Errors are computed using the bootstrap method with fully correlated fits. Since lattice artefacts are of O(a 2 ) with our discretisation, we expect discretisation errors generically of the size a 2 Λ 2 QCD . For the value of the lattice spacing a used here this amounts to ∼ 2.6% with unknown coefficient. Renormalisation.-The quark flavour-nonsinglet combinations renormalise with Z qq as For the pion, x R u−d = 0 in the isospin symmetric case, as simulated here. The quark-singlet and gluon components mix under renormalisation according to with Z s qq the quark-singlet renormalisation constant. Defining δZ qq = Z s qq − Z qq , one can solve the set of Eqs. (7) for each single flavour and gluon component: Due to lattice artefacts, renormalisation factors are different forT µν with µ = ν and µ = ν. The diagonal elements of the renormalisation matrix have been determined nonperturbatively and the off-diagonal elements perturbatively in Ref. [7], see also the appendix. Since these mixing coefficients have been determined using oneloop perturbation theory, we do not have an error estimate available. In order to account for the uncertainty, we perform the computation once including the mixing, and once excluding it, and take the spread as error estimate.
Results.-We compute R X (t) in Eq. (4) for various values of t f . Solving Eq. (3) for x and inserting Eq. (5), we then extract x X (t), where X stands for l, conn (with l ≡ u + d), l, disc, s, c, and g. For large enough t f we expect x (t) to show a plateau for t−t f /2 around 0. Thus, we fit a constant symmetrically around t − t f /2 = 0 with fit range denoted as t p to our bare data for x (t) (for plots of this bare data see the appendix). In Fig. 1 we show the result of such constant fits to the light connected contribution as a function of the source sink separation t f for different values of t p . Between t f = 4.5 fm and t f = 5.14 fm we see agreement for all values of t p . For t f = 5.78 fm the results for the smallest three t p values still agree with the previous ones. However, for the larger t p values we start to see finite size effects due to T /2 < t f , also visible in the bad χ 2 /dof values. In Fig. 2  We arrive at the final result by assigning a weight w = exp − 1 2 χ 2 − 2 dof to every fit with given χ 2 value and degrees of freedom (dof). Then we take the weighted average (see also [42]) over all constant fits in the aforementioned regions of t f values. The combined statistical and systematic error is computed by repeating this procedure on all bootstrap samples with weights corresponding to the fits on the samples. Alternatively, we have also performed fits which take explicitly into account excited state contaminations again for various t f and t p values leading to consistent results, but with a different distribution of statistical and systematic errors. Using the so extracted bare values for x X (see the appendix), we are now in the position to compute the renormalised flavour nonsinglet and singlet contributions to x . We obtain for the nonsinglet ones from Eq. (6) where we recall that x u−d = 0 due to isospin symmetry in the light-quark sector. For the singlet contributions we find, using Eqs. (7) and (8), The first error represents the combined statistical and fit range uncertatinty, the second error comes from the mixing under renormalisation. The sum of all contributions amounts to x R total = 1.20(13)( −3 ), compatible with the expected value of 1 within two sigma. This is an important result because, in contrast to phenomenological analyses where the saturation of the momentum sum rule is imposed, in our work such saturation is a result of the computation. In Fig. 3 we compare to the recent phenomenological results [43] from Ref. [44] and Ref. [45]. Our work agrees within errors with these state-of-art phenomenological results. In Table I we compile all the contributions again and compare to the literature. The only other lattice result was presented in Ref. [20], where N f = 2 + 1 flavour QCD was used and results have been extrapolated to the continuum and the physical point. However, the gluon contribution has not been computed and, thus, the mixing could not be taken into account. Therefore, the comparison of the quark singlet contributions is of limited meaningfulness. The nonsinglet contribution x R u+d−2s is better suited for a comparison. While currently there is a discrepancy between the results, we notice that a full comparison should be attempted only after the inclusion of all unaccounted systematics. Finally, the results from the momentum sum rule decomposition can be used to determine how quarks and gluons contribute to the pion mass. In principle, mass can be decomposed in QCD in different ways [46][47][48][49][50][51]. Choos- Lattice (this work) Phenomenology [43] (NLO+NLL double Mellin) Phenomenology [44] FIG. 3. Comparison to Ref. [44] and Ref. [45] at 2 GeV in the MS scheme. In the two phenomenological computations the momentum sum rule is imposed.
ing the sum rule of Ref. [51], M π,q = (3M π /4) x R quarks and M π,g = (3M π /4) x R g , which amounts to 70(5) MeV and 55(12) MeV at 2 GeV in the MS scheme, respectively. Note that the gluon contribution is the same in Ji's original mass decomposition [46,47]. The remaining contribution is split among a trace anomaly term and a term proportional to the quark mass. Summary and Outlook.-In this letter we have presented results for the complete flavour decomposition of the average momentum of quarks and gluons in the pion for the first time. The computation in N f = 2+1+1 lattice QCD are performed directly with physical values of the quark mass parameters making an extrapolation to the physical point superfluous and, thus, avoiding any systematic uncertainty from such an extrapolation. However, we work at a single value of the lattice spacing only, which does not allow us to take the continuum limit. Therefore, we have to expect lattice artefacts of O(a 2 ) which we cannot account for rigorously. The renormalisation constants have been computed nonperturbatively, while the mixing coefficients were computed in perturbation theory.
We find the momentum sum rule to be fulfilled within two sigma errors, see Fig. 3. When comparing to phenomenological determinations from Refs. [44,45] we find reasonable agreement within relatively large uncertainties. Comparing to the only other lattice QCD computation [20] including quark disconnected contributions, but not the mixing and the gluon contribution, we observe a deviation well ouside uncertainties. Future plans include determining x in the pion for two more lattice spacing values directly at the physical point. Preliminary results for the flavour non-singlet components at a finer lattice spacing show agreement within errors. Moreover, work is in progress to determine the mixing coefficients nonperturbatively. This work opens the possibility to combine lattice QCD results and experimental data in a global phenomenological analysis. Acknowledgements -15 -10 -5 0 5 10 15

APPENDIX
In this section we provide additional details on the analysis and additional figures. In Fig. 4 we show x l,conn (t − t i ) for different source sink separations t f −t i . It is clearly observable that the plateau around t − t i − (t f − t i )/2 = 0 becomes longer as t f − t i increases. However, once t f − t i > T /2, deviations become larger again.
For the disconnected contributions we show as examples x g (t − t i ) and x l,disc in Fig. 5. It can be seen that the plateau region around t − t i − (t f − t i )/2 = 0 is rather independent of the shown t f − t i values, apart from increasing statistical uncertainties with increasing t f − t i values. For t f − t i = 32 (not shown) errors are so large, that all points are compatible with zero for the fermionic disconnected contributions.
In Table II Unfortunately, Z qg is not available for µ = ν. However, since Z µ =ν gq Z µ=ν gq and Z µ=ν qg itself is small compared to the other two, we assume here that Z µ=ν qg = Z µ =ν qg , also due to the fact that the divergent part in the two is identical. Another question which needs to be discussed is the justification of treating connected and disconnected light contributions separately and even using different operators. It is based on the mixed action analysis performed by the authors of Ref. [36]. Consider the bare quark contribution to x for a twisted mass quark flavor χ not present in the simulated action with twisted quark mass µ χ . The insertionT χ leads to a purely disconnected contribution. For such an insertion, at nonzero lattice spacing the corresponding x χ (µ χ ) is well defined from spectral decomposition. Lattice rotational symmetry instead guarantees invariance of x χ (µ χ ) w.r.t. the chosen tensor component up to lattice artifacts.
We use this nonunitary setup for the strange and charm quark contribution with any value of µ s , µ c . Note that lattice artefacts are then still of O(a 2 ) [36]. In the regularised theory the so obtained x χ (µ χ ) is a smooth function of the twisted valence quark mass µ χ and we can have a well-defined limit to the unitary case lim µχ→µ l x χ µχ = x l,disc , where results from different choices of tensor components only differ by lattice artefacts. Thus in the linear decomposition of x = x l,conn + x l,disc into a connected and disconnected contribution, the disconnected contribution x l,disc is fixed up to lattice artefacts. Also the complete x l,conn + x l,disc is invariant under change of tensor components by the remnant Lorentz symmetry.
In conclusion the employed decomposition into x l,conn + x l,disc in the regularised theory is unique up to lattice artefacts. Moreover, the renormalisation pattern for diagonal and off-diagonal tensor elements is identical again up to lattice artefacts, which vanish in the continuum limit.
Note that the spectral decomposition of the relevant three-point function leads to a constant plus excited states. This is why there cannot be any delicate cancellation of exponential terms like in the case for instance of the η, η correlation functions, where connected and disconnected contribution must not be analysed separately.

Quark mass tuning
The unitary light, strange and charm quark masses for the simulated lattice action are tuned such that the physical values of the pion mass, the ratio of charm and strange quark masses and the ratio of the D s meson and its decay constant are reproduced as detailed in Ref. [31]. Within our mixed action setup, the bare mass of strange and charm quark, which enter the quark energy momentum tensor in our calculation, are tuned by matching the Ω − and the Λ c baryon mass to their physical values by the procedure in Ref. [37].
The numerical values of the bare strange and charm quark parameters read a µ s = 0.0186 , a µ c = 0.2490 .