Extended Twisted Mass Collaboration

# Quark and gluon momentum fractions in the pion and in the kaon

Constantia Alexandrou  Department of Physics, University of Cyprus, Nicosia,
Cyprus  Computation-based Science and Technology Research Center, The Cyprus
Institute, Nicosia, Cyprus  Simone Bacchio  Computation-based Science and
Technology Research Center, The Cyprus Institute, Nicosia, Cyprus  Martha
Constantinou  Department of Physics, Temple University, Philadelphia, PA 19122
- 1801, USA  Joseph Delmar  Department of Physics, Temple University,
Philadelphia, PA 19122 - 1801, USA  Jacob Finkenrath  Department of Physics,
Bergische Universität Wuppertal, Gaußstraße 20, Wuppertal, 42119, Germany
Bartosz Kostrzewa  High Performance Computing and Analytics Lab, Rheinische
Friedrich-Wilhelms-Universität Bonn, Germany  Marcus Petschlies  Helmholtz-
Institut für Strahlen- und Kernphysik, University of Bonn, Germany  Bethe
Center for Theoretical Physics, University of Bonn, Germany  Luis Alberto
Rodriguez Chacon  Computation-based Science and Technology Research Center,
The Cyprus Institute, Nicosia, Cyprus  University of Ferrara, Ferrara, Italy
Gregoris Spanoudes  Department of Physics, University of Cyprus, Nicosia,
Cyprus  Fernanda Steffens  Helmholtz-Institut für Strahlen- und Kernphysik,
University of Bonn, Germany  Bethe Center for Theoretical Physics, University
of Bonn, Germany  Carsten Urbach  Helmholtz-Institut für Strahlen- und
Kernphysik, University of Bonn, Germany  Bethe Center for Theoretical Physics,
University of Bonn, Germany  Urs Wenger  Institute for Theoretical Physics,
Albert Einstein Center for Fundamental Physics, University of Bern,
Switzerland

###### Abstract

We present the full decomposition of the momentum fraction carried by quarks
and gluons in the pion and the kaon. We employ three gauge ensembles generated
with $N_{f}=2+1+1$ Wilson twisted-mass clover-improved fermions at the
physical quark masses. For both mesons we perform a continuum extrapolation
directly at the physical pion mass, which allows us to determine for the first
time the momentum decomposition at the physical point. We find that the total
momentum fraction carried by quarks is 0.532(56) and 0.618(32) and by gluons
0.388(49) and 0.408(61) in the pion and in the kaon, respectively, in the
$\overline{\mathrm{MS}}$ scheme and at the renormalization scale of 2 GeV.
Having computed both the quark and gluon contributions in the continuum limit,
we find that the momentum sum is 0.926(68) for the pion and 1.046(90) for the
kaon, verifying the momentum sum rule.

Introduction — Quantum Chromodynamics (QCD) governs the formation of hadron
bound states through the interaction of quarks and gluons. Understanding how
these constituents combine to define the observable properties of hadrons has
been the subject of extensive research spanning the last six decades. Among
the diverse spectrum of hadrons, the light pseudo-scalar mesons hold a special
place, as they encapsulate a key property of the QCD dynamics: they are
recognized as pseudo-Goldstone bosons arising from the spontaneous breaking of
chiral symmetry in QCD, characterized by their significantly smaller masses
when compared to the more abundant nucleons. In addition, they stand out for
their valence quark-antiquark composition, in contrast to the structure of
three valence quarks in the nucleons. This inner structure is mostly studied
through the use of parton distribution functions (PDFs), which measure, in the
simplest case of unpolarized hadrons, the probability of a parton to carry a
fraction $x$ of the hadron momentum. As a result, PDFs and their averages, the
so-called first moment of the PDFs, are essential quantities to pin down the
structure of any hadron.

Among the pseudo-Goldstone bosons, pions and kaons are the simplest and the
lightest. And while pions are comprised of valence up and down quarks and
antiquarks, kaons are formed by combining either an up or a down quark
(antiquark) with a strange antiquark (quark). Thus, a significant
distinguishing feature between pions and kaons lies in the mass of the strange
quark, the mass of which exceeds the sum of the masses of the two light quarks
by more than twenty fold, while the pion mass is approximately one-fourth of
the kaon mass. Moreover, in the chiral limit their masses vanish, while the
mass of the nucleon is significantly less affected. We may then say that
comprehending the pion and the kaon structures, their similarities and
differences, allows us to gain a deeper understanding of QCD itself, possibly
providing hints on the nonperturbative mechanisms behind the generation of
hadron masses. Indeed, the future Electron-Ion Collider (EIC) Aguilar  _et
al._ (2019); Abdul Khalek  _et al._ (2022); Xie  _et al._ (2022) has among its
aims to increase our knowledge of the partonic structure of pions and kaons,
and to quantify what are the quark and gluon energy contributions to the mass
of pions and kaons.

Experimental data for the pion PDFs and their moments is scarce, restricted to
limited data from several decades ago, when they were measured in pion-induced
Drell-Yan processes Conway  _et al._ (1989). For the kaon, the situation is
even worse, with only a handful of data Badier  _et al._ (1980). These early
data has been included in modern global analysis Barry  _et al._ (2018);
Novikov  _et al._ (2020); Barry  _et al._ (2021, 2022); Kotz  _et al._ (2023)
and used to determine the pion PDFs, while the kaon distributions have been
studied mostly within models Chen  _et al._ (2016); Shi  _et al._ (2018); Lan
_et al._ (2020); Bednar  _et al._ (2020); Cui  _et al._ (2020); Han  _et al._
(2021); Roberts  _et al._ (2021); Pasquini  _et al._ (2023). On the theory
side, lattice QCD can investigate the partonic structure of hadrons either
through the direct computation of the PDFs Constantinou (2021) or through
their moments, employing the method of local operators. Based on this latter
method, the body of lattice QCD studies on the pion structure remains
relatively modest in comparison to those focused on the nucleon, a situation
that is more pronounced when considering studies on the structure of kaons.
Many of the existing studies for the pion case either omit disconnected
diagrams Martinelli and Sachrajda (1987); Best  _et al._ (1997); Guagnelli
_et al._ (2005); Capitani  _et al._ (2006); Abdel-Rehim  _et al._ (2015); Oehm
_et al._ (2019); Alexandrou  _et al._ (2021a) or only incorporate a partial
portion of the disconnected contributions Löffler  _et al._ (2022). The only
exceptions are Alexandrou  _et al._ (2021b); Hackett  _et al._ (2023). In our
previous work Alexandrou  _et al._ (2021b), we performed the first
comprehensive breakdown of the pion momentum into its constituent quark and
gluon components. The authors of Hackett  _et al._ (2023) also presented a
full momentum decomposition, albeit at larger than physical pion mass and at a
single lattice spacing. While our previous work was conducted at the physical
pion mass, it was also restricted to a singe lattice spacing. In Ref. Good
_et al._ (2023), an extrapolation to the continuum is made, but using
ensembles generated with heavier-than-physical pion masses. Moreover, their
result is restricted to the gluons only. In the case of the kaon, there is one
previous study Alexandrou  _et al._ (2021a), where only connected
contributions are considered and again a heavier than physical pion mass is
used. Going beyond the use of local operators, the last few years witnessed
the development of a variety of new methods Constantinou (2021), based on
nonlocal operators, which allow us to access the whole Bjorken $x$ dependence
of PDFs. In particular, the $x$ dependence of PDFs for pions and kaons has
been computed within the quasi-PDF Zhang  _et al._ (2019); Izubuchi  _et al._
(2019); Lin  _et al._ (2021); Gao  _et al._ (2020) and the pseudo-PDF Joó  _et
al._ (2019); Salas-Chavira  _et al._ (2022) approaches, as well as in the so
called “good lattice cross sections” method involving hadronic matrix elements
computed on the lattice that are constructed to be time independent Sufian
_et al._ (2019, 2020). Nevertheless, these works are restricted to the valence
quark distributions, with the exception of Salas-Chavira  _et al._ (2022)
where only the gluon distribution is computed using ensembles with heavier
than physical pion masses. It is clear from this discussion that the full
momentum decomposition of the pion and of the kaon is needed, both at the
physical pion mass and with an extrapolation to the continuum. Therefore, in
this work we perform the first momentum decomposition of both pions and kaons
using a lattice QCD computation at the physical point.

Lattice Setup — The parameters of the gauge field ensembles of the Extended
Twisted Mass Collaboration (ETMC) used in this work are listed in Table 1,
with the corresponding lattice spacing, lattice size and two most relevant
hadron masses for charged pion and kaon. All ensembles have $N_{f}=2+1+1$
quarks, with masses of up, down, strange and charm quark tuned to reproduce
the physical pion mass, as well as the physical value of the ratios
$M_{D_{s}}/f_{D_{s}}$ and $m_{c}/m_{s}$, see Ref. Alexandrou  _et al._ (2018)
for details.

Table 1: ETMC ensembles analyzed. $a$ is the lattice spacing and $L\,(T=2L)$ the lattice spatial (temporal) extent in fm, and $M_{\pi^{\pm}}$ and $M_{K^{\pm}}$ the charged pion and kaon mass, respectively.  Ensemble | a [fm] | L [fm] |  $M_{\pi^{\pm}}$ [MeV] |  $M_{K^{\pm}}$ [MeV]
---|---|---|---|---
cB211.072.64 (B) | $0.0796(1)$ | $5.09$ | $140.40(22)$ | $498.41(11)$
cC211.060.80 (C) | $0.0682(1)$ | $5.46$ | $136.05(30)$ | $495.27(14)$
cD211.054.96 (D) | $0.0569(1)$ | $5.46$ | $141.01(22)$ | $494.77(11)$

We extract the momentum fractions carried by the partons from the matrix
elements of the energy-momentum tensor (EMT),

$\langle
h(\mathbf{p})\,|\,\bar{T}^{X}_{\mu\nu}\,|\,h(\mathbf{p})\rangle=2\langle
x\rangle_{X}\left(p_{\mu}p_{\nu}-\delta_{\mu\nu}\frac{p^{2}}{4}\right)$ (1)

where $h(\mathbf{p})$ denotes a hadron state with mass $m_{h}$, and
$p=(E_{h}(\mathbf{p})\,,\,\,\mathbf{p})$ the on-shell four-momentum, with
$E_{h}(\mathbf{p})=\sqrt{m_{h}^{2}+\mathbf{p}^{2}}$. The index $X=q,\,g$
denotes the contributions from individual quark flavors ($q=u,\,d,\,s,\,c$)
and gluons to the hadron momentum. We compute the matrix elements appearing on
the left-hand side of Eq. (1) through the ratio of three- and two-point
functions

$R_{\mu\nu}^{X}(t_{\mathrm{ins}},t_{\mathrm{s}};\,\mathbf{p})\ =\
\frac{\langle
h(t_{\mathrm{s}},\mathbf{p})\,\bar{T}^{X}_{\mu\nu}(t_{\mathrm{ins}})\,h(0,\mathbf{p})\rangle}{\langle
h(t_{\mathrm{s}},\mathbf{p})\ h(0,\mathbf{p})\rangle}.$ (2)

If $t_{\mathrm{s}}$ and $t_{\mathrm{s}}-t_{\mathrm{ins}}$ are large enough,
then asymptotically up to contamination by excited states above the ground
state $|h(\mathbf{p})\rangle$

$R_{\mu\nu}^{X}(t,t_{f},t_{i};\mathbf{p})\to\frac{1}{2E_{h}}\,\frac{\langle
h(\mathbf{p})\,|\,\bar{T}^{X}_{\mu\nu}\,|\,h(\mathbf{p})\rangle}{1+\exp(-E_{h}(T-2\,t_{\mathrm{s}}))}\,.$
(3)

The quark part of the EMT is given by

$\bar{T}_{\mu\nu}^{q}\ =\
-\frac{(i)^{\kappa_{\mu\nu}}}{4}\,\bar{q}\,\left(\gamma_{\mu}\stackrel{{\scriptstyle\leftrightarrow}}{{D}}_{\nu}+\gamma_{\nu}\stackrel{{\scriptstyle\leftrightarrow}}{{D}}_{\mu}-\delta_{\mu\nu}\,\frac{1}{2}\,\gamma_{\rho}\,\stackrel{{\scriptstyle\leftrightarrow}}{{D}}_{\rho}\right)q\,,$
(4)

with $\kappa_{\mu\nu}=\delta_{\mu,4}\,\delta_{\nu,4}$, and
$\stackrel{{\scriptstyle\leftrightarrow}}{{D}}_{\mu}$ the symmetrized
covariant derivative. The gluon EMT is given through the gluon field-strength
tensor $F_{\mu\nu}$ by

$\bar{T}^{g}_{\mu\nu}\ =\
(i)^{\kappa_{\mu\nu}}\,\left(F_{\mu\rho}^{~{}}F_{\nu\rho}+F_{\nu\rho}^{~{}}F_{\mu\rho}-\delta_{\mu\nu}\,\frac{1}{2}\,F_{\rho\sigma}{~{}}F_{\rho\sigma}\right)\,.$
(5)

We note that both expressions have been written for a Euclidean space-time.
While quark-connected and -disconnected diagrams contribute to the 3-point
functions with quark EMT insertion, the gluonic EMT insertion produces only
purely disconnected diagrams (cf. Fig. 10 in the Supplementary Material I).

Bare momentum fractions — Equation (3), combined with Eq. (1), allows us to
determine $\langle x\rangle$ from different components of the EMT. We compute
the connected contributions using the tensor element $\bar{T}_{44}$ with
hadron momentum $\mathbf{p}=0$, while for the disconnected contributions we
employ the combination of components $\bar{T}_{4k}$, $k=1,2,3$, with nonzero
momentum $|\mathbf{p}|=2\pi/L$. We thereby optimize the signal-to-noise ratio
for each type of contribution: the temporal-spatial EMT elements
$\bar{T}_{4k}$ do not require additive renormalization by trace subtraction
and with minimal nonzero hadron momentum allow for most precise projection on
the hadron ground state, necessary for the weak signal in the disconnected
diagrams. On the other hand, statistical noise from trace subtraction is
negligible for the connected diagrams from $\bar{T}_{44}$, while zero hadron
momentum renders the ground state projection most precise. By Eq. (1) and
Lorentz-covariance, the matrix elements from different tensor components
become equivalent in the continuum limit. As an illustration of the data, we
show the ratios for connected and disconnected contributions for the light
quarks in the kaon in Fig. 1 for the ensemble C in Table 1.

Figure 1: Bare ratios of the light-quark momentum fraction in the kaon
dependence on the sink-source time separation $t_{\mathrm{s}}$ for the C
ensemble. In the upper panel, we show the ratio for the connected contribution
and in the lower panel, the ratio for the disconnected contribution. We show
results for several values of $t_{\mathrm{s}}$, from $t_{\mathrm{s}}/a=30$ to
$t_{\mathrm{s}}/a=80$. The red bands show the result after model averaging of
fits to a constant when varying the ranges of $t_{\mathrm{s}}$ used.

The ratios were computed for several source-sink separations $t_{\mathrm{s}}$,
and we see in Fig. 1 the formation of a plateau for large values of
$t_{\mathrm{s}}$. We fitted the plateau using different combinations of
$t_{\mathrm{s}}$ and of the insertion time $t_{\mathrm{ins}}$. Model averaging
Jay and Neil (2021) is used on the results from the plateau fit in order to
obtain $\langle x\rangle$. The same procedure is employed for the ratios
involving the gluon EMT insertion. To optimize the signal-to-noise ratio we
apply four-dimensional stout smearing Morningstar and Peardon (2004) to the
gauge field in the lattice gluon field-strength tensor when constructing
$\bar{T}^{g}_{4k}$. In Fig. 2 we show the ratio for several values of
$t_{\mathrm{s}}$ for $N_{\mathrm{stout}}=10$ smearing steps, for both pion and
kaon case.

We probe the systematic dependence on the level of stout smearing by comparing
the average momentum fraction for the range of stout-smearing steps $5\leq
N_{\mathrm{stout}}\leq 10$, weighted by the renormalization factor $Z_{gg}$,
which is determined for the corresponding values of $N_{\mathrm{stout}}$. As
we find agreement among those individual results, we then perform a
simultaneous fit to the $Z_{gg}$-weighted ratios from this set of stout-
smearing steps.

Figure 2: Examples of the bare ratios for the gluon contribution, for six
values of the time sink-source separation $t_{s}$ using 10 stout-smearing
steps. The top panel shows the pion and the bottom panel the kaon case. The
horizontal bands show the model-averaged fit results.

All the plateau fits, as well as the results after model averaging, for each
ensemble and flavor, are shown in the Supplemental Material I.

Renormalization and continuum limit — We proceed with the renormalization
procedure for the flavor-singlet and flavor-nonsinglet combinations. The
renormalized results are given in the $\overline{\mathrm{MS}}$ scheme at scale
2 GeV. The flavor-nonsinglet combinations are renormalized with the
nonperturbatively determined renormalization factor $Z_{qq}$ as

$\begin{split}\langle x\rangle_{u+d-2s,\mathrm{R}}&=Z_{qq}(\langle
x\rangle_{u}+\langle x\rangle_{d}-2\langle x\rangle_{s})\,,\\\ \langle
x\rangle_{u+d+s-3c,\mathrm{R}}&=Z_{qq}(\langle x\rangle_{u}+\langle
x\rangle_{d}+\langle x\rangle_{s}-3\langle x\rangle_{c})\,.\\\ \end{split}$
(6)

Compared to our previous analysis Alexandrou  _et al._ (2021b), we now use
$\mathcal{O}(\alpha_{s}^{4})$ corrections to the anomalous dimension, thus
improving the scale evolution of the data to the common $2$ GeV scale used
throughout this work. The quark-singlet and gluon components mix under
renormalization according to

$\begin{pmatrix}\langle x\rangle_{q,\mathrm{R}}\\\ \langle
x\rangle_{g,\mathrm{R}}\end{pmatrix}=\begin{pmatrix}Z_{qq}^{s}&Z_{qg}\\\
Z_{gq}&Z_{gg}\\\ \end{pmatrix}\begin{pmatrix}\langle x\rangle_{q}\\\ \langle
x\rangle_{g}\end{pmatrix}$ (7)

with $Z_{qq}^{s}$ the quark-singlet renormalization constant and the total
bare (renormalized) quark contribution is denoted by

$\displaystyle\langle x\rangle_{q,(\mathrm{R})}$
$\displaystyle=\sum\limits_{f=u,d,s,c}\,\langle x\rangle_{f,(\mathrm{R})}\,.$
(8)

Defining $\delta Z_{qq}=Z_{qq}^{s}-Z_{qq}$, one can solve Eq. (7) for each
single flavor and the gluon component:

$\begin{split}\langle x\rangle_{f,\mathrm{R}}&=Z_{qq}\,\langle
x\rangle_{f}+\frac{\delta Z_{qq}}{N_{f}}\,\sum_{f^{\prime}}\,\langle
x\rangle_{f^{\prime}}+\frac{Z_{qg}}{N_{f}}\,\langle x\rangle_{g}\,,\\\ \langle
x\rangle_{g,\mathrm{R}}&=Z_{gg}\,\langle
x\rangle_{g}+Z_{gq}\,\sum_{f^{\prime}}\langle x\rangle_{f^{\prime}}\,.\\\
\end{split}$ (9)

The calculation of $Z_{gq}$ and $Z_{qg}$ have also been improved, with the
details presented in the Supplemental Material.

For the twisted-mass lattice action at maximal twist with automatic
$\mathcal{O}(a)$ improvement we expect leading discretization artefacts at
second order in the lattice spacing $a$. We, thus, take the continuum limit
using a constant or a linear fit in $a^{2}$, for both the renormalized quark
and gluon contributions, for the pion and the kaon cases. Our final value is
the result of the model average among these fits. In Fig. (3), we show the
linear extrapolation of the total renormalized quark contribution $\langle
x\rangle_{q,\mathrm{R}}$, the gluon contribution $\langle
x\rangle_{g,\mathrm{R}}$, and in addition their sum. We note that the total
quark contribution to $\langle x\rangle^{K}$ tends to be higher than its
analog in the pion, as also observed in Ref Bednar  _et al._ (2020). The gluon
contribution to the pion and to the kaon momenta, however, are compatible
within errors, and the momentum sum rule

$\displaystyle\langle x\rangle_{q,\mathrm{R}}+\langle
x\rangle_{g,\mathrm{R}}=1$ (10)

is satisfied, in both cases, within a standard deviation.

Figure 3: Continnum limit extrapolation for the pion (left panel) and the kaon
(right panel). We present our results for the total quark and gluon
contributions, as well as the momentum sum rule. The blue filled circles are
results for ensemble B Alexandrou  _et al._ (2021b), the orange upwards
triangles for ensemble C and the green downwards triangles for ensemble D. The
open symbol is the result after model averaging between constant and linear
fit. Further details can be found in the Supplemental Material.

The full list of the renormalized results in the continuum limit, for both
pion and kaon, is presented in Table 2 at the scale of $\mu=2$ GeV in the
$\overline{\mathrm{MS}}$ scheme.

In Fig. 4, we show a comparison of our results with those from other groups.
We include recent results from phenomenological analyses Novikov  _et al._
(2020); Barry  _et al._ (2021), from the Dyson Schwinger equations (DSE)
Bednar  _et al._ (2020), from a calculation using the light-front wave
function (LFWF) approach Cui  _et al._ (2020) and from lattice QCD
computations Löffler  _et al._ (2022); Good  _et al._ (2023), limiting
ourselves to those which are extrapolated to the physical pion mass and/or the
continuum limit. The lattice QCD results shown are computed using ensembles
simulated with pion masses larger than physical and then extrapolated to the
physical pion mass. Furthermore, they both considered only a partial sum of
disconnected contributions, namely only the quark disconnected is considered
in Ref. Löffler  _et al._ (2022) but not the gluon contribution, whereas only
the gluon is considered in Ref. Good  _et al._ (2023). For the cases where all
contributions are included, the data is restricted to only one lattice spacing
Alexandrou  _et al._ (2021b) at physical pion mass or using larger than
physical pion masses Hackett  _et al._ (2023); Alexandrou  _et al._ (2021a).
Our result for $\langle x\rangle_{g,R}^{\pi}$ is in agreement with the one
from Ref. Good  _et al._ (2023) as is our result for $\langle
x\rangle_{q,R}^{\pi}$ with RQCD Löffler  _et al._ (2022) that carries however
a large error. Compared to our previous work Alexandrou  _et al._ (2021b)
using only the B ensemble, taking the continuum extrapolation increases the
error, as can be seen in Fig. 3. For the kaon, the only available lattice QCD
data Alexandrou  _et al._ (2021a) is for one lattice spacing using an ensemble
with a heavier-than-physical pion mass and thus one cannot directly compare to
the present work.

Table 2: Compilation of results for the pion and for the kaon in the continuum limit. All quantities are presented at the scale $2\ \mathrm{GeV}$ in the $\overline{\mathrm{MS}}$ scheme. | $\pi$ | K
---|---|---
$\langle x\rangle_{l,\mathrm{R}}$ | $0.448(34)$ | $0.260(09)$
$\langle x\rangle_{s,\mathrm{R}}$ | $0.043(15)$ | $0.333(11)$
$\langle x\rangle_{c,\mathrm{R}}$ | $0.019(17)$ | $0.024(17)$
$\langle x\rangle_{g,\mathrm{R}}$ | $0.388(49)$ | $0.408(61)$
$\langle x\rangle_{q,\mathrm{R}}$ | $0.532(56)$ | $0.618(32)$
$\langle x\rangle_{u+d-2s,\mathrm{R}}$ | $0.382(17)$ | $-0.409(16)$
$\langle x\rangle_{u+d+s-3c,\mathrm{R}}$ | $0.445(48)$ | $0.487(39)$
Figure 4: Comparison of the results of this work, with other available data,
both from phenomenology and from lattice QCD. All results are given in the
$\overline{\mathrm{MS}}$ scheme at the scale of $\mu=2\ \mathrm{GeV}$. The
upper panels show the results for the pion gluon (left) and quark (right)
momentum fractions, $\langle x\rangle_{\textrm{g,R}}^{\pi}$ and $\langle
x\rangle_{\textrm{q,R}}^{\pi}$, respectively. The lower panels show the
corresponding results for the kaon. The red filled squares show the results of
this work with the red band the associated error band. Recent results from
phenomenological analyses of PDFs data are given by open symbols: left green
triangle (Novikov et al. Novikov  _et al._ (2020)) and right orange triangle
(JAM Collaboration Barry  _et al._ (2021)). The result while the result based
on the LFWF is represented by the down blue triangle (Cui et al. Cui  _et al._
(2020)), while using the DSE Bednar  _et al._ (2020) by the black cross, where
no error is provided. Recent lattice QCD results extrapolated to the continuum
limit are given by the brown filled circle (RQCD Löffler  _et al._ (2022)) and
the gray pentagon (MSULat Good  _et al._ (2023)).

Conclusions and outlook — We present in this work the first complete momentum
decomposition for both the pion and the kaon in terms of their quark and gluon
consituents, performed within lattice QCD at the physical point. Three
ensembles with $N_{f}=2+1+1$ quark flavors with their masses tuned to
reproduce the physical light, strange and charm quark masses are analyzed. We
use a model average for the continuum extrapolation performed using both a
constant and a linear fit in the squared lattice spacing, $a^{2}$. Our results
for $\langle x\rangle_{g}^{\pi,K}$ indicate a similar momentum fraction
carried by gluons in the kaon and the pion. However, they have larger errors
than the corresponding quark momentum fraction $\langle x\rangle_{q}^{\pi,K}$,
which tends to be smaller in the pion. This indicates a possible larger
momentum fraction carried by gluons in the pion as compared to the kaon. Thus,
among possible improvements, in future work both the statistical and
systematic errors, in particular for $\langle x\rangle_{g}$, will be reduced.
This can be accomplished by increasing the number of gauge configurations and
source positions used for the two-point functions, Eq. (3), improving the
statistical accuracy in the determination of the gluon momentum fraction. In
addition, we plan to analyze an additional ensemble closer to the continuum
limit in order to help a better error control when performing the continuum
extrapolation. However, the current work presents a remarkable achievement as
we now have a fully consistent theoretical calculation that allows us to
directly compare the pion and the kaon at the level of their quark and gluon
structure. In particular, the momentum sum rule can be tested by computing all
components from first principles. We stress that prior to the current work,
there was no such decomposition into quark and gluon parts available of
$\langle x\rangle^{K}$ using a first-principles calculation. Furthermore, the
current work demonstrates clearly that the momentum fraction carried by gluons
and sea quarks (the disconnected contributions) are important components of
the two lightest pseudo-Goldstone bosons, and added together with the valence
contributions verify the momentum sum rule for both cases.

###### Acknowledgements.

G.S. acknowledges financial support from the European Regional Development
Fund and the Republic of Cyprus through the Cyprus Research and Innovation
Foundation under contract number EXCELLENCE/0421/0195. This project is partly
funded by the European Union’s Horizon 2020 Research and Innovation Programme
ENGAGE under the Marie Sklodowska-Curie COFUND scheme with grant agreement No.
101034267. C.A. acknowledges partial support from the European Regional
Development Fund and the Republic of Cyprus through the Cyprus Research and
Innovation Foundation under contract number EXCELLENCE/0421/0043 and the
European Joint Doctorate project AQTIVATE funded by the European Commission
under the Marie Sklodowska-Curie Doctoral Networks action and Grant Agreement
No 101072344. This work is supported by the Swiss National Science Foundation
(SNSF) through grant No. 200021_175761, 200021_208222, and 200020_200424, as
well as by the DFG and the NFSC as part of the Sino-German Collaborative
Research Center CRC110 “Symmetries and the emergence of structure”. J. D. and
M. C. acknowledge financial support from the U.S. Department of Energy, Office
of Nuclear Physics, Early Career Award under Grant No. DE-SC0020405. We
gratefully acknowledge computing time granted on Piz Daint at Centro Svizzero
di Calcolo Scientifico (CSCS) via the projects s849, s982, s1045, s1133 and
s1197 and the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for
funding this project by providing computing time on the GCS Supercomputers
SuperMUC-NG at Leibniz Supercomputing Centre and JUWELS Jülich Supercomputing
Centre (2021) at Juelich Supercomputing Centre. The authors acknowledge the
Texas Advanced Computing Center (TACC) at The University of Texas at Austin
for providing HPC resources (Project ID PHY21001). J.F. is supported by the
DFG research unit FOR5269 ”Future methods for studying confined gluons in
QCD”. The authors gratefully acknowledge PRACE for awarding access to HAWK at
HLRS within the project with Id Acid 4886, and the Swiss National
Supercomputing Centre (CSCS) and the EuroHPC Joint Undertaking for awarding
this project access to the LUMI supercomputer, owned by the EuroHPC Joint
Undertaking, hosted by CSC (Finland) and the LUMI consortium through the
Chronos programme under project IDs CH17-CSCS-CYP and CH21-CSCS-UNIBE as well
as the EuroHPC Regular Access Mode under project ID EHPC-REG-2021R0095. The
open source software packages tmLQCD Jansen and Urbach (2009); Abdel-Rehim
_et al._ (2014); Deuzeman  _et al._ (2013), Lemon Deuzeman  _et al._ (2012),
QUDA Clark  _et al._ (2010); Babich  _et al._ (2011); Clark  _et al._ (2016),
R R Core Team (2019), cvc Petschlies et al. (2024) and plegma have been used.
Finally, We thank the authors of Cui  _et al._ (2020) and Bednar  _et al._
(2020) for kindly sending us their results evolved to the scale used in our
computation, $2\ \mathrm{GeV}$.

## References

  * Aguilar  _et al._ (2019) Arlene 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] .
  * Abdul Khalek  _et al._ (2022) R. Abdul Khalek _et al._ , “Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report,” Nucl. Phys. A 1026, 122447 (2022), arXiv:2103.05419 [physics.ins-det] .
  * Xie  _et al._ (2022) Gang Xie, Chengdong Han, Rong Wang, and Xurong Chen, “Tackling the kaon structure function at EicC *,” Chin. Phys. C 46, 064107 (2022), arXiv:2109.08483 [hep-ph] .
  * 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–122 (1989).
  * Badier  _et al._ (1980) J. Badier _et al._ (Saclay-CERN-College de France-Ecole Poly-Orsay), “Measurement of the $K^{-}/\pi^{-}$ Structure Function Ratio Using the Drell-Yan Process,” Phys. Lett. B 93, 354–356 (1980).
  * Barry  _et al._ (2018) P. C. Barry, N. Sato, W. Melnitchouk, and Chueng-Ryong Ji, “First Monte Carlo Global QCD Analysis of Pion Parton Distributions,” Phys. Rev. Lett. 121, 152001 (2018), arXiv:1804.01965 [hep-ph] .
  * Novikov  _et al._ (2020) Ivan Novikov _et al._ , “Parton Distribution Functions of the Charged Pion Within The xFitter Framework,” Phys. Rev. D 102, 014040 (2020), arXiv:2002.02902 [hep-ph] .
  * Barry  _et al._ (2021) P. C. Barry, Chueng-Ryong Ji, N. Sato, and W. Melnitchouk (Jefferson Lab Angular Momentum (JAM)), “Global QCD Analysis of Pion Parton Distributions with Threshold Resummation,” Phys. Rev. Lett. 127, 232001 (2021), arXiv:2108.05822 [hep-ph] .
  * Barry  _et al._ (2022) P. C. Barry _et al._ (Jefferson Lab Angular Momentum (JAM), HadStruc), “Complementarity of experimental and lattice QCD data on pion parton distributions,” Phys. Rev. D 105, 114051 (2022), arXiv:2204.00543 [hep-ph] .
  * Kotz  _et al._ (2023) Lucas Kotz, Aurore Courtoy, Pavel Nadolsky, Fredrick Olness, and Maximiliano Ponce-Chavez, “An analysis of parton distributions in a pion with Bézier parametrizations,” (2023), arXiv:2311.08447 [hep-ph] .
  * Chen  _et al._ (2016) Chen Chen, Lei Chang, Craig D. Roberts, Shaolong Wan, and Hong-Shi Zong, “Valence-quark distribution functions in the kaon and pion,” Phys. Rev. D 93, 074021 (2016), arXiv:1602.01502 [nucl-th] .
  * Shi  _et al._ (2018) Chao Shi, Cédric Mezrag, and Hong-shi Zong, “Pion and kaon valence quark distribution functions from Dyson-Schwinger equations,” Phys. Rev. D 98, 054029 (2018), arXiv:1806.10232 [nucl-th] .
  * Lan  _et al._ (2020) Jiangshan Lan, Chandan Mondal, Shaoyang Jia, Xingbo Zhao, and James P. Vary, “Pion and kaon parton distribution functions from basis light front quantization and QCD evolution,” Phys. Rev. D 101, 034024 (2020), arXiv:1907.01509 [nucl-th] .
  * Bednar  _et al._ (2020) Kyle D. Bednar, Ian C. Cloët, and Peter C. Tandy, “Distinguishing Quarks and Gluons in Pion and Kaon Parton Distribution Functions,” Phys. Rev. Lett. 124, 042002 (2020), arXiv:1811.12310 [nucl-th] .
  * Cui  _et al._ (2020) Zhu-Fang Cui, Minghui Ding, Fei Gao, Khépani Raya, Daniele Binosi, Lei Chang, Craig D Roberts, Jose Rodríguez-Quintero, and Sebastian M Schmidt, “Kaon and pion parton distributions,” Eur. Phys. J. C 80, 1064 (2020).
  * Han  _et al._ (2021) Chengdong Han, Gang Xie, Rong Wang, and Xurong Chen, “An Analysis of Parton Distribution Functions of the Pion and the Kaon with the Maximum Entropy Input,” Eur. Phys. J. C 81, 302 (2021), arXiv:2010.14284 [hep-ph] .
  * Roberts  _et al._ (2021) Craig D. Roberts, David G. Richards, Tanja Horn, and Lei Chang, “Insights into the emergence of mass from studies of pion and kaon structure,” Prog. Part. Nucl. Phys. 120, 103883 (2021), arXiv:2102.01765 [hep-ph] .
  * Pasquini  _et al._ (2023) Barbara Pasquini, Simone Rodini, and Simone Venturini (MAP (Multi-dimensional Analyses of Partonic distributions)), “Valence quark, sea, and gluon content of the pion from the parton distribution functions and the electromagnetic form factor,” Phys. Rev. D 107, 114023 (2023), arXiv:2303.01789 [hep-ph] .
  * Constantinou (2021) Martha Constantinou, “The x-dependence of hadronic parton distributions: A review on the progress of lattice QCD,” Eur. Phys. J. A 57, 77 (2021), arXiv:2010.02445 [hep-lat] .
  * Martinelli and Sachrajda (1987) G. Martinelli and Christopher T. Sachrajda, “Pion Structure Functions From Lattice QCD,” Phys. Lett. B 196, 184–190 (1987).
  * Best  _et al._ (1997) C. Best, M. Göckeler, R. Horsley, Ernst-Michael Ilgenfritz, H. Perlt, Paul E. L. Rakow, A. Schäfer, G. Schierholz, A. Schiller, and S. Schramm, “Pion and rho structure functions from lattice QCD,” Phys. Rev. D 56, 2743–2754 (1997), arXiv:hep-lat/9703014 .
  * Guagnelli  _et al._ (2005) M. Guagnelli, K. Jansen, F. Palombi, R. Petronzio, A. Shindler, and I. Wetzorke (Zeuthen-Rome (ZeRo)), “Non-perturbative pion matrix element of a twist-2 operator from the lattice,” Eur. Phys. J. C 40, 69–80 (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–526 (2006), arXiv:hep-lat/0511013 .
  * 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 (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._ (2021a) Constantia Alexandrou, Simone Bacchio, Ian Cloet, Martha Constantinou, Kyriakos Hadjiyiannakou, Giannis Koutsou, and Colin Lauer (ETM), “Mellin moments $\langle x\rangle$ and $\langle x^{2}\rangle$ for the pion and kaon from lattice QCD,” Phys. Rev. D 103, 014508 (2021a), arXiv:2010.03495 [hep-lat] .
  * Löffler  _et al._ (2022) Marius Löffler, Philipp Wein, Thomas Wurm, Simon Weishäupl, Daniel Jenkins, Rudolf Rödl, Andreas Schäfer, and Lisa Walter (RQCD), “Mellin moments of spin dependent and independent PDFs of the pion and rho meson,” Phys. Rev. D 105, 014505 (2022), arXiv:2108.07544 [hep-lat] .
  * Alexandrou  _et al._ (2021b) Constantia Alexandrou _et al._ (Extended Twisted Mass), “Quark and Gluon Momentum Fractions in the Pion from Nf=2+1+1 Lattice QCD,” Phys. Rev. Lett. 127, 252001 (2021b), arXiv:2109.10692 [hep-lat] .
  * Hackett  _et al._ (2023) Daniel C. Hackett, Patrick R. Oare, Dimitra A. Pefkou, and Phiala E. Shanahan, “Gravitational form factors of the pion from lattice QCD,” (2023), arXiv:2307.11707 [hep-lat] .
  * Good  _et al._ (2023) William Good, Kinza Hasan, Allison Chevis, and Huey-Wen Lin, “The Gluon Moment and Parton Distribution Function of the Pion from $N_{f}=2+1+1$ Lattice QCD,” (2023), arXiv:2310.12034 [hep-lat] .
  * Zhang  _et al._ (2019) Jian-Hui Zhang, Jiunn-Wei Chen, Luchang Jin, Huey-Wen Lin, Andreas Schäfer, and Yong 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) Taku Izubuchi, Luchang Jin, Christos Kallidonis, Nikhil Karthik, Swagato Mukherjee, Peter Petreczky, Charles Shugert, and Sergey Syritsyn, “Valence parton distribution function of pion from fine lattice,” Phys. Rev. D 100, 034516 (2019), arXiv:1905.06349 [hep-lat] .
  * Lin  _et al._ (2021) Huey-Wen Lin, Jiunn-Wei Chen, Zhouyou Fan, Jian-Hui Zhang, and Rui Zhang, “Valence-Quark Distribution of the Kaon and Pion from Lattice QCD,” Phys. Rev. D 103, 014516 (2021), arXiv:2003.14128 [hep-lat] .
  * Gao  _et al._ (2020) Xiang Gao, Luchang Jin, Christos Kallidonis, Nikhil Karthik, Swagato Mukherjee, Peter Petreczky, Charles Shugert, Sergey Syritsyn, and Yong Zhao, “Valence parton distribution of the pion from lattice QCD: Approaching the continuum limit,” Phys. Rev. D 102, 094513 (2020), arXiv:2007.06590 [hep-lat] .
  * Joó  _et al._ (2019) Bálint Joó, Joseph Karpie, Kostas Orginos, Anatoly V. Radyushkin, David G. Richards, Raza Sabbir Sufian, and Savvas Zafeiropoulos, “Pion valence structure from Ioffe-time parton pseudodistribution functions,” Phys. Rev. D 100, 114512 (2019), arXiv:1909.08517 [hep-lat] .
  * Salas-Chavira  _et al._ (2022) Alejandro Salas-Chavira, Zhouyou Fan, and Huey-Wen Lin, “First glimpse into the kaon gluon parton distribution using lattice QCD,” Phys. Rev. D 106, 094510 (2022), arXiv:2112.03124 [hep-lat] .
  * Sufian  _et al._ (2019) Raza Sabbir Sufian, Joseph Karpie, Colin Egerer, Kostas Orginos, Jian-Wei Qiu, and David 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) Raza Sabbir Sufian, Colin Egerer, Joseph Karpie, Robert G. Edwards, Bálint Joó, Yan-Qing Ma, Kostas Orginos, Jian-Wei Qiu, and David G. Richards, “Pion Valence Quark Distribution from Current-Current Correlation in Lattice QCD,” Phys. Rev. D 102, 054508 (2020), arXiv:2001.04960 [hep-lat] .
  * Alexandrou  _et al._ (2018) Constantia Alexandrou _et al._ , “Simulating twisted mass fermions at physical light, strange and charm quark masses,” Phys. Rev. D 98, 054518 (2018), arXiv:1807.00495 [hep-lat] .
  * Jay and Neil (2021) William I. Jay and Ethan T. Neil, “Bayesian model averaging for analysis of lattice field theory results,” Phys. Rev. D 103, 114502 (2021), arXiv:2008.01069 [stat.ME] .
  * Morningstar and Peardon (2004) Colin Morningstar and Mike J. Peardon, “Analytic smearing of SU(3) link variables in lattice QCD,” Phys. Rev. D 69, 054501 (2004), arXiv:hep-lat/0311018 .
  * Jülich Supercomputing Centre (2021) Jülich Supercomputing Centre, “JUWELS Cluster and Booster: Exascale Pathfinder with Modular Supercomputing Architecture at Juelich Supercomputing Centre,” Journal of large-scale research facilities 7 (2021), 10.17815/jlsrf-7-183.
  * Jansen and Urbach (2009) K. Jansen and C. Urbach, “tmLQCD: A Program suite to simulate Wilson Twisted mass Lattice QCD,” Comput.Phys.Commun. 180, 2717–2738 (2009), arXiv:0905.3331 [hep-lat] .
  * Abdel-Rehim  _et al._ (2014) Abdou Abdel-Rehim, Florian Burger, Alber Deuzeman, Karl Jansen, Bartosz Kostrzewa, Luigi Scorzato, and Carsten Urbach, “Recent developments in the tmLQCD software suite,” PoS LATTICE2013, 414 (2014), arXiv:1311.5495 [hep-lat] .
  * Deuzeman  _et al._ (2013) A. Deuzeman, K. Jansen, B. Kostrzewa, and C. Urbach, “Experiences with OpenMP in tmLQCD,” PoS LATTICE2013, 416 (2013), arXiv:1311.4521 [hep-lat] .
  * Deuzeman  _et al._ (2012) Albert Deuzeman, Siebren Reker, and Carsten Urbach (ETM), “Lemon: an MPI parallel I/O library for data encapsulation using LIME,” Comput. Phys. Commun. 183, 1321–1335 (2012), arXiv:1106.4177 [hep-lat] .
  * Clark  _et al._ (2010) M. A. Clark, R. Babich, K. Barros, R. C. Brower, and C. Rebbi, “Solving Lattice QCD systems of equations using mixed precision solvers on GPUs,” Comput. Phys. Commun. 181, 1517–1528 (2010), arXiv:0911.3191 [hep-lat] .
  * Babich  _et al._ (2011) R. Babich, M. A. Clark, B. Joo, G. Shi, R. C. Brower, and S. Gottlieb, “Scaling Lattice QCD beyond 100 GPUs,” in  _SC11 International Conference for High Performance Computing, Networking, Storage and Analysis Seattle, Washington, November 12-18, 2011_ (2011) arXiv:1109.2935 [hep-lat] .
  * Clark  _et al._ (2016) M. A. Clark, Bálint Joó, Alexei Strelchenko, Michael Cheng, Arjun Gambhir, and Richard Brower, “Accelerating Lattice QCD Multigrid on GPUs Using Fine-Grained Parallelization,” (2016), arXiv:1612.07873 [hep-lat] .
  * R Core Team (2019) R Core Team,  _R: A Language and Environment for Statistical Computing_, R Foundation for Statistical Computing, Vienna, Austria (2019).
  * Petschlies et al. (2024) M. Petschlies et al., “cvc package,” https://github.com/marcuspetschlies/cvc (2024).
  * Alexandrou  _et al._ (2017) Constantia Alexandrou, Martha Constantinou, and Haralambos Panagopoulos (ETM), “Renormalization functions for Nf=2 and Nf=4 twisted mass fermions,” Phys. Rev. D 95, 034505 (2017), arXiv:1509.00213 [hep-lat] .
  * Caracciolo  _et al._ (1992) Sergio Caracciolo, Pietro Menotti, and Andrea Pelissetto, “One loop analytic computation of the energy momentum tensor for lattice gauge theories,” Nucl. Phys. B 375, 195–239 (1992).
  * Panagopoulos  _et al._ (2021) George Panagopoulos, Haralambos Panagopoulos, and Gregoris Spanoudes, “Two-loop renormalization and mixing of gluon and quark energy-momentum tensor operators,” Phys. Rev. D 103, 014515 (2021), arXiv:2010.02062 [hep-lat] .
  * Alexandrou  _et al._ (2020) Constantia Alexandrou, Andreas Athenodorou, Krzysztof Cichy, Arthur Dromard, Elena Garcia-Ramos, Karl Jansen, Urs Wenger, and Falk Zimmermann, “Comparison of topological charge definitions in Lattice QCD,” Eur. Phys. J. C 80, 424 (2020), arXiv:1708.00696 [hep-lat] .

##  I Supplemental Material

###  I.1 Results for the bare and renormalized momentum fractions

The bare momentum fractions can be extracted from the ratios in Eq. (2) by
performing a simultaneous fit for all $t_{\mathrm{s}}>t_{s,\textrm{low}}$. In
this way, from the set of all available $t_{s}$, $t_{s,\textrm{low}}$ marks
the lowest value within this set from where we choose to fit. Figs. 5 and 6
show the overlay of the raw ratio data for individual source-sink separations
$t_{\mathrm{s}}$ in the left-hand column and the fit for individual values of
$t_{\mathrm{s}}$ in the center for the pion and the kaon, respectively, using
the B ensemble as defined in Table 1. Finally, the right-hand column shows the
stability of the fit as we vary the lower bound $t_{s,\textrm{low}}$ of data
sets included in the simultaneous fit. In Tables 4–6 we compile all the bare
quantities extracted using this procedure for each ensemble. We show connected
contributions as well as disconnected contributions and the component of the
EMT used to extract those values for both the pion and the kaon.

In the case of the gluon contribution, we can perform two-state fits to the
two- and three-point functions and compare them with the plateau fit results.
In Fig. 7 we show the curves for the two-state fit and in Fig. 8 we show the
comparison between the values extracted from the two-state fit and the plateau
fit, and they are in agreement within errors. Because the errors increase very
fast for the two-state fit, we decided to use only the plateau fit results in
our analysis.

Figure 5: Complete decomposition of the bare momentum fraction carried by
quarks for the pion using the C ensemble. The two upper rows give the momentum
fraction $\langle x\rangle_{q}$ carried by $u+d$ (top for connected and second
row from top for disconnected), the 3rd row for the strange quark and the
lower row for the charm quark. In the first column we show the bare ratio for
different values of $t_{s}$, in the second column the individual fits per
$t_{s}$, and in the last column the plateau fit as a function of the lowest
value of $t_{s}$, $t_{s,\textrm{low}}$, used in the fits which demonstrates
the convergence as a function of $t_{s,\textrm{low}}$ varies. The band shows
the value after model averaging and the open symbol on the right panel
corresponds to the fit with the highest probability.  Figure 6: Complete
decomposition of the momentum carried by quarks in the kaon using the C
ensemble. The notation is the same as that of Fig.5.  Figure 7: Two state fit
for both pion and kaon at 10 stout-smearing steps. The bands in colors are
produced with the values of the fits to the correlators. The horizontal band
shows the value after model averaging. Figure 8: Plateau and two-state fit for
$\langle x\rangle_{g}$ for the pion(upper panel) and the kaon(lower panel).
The black squares represent the values per $t_{s,low}$ for the plateau fit and
the green squares the values for the two-state fit. Both fits agree, while the
plateau shows a good stability at higher values of $t_{s,low}$ the two-state
fit gains error faster. The open symbol shows the values with the highest
probability for each fit.

The renormalized results for the three ensembles are presented in Tables 7, 8,
and 9. The nominal difference of the renormalized results for the pion with
our previous analysis on ensemble B (cB211.072.64) in Ref. Alexandrou  _et
al._ (2021b) is rooted in the improvement of available statistics, especially
for the gluon operator insertion, as well as the improvement in the
determination of the renormalization and mixing coefficients.

In Fig. 9 we show the individual quark and gluon momentum fractions in the
pion and kaon at the continuum limit using the numbers of Table 2. The
continuum limit has been calculated by fitting the results for each ensemble
to a constant (including the three points and also the two points closer to
the continuum) and the linear fit in $a^{2}$. In Table 3 we present the result
of these fits and the probabilities associated to each one.

Table 3: Weights and values for each fit contributing to the continuum limit. Each row represents a different contribution to the momentum, from quarks or gluons, and then their sum for the pion and for the kaon. Each column represents the weight and value for each fit taken, with the first one the constant fit taking only two points, the second one the constant fit with three points and the last one the linear fit taking all three points. | Plateau 2 points | Plateau 3 points | Linear fit
---|---|---|---
$\langle x\rangle^{\pi}_{q,R}$ | $0.218,0.539(30)$ | $0.523,0.552(15)$ | $0.259,0.486(88)$
$\langle x\rangle^{\pi}_{g,R}$ | $0.304,0.403(35)$ | $0.504,0.375(21)$ | $0.191,0.398(94)$
$\langle x\rangle^{\pi}_{R}$ | $0.237,0.948(47)$ | $0.555,0.926(26)$ | $0.208,0.90(13)$
$\langle x\rangle^{K}_{q,R}$ | $0.208,0.609(18)$ | $0.554,0.612(10)$ | $0.238,0.642(55)$
$\langle x\rangle^{K}_{g,R}$ | $0.317,0.386(22)$ | $0.296,0.361(13)$ | $0.387,0.461(64)$
$\langle x\rangle^{K}_{R}$ | $0.211,0.998(28)$ | $0.310,0.973(17)$ | $0.478,1.115(85)$
Figure 9: The quark and gluon momentum fractions for the pion (upper panel) and kaon (lower panel) using the numbers in Table 2. Inner bars represent only the connected contributions, while the outer bars show the total, including disconnected contributions. Table 4: Bare values for the different contributions for the B (cB211.072.64) ensemble. Contribution | Operator | $\pi$ | $K$
---|---|---|---
$\langle x\rangle_{l}^{\mathrm{conn}}$ | $\bar{T}_{44}$ | 0.3929(44) | 0.1926(13)
$\langle x\rangle_{s}^{\mathrm{conn}}$ | $\bar{T}_{44}$ | \- - - | 0.2694(06)
$\langle x\rangle_{l}^{\mathrm{disc}}$ | $\bar{T}_{4k}$ | 0.0612(59) | 0.0510(34)
$\langle x\rangle_{s}^{\mathrm{disc}}$ | $\bar{T}_{4k}$ | 0.0441(58) | 0.0385(32)
$\langle x\rangle_{c}^{\mathrm{disc}}$ | $\bar{T}_{4k}$ | 0.0150(46) | 0.0137(23)
Table 5: Bare values for the different contributions to $\langle x\rangle$ for the C (cC211.060.80) ensemble. Contribution | Operator | $\pi$ | $K$
---|---|---|---
$\langle x\rangle_{l}^{\mathrm{conn}}$ | $\bar{T}_{44}$ | 0.369(06) | 0.183(02)
$\langle x\rangle_{s}^{\mathrm{conn}}$ | $\bar{T}_{44}$ | \- - - | 0.260(02)
$\langle x\rangle_{l}^{\mathrm{disc}}$ | $\bar{T}_{4k}$ | 0.066(13) | 0.044(07)
$\langle x\rangle_{s}^{\mathrm{disc}}$ | $\bar{T}_{4k}$ | 0.055(12) | 0.036(07)
$\langle x\rangle_{c}^{\mathrm{disc}}$ | $\bar{T}_{4k}$ | 0.023(10) | 0.011(05)
Table 6: Bare values for the different contributions to $\langle x\rangle$ for the D (cD211.054.96) ensemble. Contribution | Operator | $\pi$ | $K$
---|---|---|---
$\langle x\rangle_{l}^{\mathrm{conn}}$ | $\bar{T}_{44}$ | 0.357(03) | 0.180(02)
$\langle x\rangle_{s}^{\mathrm{conn}}$ | $\bar{T}_{44}$ | \- - - | 0.249(01)
$\langle x\rangle_{l}^{\mathrm{disc}}$ | $\bar{T}_{4k}$ | 0.041(11) | 0.055(07)
$\langle x\rangle_{s}^{\mathrm{disc}}$ | $\bar{T}_{4k}$ | 0.029(11) | 0.047(07)
$\langle x\rangle_{c}^{\mathrm{disc}}$ | $\bar{T}_{4k}$ | 0.010(19) | 0.045(11)
Table 7: Renormalized values for the different contributions for B (cB211.072.64) ensemble in the $\overline{\mathrm{MS}}$ scheme at a scale of 2 GeV. Contribution | $\pi$ | $K$
---|---|---
$\langle x\rangle_{l,R}$ | $0.495(09)$ | $0.264(05)$
$\langle x\rangle_{s,R}$ | $0.047(07)$ | $0.337(07)$
$\langle x\rangle_{c,R}$ | $0.014(05)$ | $0.013(03)$
$\langle x\rangle_{q,R}$ | $0.557(18)$ | $0.613(13)$
$\langle x\rangle_{g,R}$ | $0.360(25)$ | $0.346(16)$
Table 8: Renormalized values for the different contributions for C (cC211.060.80) ensemble in the $\overline{\mathrm{MS}}$ scheme at a scale of 2 GeV. Contribution | $\pi$ | $K$
---|---|---
$\langle x\rangle_{l,R}$ | $0.479(15)$ | $0.249(08)$
$\langle x\rangle_{s,R}$ | $0.059(14)$ | $0.327(08)$
$\langle x\rangle_{c,R}$ | $0.022(11)$ | $0.009(06)$
$\langle x\rangle_{q,R}$ | $0.560(38)$ | $0.585(21)$
$\langle x\rangle_{g,R}$ | $0.478(51)$ | $0.375(28)$
Table 9: Renormalized values for the different contributions for D (cD211.054.96) ensemble in the $\overline{\mathrm{MS}}$ scheme at a scale of 2 GeV. Contribution | $\pi$ | $K$
---|---|---
$\langle x\rangle_{l,R}$ | $0.459(15)$ | $0.270(09)$
$\langle x\rangle_{s,R}$ | $0.033(13)$ | $0.342(12)$
$\langle x\rangle_{c,R}$ | $0.012(22)$ | $0.052(13)$
$\langle x\rangle_{q,R}$ | $0.504(50)$ | $0.664(32)$
$\langle x\rangle_{g,R}$ | $0.337(48)$ | $0.404(35)$

###  I.2 Renormalization in the $\overline{\mathrm{MS}}$ scheme

The bare quark and gluon average momentum fractions have been renormalized
nonperturbatively by using the RI′/MOM scheme followed by perturbative
conversion to $\overline{\rm MS}$ at the scale $\mu=2$ GeV. Based on global
symmetries, the singlet quark and gluon EMTs in Eqs. (4 – 5) mix under
renormalization according to Eq. (7). Due to the breaking of rotational
symmetry on the lattice, $\bar{T}^{q(g)}_{\mu=\nu}$ and
$\bar{T}^{q(g)}_{\mu\neq\nu}$ belong to different representations of the
hypercubic group, and they are renormalized differently. The RI′/MOM
renormalization matrix for each case is calculated by the following
conditions:

${\begin{pmatrix}Z_{qq}^{s,{\rm RI}^{\prime}}&Z_{qg}^{{\rm RI}^{\prime}}\\\
Z_{gq}^{{\rm RI}^{\prime}}&Z_{gg}^{{\rm
RI}^{\prime}}\end{pmatrix}}^{-1}=\begin{pmatrix}\frac{1}{12}{\rm
Tr}[V_{qq}(V_{qq}^{\rm tree})^{-1}]&V_{qg}(V_{gg}^{\rm tree})^{-1}\\\
\frac{1}{12}{\rm Tr}[V_{gq}(V_{qq}^{\rm tree})^{-1}]&V_{gg}(V_{gg}^{\rm
tree})^{-1}\end{pmatrix},$ (11)

where

$\displaystyle V_{qq}$ $\displaystyle=$ $\displaystyle\frac{\langle
q(p)\bar{T}^{q}_{\mu\nu}\bar{q}(p)\rangle_{\rm amp}}{Z_{q}^{{\rm
RI}^{\prime}}},$ (12) $\displaystyle V_{qg}$ $\displaystyle=$
$\displaystyle\frac{\langle{\rm Tr}[A_{\rho}(p)A_{\rho}(-p)]\
\bar{T}^{q}_{\mu\nu}\rangle}{\langle{\rm
Tr}[A_{\rho}(p)A_{\rho}(-p)]\rangle},$ (13) $\displaystyle V_{gq}$
$\displaystyle=$ $\displaystyle\frac{\langle
q(p)\bar{T}^{g}_{\mu\nu}\bar{q}(p)\rangle_{\rm amp}}{Z_{q}^{{\rm
RI}^{\prime}}},$ (14) $\displaystyle V_{gg}$ $\displaystyle=$
$\displaystyle\frac{\langle{\rm Tr}[A_{\rho}(p)A_{\rho}(-p)]\
\bar{T}^{g}_{\mu\nu}\rangle}{\langle{\rm
Tr}[A_{\rho}(p)A_{\rho}(-p)]\rangle},$ (15)

and $V_{qq}^{\rm tree}$ ($V_{gg}^{\rm tree}$) is the tree-level value of
$V_{qq}$ ($V_{gg}$). The renormalization factor of the quark field
$Z_{q}^{{\rm RI}^{\prime}}$ is calculated in standard way through the quark
propagator Alexandrou  _et al._ (2017). Note that in $V_{gg}$ and $V_{qg}$, we
set $\rho\neq(\mu,\nu)$ and $p_{\rho}=0$ in order to cancel the mixing of EMTs
with gauge-noninvariant operators Caracciolo  _et al._ (1992), at least, to
the leading order. For the same reason, in the case $\mu=\nu$,
$\bar{T}_{\mu\nu}^{q(g)}$ is replaced by
$\bar{T}_{\mu\mu}^{q(g)}-\bar{T}_{\sigma\sigma}^{q(g)}$,
$\sigma\neq(\mu,\rho)$. The corresponding $\overline{\rm MS}$ renormalization
matrix is extracted through:

$\displaystyle{\begin{pmatrix}Z_{qq}^{s,\overline{\rm
MS}}&Z_{qg}^{\overline{\rm MS}}\\\ Z_{gq}^{\overline{\rm
MS}}&Z_{gg}^{\overline{\rm MS}}\end{pmatrix}}^{-1}$ $\displaystyle=$
$\displaystyle{\begin{pmatrix}Z_{qq}^{s,{\rm RI}^{\prime}}&Z_{qg}^{{\rm
RI}^{\prime}}\\\ Z_{gq}^{{\rm RI}^{\prime}}&Z_{gg}^{{\rm
RI}^{\prime}}\end{pmatrix}}^{-1}$ (16) $\displaystyle\times$
$\displaystyle{\begin{pmatrix}C_{qq}^{\overline{\rm MS},{\rm
RI}^{\prime}}&C_{qg}^{\overline{\rm MS},{\rm RI}^{\prime}}\\\
C_{gq}^{\overline{\rm MS},{\rm RI}^{\prime}}&C_{gg}^{\overline{\rm MS},{\rm
RI}^{\prime}}\end{pmatrix}}^{-1},$

where $C_{ij}^{\overline{\rm MS},{\rm RI}^{\prime}}$ are the conversion-matrix
coefficients calculated up to two-loop order Panagopoulos  _et al._ (2021).
For the renormalization of flavor-nonsinglet combinations $\langle
x\rangle_{u+d-2s}$ and $\langle x\rangle_{u+d+s-3c}$, the nonsinglet
renormalization factor $Z_{qq}$ is used, which is calculated from the
connected part of $V_{qq}$ (see Alexandrou  _et al._ (2017)).

In this work, the diagonal elements $Z_{qq}^{s,\overline{\rm MS}}$ and
$Z_{gg}^{\overline{\rm MS}}$ are evaluated non-perturbatively, while the off-
diagonal $Z_{qg}^{\overline{\rm MS}}$ and $Z_{gq}^{\overline{\rm MS}}$ are
calculated in one-loop lattice perturbation theory. For the non-perturbative
calculations, we have employed high-statistical ($\sim 10000-30000$
configurations) $N_{f}=4$ mass-degenerate ensembles generated by ETMC at the
same lattice spacings of the three physical ensembles (see Table 1). The
parameters of the $N_{f}=4$ ensembles are given in Table 10.

Table 10: $N_{f}=4$ ensembles and their parameters used for the renormalization: $\beta=(2N/g^{2})$, lattice volume ($L^{3}\times T$), twisted-mass parameter $(a\mu)$, hopping parameter $\kappa$ and clover coefficient $c_{\rm SW}$. Ensemble |  |  |  |  |
---|---|---|---|---|---
Label | $\beta$ | $L^{3}\times T$ | $a\mu$ | $\kappa$ | $c_{\rm SW}$
B | 1.778 | $12^{3}\times 24$ | 0.006 | 0.1393050 | 1.6900
C | 1.836 | $24^{3}\times 48$ | 0.005 | 0.1386735 | 1.6452
D | 1.900 | $24^{3}\times 48$ | 0.004 | 0.1379346 | 1.6112

An important aspect of our renormalization program is the improvement of the
nonperturbative estimates by subtracting tree-level or one-loop discretization
errors Alexandrou  _et al._ (2017). Stout smearing is employed in the gluon
links of $\bar{T}_{\mu\nu}^{g}$ with the same set of smearing levels as used
during the computation of the matrix elements for ensembles B, C and D. In the
perturbative calculation of $Z_{gq}^{\overline{\rm MS}}$, we restrict the
stout smearing to three steps to manage the rapid growth of algebraic
expressions, which can result in millions of terms when more stout-smearing
steps are included. We have also employed Wilson flow (instead of stout
smearing) in the nonperturbative calculations and we have confirmed the
equivalence of the two smearing methods Alexandrou  _et al._ (2020) when
$\tau=\rho_{\mathrm{stout}}\,N_{\mathrm{stout}}$ ($\tau$ is the flow time, and
$\rho_{\mathrm{stout}}$, $N_{\mathrm{stout}}$ are the stout parameter and
number of stout-smearing steps, respectively). In table 11, we collect our
results for the renormalization coefficients that enter our analysis.

Table 11: Renormalization coefficients in the $\overline{\rm MS}$ scheme at 2 GeV for the three ensembles of Table 10 and for different tensor structures. For the computation of $Z_{gg}$, 5-10 steps of stout smearing are employed. The error quoted in the parenthesis is statistical; $Z_{qg}$ and $Z_{gq}$ have no statistical error because they are calculated in perturbation theory. |  Ensemble |  Ensemble |  Ensemble
---|---|---|---
|  B |  C |  D
$Z_{qq}\ (\mu=\nu)$ |  1.0982(3) |  1.1164(3) |  1.1564(1)
$Z_{qq}\ (\mu\neq\nu)$ |  1.1228(2) |  1.1481(3) |  1.1825(1)
$Z_{qq}^{\rm s}\ (\mu=\nu)$ |  1.0888(72) |  1.0987(99) |  1.1636(88)
$Z_{qq}^{\rm s}\ (\mu\neq\nu)$ |  1.0996(15) |  1.1323(12) |  1.1740(12)
$Z_{qg}\ (\mu\neq\nu)$ |  -0.0106 |  -0.0145 |  -0.0188
$Z_{gq}\ (\mu=\nu)$ |  0.0658 |  0.0410 |  0.0139
$Z_{gq}\ (\mu\neq\nu)$ |  0.0772 |  0.0521 |  0.0246
$Z_{gg}^{n_{\rm st}=5}\ (\mu\neq\nu)$ |  0.739(8) |  0.727(21) |  0.792(35)
$Z_{gg}^{n_{\rm st}=6}\ (\mu\neq\nu)$ |  0.744(7) |  0.738(19) |  0.799(32)
$Z_{gg}^{n_{\rm st}=7}\ (\mu\neq\nu)$ |  0.755(7) |  0.753(18) |  0.812(30)
$Z_{gg}^{n_{\rm st}=8}\ (\mu\neq\nu)$ |  0.767(7) |  0.769(17) |  0.826(28)
$Z_{gg}^{n_{\rm st}=9}\ (\mu\neq\nu)$ |  0.778(6) |  0.784(16) |  0.839(26)
$Z_{gg}^{n_{\rm st}=10}\ (\mu\neq\nu)$ |  0.790(6) |  0.797(16) |  0.850(25)

###  I.3 Lattice setup

We provide additional information about the ETMC ensembles Alexandrou  _et
al._ (2018) in Tab. 12, together with volume of statistics used in this work.

Table 12: Statistics per ensemble used in this work. $n_{\mathrm{2pt}}$ and $n_{\mathrm{3pt}}$ give the number of stochastic samples produced per gauge configuration for two-point and (connected) three-point functions, respectively. In both cases stats is the total statistics for the given ensemble. Label | Ensemble name | $n_{\mathrm{conf}}$ | $n_{\mathrm{2pt}}$ | stats | $n_{\mathrm{3pt}}$ | stats
---|---|---|---|---|---|---
B | cB211.072.64 | 755 | 400 | 302,000 | 8 | 6040
C | cC211.060.80 | 400 | 240 | 96,000 | 8 | 3200
D | cD211.054.96 | 499 | 600 | 299,400 | 8 | 3992

The three different types of diagrams generated by the EMT insertions
$\bar{T}^{q}$ and $\bar{T}^{g}$ into the pion and kaon correlation function
are depicted in Fig. 10. We distinguish the (quark-)connected diagram $(I)$,
and the disconnected diagrams with a quark loop $(II)$ from inserting
$\bar{T}^{q}$ and with the gluonic operator insertion $\bar{T}^{g}$ in
$(III)$.

All quark-line (sub-)diagrams for three-point and two-point correlation
functions are produced with stochastic timeslice-to-all propagators with
randomly chosen source timeslices.

Figure 10: Quark-flow diagrams for three-point function with quark and gluon
energy momentum tensor operator insertion: $(I)$ connected and $(II)$
disconnected quark EMT insertion diagram, $(III)$ disconnected gluon EMT
insertion diagram.

In case of the connected three-point function diagram $(I)$ we apply another
sequential inversion through the hadron vertex at sink time.

To compute the diagrams involving strange or charm quark flow we use a mixed
action for the strange and charm quark: to avoid the technical complication of
strange and charm quark mixing inherent in the twisted mass sea quark action
we introduce additional valence quark doublets for both strange and charm
flavor in analogy to the light quark doublet with bare twisted masses of
opposite sign within each doublet. These bare twisted valence quark masses,
$\mu_{s}$ and $\mu_{c}$, are tuned such that the $\Omega^{-}$ and
$\Lambda_{c}$ baryons take their physical values in the continuum limit,
leading to the parameters given in Table 13.

Table 13: The tuned strange $\mu_{s}$ and charm $\mu_{c}$ bare masses for the three ensembles | B | C | D
---|---|---|---
$a\mu_{s}$ | $0.0186$ | $0.01615$ | $0.0136$
$a\mu_{c}$ | $0.2490$ | $0.2060$ | $0.1660$

Observables derived from the mixed action retain automatic $\mathcal{O}(a)$
improvement with unitarity violations scaling equally at last at second order
in the lattice spacing in the continuum limit.

Tables 12 and 14 list the number of stochastic samples for two- and three-
point functions per gauge configuration together with the total achieved
statistics.

The quark loops for diagram $(II)$ with covariant derivative insertion are
produced with hierarchical probing and full dilution in spin and color degrees
of freedom. For the light quark loops of ensembles B and C additional data
sets were produced using low-mode deflation on top of hierarchical probing.
The parameters for the loop computation per gauge configuration are given in
Table 14.

Table 14: For each quark flavor we give in this order the number of deflated low modes $N_{\rm defl}$, the number of stochastic sources $N_{\rm srcs}$ and the probing distance $N_{col}$. We note that for the B and C ensembles, we have two sets of loops for the light quarks, which are averaged at the level of the analysis. Ensemble | Light loops | Strange loops | Charm loops
---|---|---|---
$N_{\rm defl}$ | $N_{\rm srcs}$ | $N_{\rm col}$ | $N_{\rm defl}$ | $N_{\rm srcs}$ | $N_{\rm col}$ | $N_{\rm defl}$ | $N_{\rm srcs}$ | $N_{\rm col}$
B | 0 | 2 | 8 | 0 | 2 | 8 | 0 | 12 | 4
| 200 | 1 | 8 |  |  |  |  |  |
C | 0 | 2 | 8 | 0 | 4 | 8 | 0 | 1 | 8
| 450 | 1 | 8 |  |  |  |  |  |
D | 0 | 8 | 8 | 0 | 4 | 8 | 0 | 1 | 8

For completeness we give the connection between the ratios of two- and three-
point functions with EMT insertion $\bar{T}^{q},\,\bar{T}^{g}$ for tensor
components used in this work and the bare $\langle x\rangle$ values in Eq.
(17).

$\displaystyle R_{44}^{X}(t_{\mathrm{ins}},t_{\mathrm{s}};\,\mathbf{0})$
$\displaystyle\overset{t_{\mathrm{s}},t_{\mathrm{ins}}\mathrm{~{}large}}{\longrightarrow}-\frac{3}{4}\,\langle
x\rangle_{X}\,m_{h}\,,$ $\displaystyle
R_{4k}^{X}(t_{\mathrm{ins}},t_{\mathrm{s}};\,\mathbf{p})$
$\displaystyle\overset{t_{\mathrm{s}},t_{\mathrm{ins}}\mathrm{~{}large}}{\longrightarrow}\langle
x\rangle_{X}\,p_{k}\,.$ (17)

In both cases the relation holds up to excited state contamination and up to
the finite-time wrapping correction factor
$[1+\exp(-E_{h}(T-2\,t_{\mathrm{s}}))]^{-1}$ as given in Eq. (3).