# Nonlocal chiral contributions to generalized parton distributions of the
proton at nonzero skewness

Zhengyang Gao  Institute of High Energy Physics, Chinese Academy of Sciences,
Beijing 100049, China  College of Physics Sciences, University of Chinese
Academy of Sciences, Beijing 100049, China Fangcheng He  Department of Physics
and Astronomy, Stony Brook University, Stony Brook, New York 11790, USA
Chueng-Ryong Ji  Department of Physics, North Carolina State University,
Raleigh, North Carolina 27695, USA  W. Melnitchouk  Jefferson Lab, Newport
News, Virginia 23606, USA  Y. Salamu  School of Physics and Electrical
Engineering, Kashi University, Kashi, 844000, Xinjiang, China  P. Wang
Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049,
China  College of Physics Sciences, University of Chinese Academy of Sciences,
Beijing 100049, China

###### Abstract

We compute the one-loop contributions to spin-averaged generalized parton
distributions (GPDs) in the proton from pseudoscalar mesons with intermediate
octet and decuplet baryon states at nonzero skewness. Our framework is based
on nonlocal covariant chiral effective theory, with ultraviolet divergences
regularized by introducing a relativistic regulator derived consistently from
the nonlocal Lagrangian. Using the splitting functions calculated from the
nonlocal Lagrangian, we find the nonzero skewness GPDs from meson loops by
convoluting with the phenomenological pion GPD and the generalized
distribution amplitude, and verify that these satisfy the correct
polynomiality properties. We also compute the lowest two moments of GPDs to
quantify the meson loop effects on the Dirac, Pauli and gravitational form
factors of the proton.

††preprint: JLAB-THY-24-4069

##  I Introduction

Generalized parton distributions (GPDs) contain unique information about the
3-dimensional structure of hadrons, in terms of their fundamental quark and
gluon (or parton) constituents. They describe the distributions in momentum
and position space of partons carrying a specific fraction $x$ of the hadron’s
light-front momentum and squared four-momentum transfer $t$, and interpolate
between collinear parton distribution functions (PDFs) in the forward limit
and elastic form factors when integrated over $x$. The latter include the
Dirac and Pauli form factors, from the lowest moments of the GPDs, and
gravitational form factors from the $x$-weighted GPD moments Ji (1998) (for
reviews of GPDs see, e.g., Refs. Goeke  _et al._ (2001); Belitsky and
Radyushkin (2005)). A further property of GPDs is that they obey different
evolution equations in different kinematic regions: either the DGLAP
(Dokshitzer-Gribov-Lipatov-Altarelli-Parisi) region as for collinear PDFs,
where both the struck and returning partons carry positive momentum fractions
of the initial proton momentum Altarelli and Parisi (1977); Gribov and Lipatov
(1972); Dokshitzer (1977), or the ERBL (Efremov-Radyushkin-Brodsky-Lepage)
region, which can be interpreted as describing a quark-antiquark pair emerging
from the proton, as appropriate for distribution amplitudes (DAs) Efremov and
Radyushkin (1980); Lepage and Brodsky (1980, 1979).

On the experimental side, data from processes such as deeply-virtual Compton
scattering (DVCS) Ji (1997) and hard exclusive meson production (HEMP)
Radyushkin (1996); Mankiewicz  _et al._ (1998); Collins and Freund (1999),
including from the H1 Adloff  _et al._ (2001); Aktas  _et al._ (2005); Aaron
_et al._ (2009) and ZEUS Breitweg  _et al._ (1999); Chekanov  _et al._ (2003)
collaborations at HERA, the HERMES fixed target experiment Airapetian  _et
al._ (2001, 2011, 2012), the COMPASS Collaboration at CERN d’Hose  _et al._
(2004); Fuchey (2015), and from the CLAS Stepanyan  _et al._ (2001);
Bedlinskiy  _et al._ (2012); Seder  _et al._ (2015); Jo  _et al._ (2015);
Hattawy  _et al._ (2019); Burkert  _et al._ (2021) and Hall A Defurne  _et
al._ (2015); Georges  _et al._ (2022) collaborations at Jefferson Lab, have
been used to provide indirect information on GPDs through the Compton form
factors. The latter are given as convolutions of GPDs with hard scattering
kernels derived within QCD factorization Collins and Freund (1999); Belitsky
_et al._ (2002). The extraction of GPDs and their moments has also been a
motivation for planned experimental programs at future facilities such as the
Electron-Ion Collider (EIC) Accardi  _et al._ (2016); Avakian (2019).

The extraction of GPDs directly from experimental data is of course the
cornerstone of the quest to determine the 3-dimensional structure of the
nucleon. Pioneering studies have already been made by several groups worldwide
on this effort Diehl  _et al._ (2005); Kumerički and Müller (2010); Goldstein
_et al._ (2011); Guo  _et al._ (2022, 2023), although direct, model-
independent extractions remain a formidable challenge. Among some early
studies, Diehl et al. Diehl  _et al._ (2005) used a simple empirical
parametrization of the $x$ and $t$ dependence of GPDs at zero skewness, with
forward collinear PDFs as input. Kumerički and Müller Kumerički and Müller
(2010) studied DVCS at small values of the Bjorken-$x$ variable, using
conformal integral GPD moments, performing a first model-dependent extraction
of the unpolarized GPD $H$ from HERA and Jefferson Lab DVCS data. Goldstein et
al. Goldstein  _et al._ (2011) carried out a global analysis of DVCS
observables, together with nucleon elastic form factors and deep-inelastic
scattering measurements, using a flexible parametrization of GPDs inspired by
a hybrid model of the nucleon as a quark-diquark system with Regge behavior.
Most recently, Guo et al. Guo  _et al._ (2022, 2023) employed the GUMP moment
parametrization, motivated by the complex conformal spin partial wave
expansion method of Müller and Schäfer Müller and Schäfer (2006), to perform a
global analysis of GPDs from DVCS data, constrained by input on PDFs and
elastic form factors, as well as recent lattice QCD calculations Alexandrou
_et al._ (2020, 2022).

Notwithstanding these important developments, in a seminal paper Bertone et
al. Bertone  _et al._ (2021) recently pointed out, within a next-to-leading
order QCD analysis, that a critical limitation of processes such as DVCS is
their inability to uniquely determine the $x$ dependence of the GPDs due to
the presence of so-called “shadow GPDs.” The shadow GPDs are a set of
solutions to the inverse problem of extracting GPDs from DVCS data, which
depend on the process and do not possess a unique solution at a fixed scale
$Q^{2}$. Moffat et al. Moffat  _et al._ (2023) further investigated the extent
to which QCD evolution can provide constraints on the shadow GPDs, observing
that given empirical data over a sufficiently wide range of skewness $\xi$ and
scale $Q^{2}$, the $Q^{2}$ evolution may in principle constrain the shadow
GPDs, but over a rather limited range of $x$ and $\xi$.

More recently, Qiu and Yu Qiu and Yu (2023, 2022) explored a novel new set of
exclusive processes, including single diffractive hard exclusive
photoproduction, as a means of more directly constraining the $x$ dependence
of GPDs. Such processes can be shown to be factorized into process-independent
GPDs and perturbatively calculable, infrared safe hard coefficients Qiu and Yu
(2023). They can be obtained from polarized photoproduction cross sections, as
well as asymmetries constructed from photon polarization and hadron spin,
which could in future be measured, for example, at Jefferson Lab Hall D. Other
related processes, such as exclusive production of high transverse momentum
photons in $\pi N$ scattering, can be factorized into GPDs and pion DAs
convoluted with short distance hard kernels Qiu and Yu (2022), if the photons’
transverse momentum $q$ is $\gg\Lambda_{\rm QCD}$, with corrections suppressed
by powers of $1/q$.

Along with the empirical determinations of GPDs from data, developments in
lattice QCD have allowed access to nucleon structures from first principles
calculations, providing complementary information that is often difficult to
obtain from experiment Lin (2021); Bhattacharya  _et al._ (2022); Alexandrou
_et al._ (2020, 2006). In particular, considerable effort has been devoted to
directly computing the parton momentum fraction dependence of distributions,
using either the quasi-parton distribution Ji (2013), pseudo-distribution
Orginos  _et al._ (2017), or lattice good cross sections Ma and Qiu (2018a, b)
approaches.

While the lattice continues to make progress with improved control over
systematic uncertainties in the simulations, it can also be useful to explore
possible insights from model calculations, many of which have been performed
over the past 25 years Ji  _et al._ (1997); Petrov  _et al._ (1998); Penttinen
_et al._ (2000); Goeke  _et al._ (2001); Choi  _et al._ (2001, 2002); Tiburzi
and Miller (2002); Theussl  _et al._ (2004); Boffi  _et al._ (2003);
Schweitzer  _et al._ (2002); Scopetta and Vento (2004); Ossmann  _et al._
(2005); Wakamatsu and Tsujimoto (2005); Mineo  _et al._ (2005); Pasquini and
Boffi (2007); Goeke  _et al._ (2008). Even though the direct connection to QCD
is not always transparent in such models, they can nevertheless sometimes
provide glimpses into some of the qualitative features and characteristics of
GPDs.

In between the phenomenological models and lattice QCD simulations,
constraints from chiral effective field theory (EFT) have been used to make
predictions for various sea quark flavor asymmetries in the nucleon Thomas
(1983); Thomas  _et al._ (2000); Chen and Ji (2001); Salamu  _et al._ (2015);
Luan and Lu (2023). These are based on the observation that the long-range
structure of the nucleon has contributions from the pseudoscalar meson cloud,
associated with the model-independent leading nonanalytic behavior of
observables expanded in a series of $m_{\pi}$. Operationally, the sea quark
distributions can be obtained through the convolution of a probability of the
nucleon to split into a virtual meson and recoil baryon (“splitting
function”), with the valence quark distribution of the meson. The splitting
function can be calculated from EFT, in which its moments are expanded as a
series in the pseudoscalar meson mass ${\cal O}(m_{\phi}/\Lambda_{\chi})$ or
small external momentum ${\cal O}(q/\Lambda_{\chi})$. Traditional EFT
calculations based on dimensional regularization have found it challenging to
describe lattice results in the large-$Q^{2}$ region, or with pion mass 400
MeV Young  _et al._ (2003); Leinweber  _et al._ (2004), since the short-
distance effect arsing from the loop integral leads to a poor convergence of
the chiral expansion Donoghue  _et al._ (1999); Leinweber  _et al._ (2004). To
improve the convergence, finite range regularization has been proposed and
applied for the extrapolation of various quantities calculated by lattice QCD,
including the vector meson mass, magnetic moments, magnetic and strange form
factors, charge radii, and moments of PDFs and GPDs, from large pion masses to
the physical mass Young  _et al._ (2003); Leinweber  _et al._ (2004); Wang
_et al._ (2007); Allton  _et al._ (2005); Wang  _et al._ (2009, 2012); Hall
_et al._ (2014); Shanahan  _et al._ (2013, 2014). Here the three dimensional
regulator is included in the loop integral and can be regarded as partially
resumming the contribution from the higher order interaction. In a similar
vein, a nonlocal chiral effective Lagrangian has also been proposed Wang
(2011); He and Wang (2018a, b), in which the covariant regulator is
automatically generated from the Lagrangian. One can obtain reasonable
descriptions of the nucleon electromagnetic and strange form factors at large
$Q^{2}$ using this method, which has also been applied to calculate the
strange–antistrange PDF asymmetry $s-\bar{s}$ Salamu  _et al._ (2019a, b), the
sea quark transverse momentum dependent (TMD) Sivers function He and Wang
(2019), and zero-skewness GPDs He  _et al._ (2022) in the proton. More details
about the nonlocal chiral effective theory can be found in Ref. Wang  _et al._
(2023). Recently, Copeland and Mehen Copeland and Mehen (2024) also discussed
a framework for matching TMD PDFs onto chiral effective theory operators in
terms of TMD hadronic distribution functions calculated in chiral perturbation
theory.

In this paper, we extend our previous analysis He  _et al._ (2022) of zero-
skewness GPDs in the proton within the nonlocal chiral effective theory to the
nonzero skewness case. In Sec. II we present the basic theoretical framework
used in our analysis, including the definition of unpolarized GPDs (Sec. II.1)
and the nonlocal meson-baryon interaction that underlies our calculation (Sec.
II.2). The one-loop nucleon splitting into meson plus octet and decuplet
baryon splitting functions are given in Sec. II.3 from the set of rainbow,
tadpole, and bubble diagrams that contribute to the antiquark GPDs in the
nucleon. The convolution formalism is derived in Sec. III, where we outline
the differences with the zero skewness GPD case. In particular, we separate
the contribution to DGLAP and ERBL region, and prove that such a convolution
formula can generate the skewness-independent form factors after integrating
over $x$. Numerical results for the two and three dimensions splitting
functions and light antiquark GPDs using the convolution formulas are
presented in Sec. IV. Finally, Sec. V summarises our results and anticipates
future extensions of this analysis. In Appendix A, we also demonstrate the
GPDs’ $n$-th moments obtained with convolution formula satisfy the general
polynomial property. In Appendix B, we summarize the explicit expressions for
the splitting functions and the coupling constants presented in this study.

##  II Theoretical framework

In this section we summarize the theoretical framework on which our analysis
is based. We begin with a summary of the definitions of the spin-averaged GPDs
that are the focus of this analysis, followed by a discussion of the
contributions to the GPDs from pseudoscalar meson loops, formulated within
chiral effective theory. We summarize the pertinent aspects of the chiral
Lagrangian describing the meson-nucleon interaction, and present results for
the nucleon to baryon $+$ meson splitting functions at nonzero skewness.

###  II.1 Unpolarized GPDs

The spin-averaged GPDs for a quark of flavor $q$ in a nucleon are defined in
terms of the Fourier transform of the nonforward matrix elements of quark
bilocal field operators, taken between nucleon states with initial momentum
$p$ and final momentum $p^{\prime}$ Ji (1997),

$\displaystyle\int_{-\infty}^{\infty}\frac{\differential{\lambda}}{2\pi}e^{-ix\lambda}\langle
p^{\prime}|\bar{\psi}_{q}(\tfrac{1}{2}\lambda
n)\not{n}\,\psi_{q}(-\tfrac{1}{2}\lambda n)|p\rangle$ $\displaystyle=$
$\displaystyle\bar{u}(p^{\prime})\Big{[}\not{n}H_{q}(x,\xi,t)+\frac{i\sigma^{\mu\nu}n_{\mu}\Delta_{\nu}}{2M}\,E_{q}(x,\xi,t)\Big{]}u(p),$

where the light-cone vector $n_{\mu}$ projects the “plus” component of
momentum, and $\lambda$ is a dimensionless parameter. Invariance under Lorentz
transformations requires the Dirac ($H_{q}$) and Pauli ($E_{q}$) GPDs to be
functions of the light-cone momentum fraction $x$ of the nucleon carried by
the initial quark with momentum $k_{q}$,

$x\,\equiv\,\frac{k_{q}^{+}}{P^{+}},$ (2)

and the skewness parameter $\xi$,

$\xi\,\equiv\,-\frac{\Delta^{+}}{2P^{+}},$ (3)

where

$\displaystyle P$ $\displaystyle=$ $\displaystyle\frac{1}{2}(p+p^{\prime}),$
(4a) $\displaystyle\Delta$ $\displaystyle=$ $\displaystyle p^{\prime}-p$ (4b)

are the average and difference of the initial and final nucleon momenta,
respectively. We define the light-cone components of a four-vector $k^{\mu}$
by $k^{\pm}=\frac{1}{\sqrt{2}}(k^{0}\pm k^{3})$, and for convenience choose a
symmetric frame of reference for the initial and final nucleon momenta,

$\displaystyle p^{\mu}$ $\displaystyle=$
$\displaystyle\Big{(}P^{+}-\frac{\Delta^{+}}{2},P^{-}-\frac{\Delta^{-}}{2},-\frac{\bm{\Delta}_{\perp}}{2}\Big{)},$
(5a) $\displaystyle p^{\prime\mu}\ $ $\displaystyle=$
$\displaystyle\Big{(}P^{+}+\frac{\Delta^{+}}{2},P^{-}+\frac{\Delta^{-}}{2},\quad\frac{\bm{\Delta}_{\perp}}{2}\Big{)}.$
(5b)

In addition to being functions of $x$ and $\xi$, the GPDs also depend on the
hadronic four-momentum transfer squared, $t\equiv\Delta^{2}$, as well as on
the scale $Q^{2}$, which is typically taken to be the four-momentum transfer
squared from the incident lepton. (The dependence on the latter is suppressed
for notational convenience.) The momentum transfer should satisfy the
kinematic constraint $-\Delta^{2}>4\xi^{2}M^{2}/(1-\xi^{2})$ in order to
guarantee positive transferred momentum squared, $\bm{\Delta_{\perp}}^{2}>0$.

Integrating the GPDs $H_{q}$ and $E_{q}$ over the quark light-cone momentum
fraction $x$, the Dirac and Pauli form factors for a quark flavor $q$ can be
expressed as

$\displaystyle\int_{-1}^{1}\,\differential{x}H_{q}(x,\xi,t)$ $\displaystyle=$
$\displaystyle F_{1}^{q}(t),$ (6a)
$\displaystyle\int_{-1}^{1}\,\differential{x}E_{q}(x,\xi,t)$ $\displaystyle=$
$\displaystyle F_{2}^{q}(t).$ (6b)

Note that the $\xi$ dependence vanishes in the integrals of the GPDs over $x$.
Furthermore, the gravitational form factors $A^{q}$, $B^{q}$ and $C^{q}$ can
be obtained from the first moments of the $H_{q}$ and $E_{q}$ GPDs as

$\displaystyle\int_{-1}^{1}\,\differential{x}x\,H_{q}(x,\xi,t)$
$\displaystyle=$ $\displaystyle A^{q}(t)+(2\xi)^{2}\,C^{q}(t),$ (7a)
$\displaystyle\int_{-1}^{1}\,\differential{x}x\,E_{q}(x,\xi,t)$
$\displaystyle=$ $\displaystyle B^{q}(t)-(2\xi)^{2}\,C^{q}(t),$ (7b)

where now explicit dependence on the skewness $\xi$ now appears on the right
hand side.

###  II.2 Chiral loop contributions to GPDs

In our analysis we consider the contributions to the sea quark and antiquark
GPDs in the proton arising from the virtual pseudoscalar meson cloud dressing
of the bare baryon, as a first step towards a complete calculation to one-loop
order. Typically, in calculations of meson loop contributions to sea quark and
antiquark asymmetries, assuming that the undressed proton has a flavor
symmetric sea Salamu  _et al._ (2015), the meson coupling diagrams in Fig. 1
are the dominant source of differences between sea quark and antiquark PDFs
and GPDs.

Figure 1: One-loop diagrams for the proton to pseudoscalar meson (dashed
lines) and octet baryon (solid lines) or decuplet baryon (double solid lines)
splitting functions. The crossed circles ($\otimes$) represent the interaction
with the external vector field from the minimal substitution, and the gray
square (${{\color[rgb]{.5,.5,.5}\blacksquare}}$) denotes the last interaction
term in Eq. (10).

In our previous analysis He  _et al._ (2022) we considered the nonforward
nucleon to to virtual meson and recoil baryon splitting functions at zero
skewness. Generalizing to nonzero $\xi$, the splitting functions can be
written as

$\bar{u}(p^{\prime})\Big{[}\gamma^{+}f(y,\xi,t)+\frac{i\sigma^{+\nu}\Delta_{\nu}}{2M}\,g(y,\xi,t)\Big{]}u(p)=\int\differential[4]{k}\widetilde{\Gamma}^{+}(k)\,\delta\Big{(}y+\xi-\frac{k^{+}}{P^{+}}\Big{)}\,\equiv\,\Gamma^{+},$
(8)

where $f(y,\xi,t)$ and $g(y,\xi,t)$ are the corresponding Dirac- and Pauli-
like splitting functions, and $y$ is the light-cone momentum fraction of the
nucleon carried by the intermediate state hadron. The operator
$\widetilde{\Gamma}^{\mu}(k)$ is related to the matrix elements of the
hadronic current $J^{\mu}$,

$\langle
N(p^{\prime})|J^{\mu}|N(p)\rangle\,=\,\bar{u}(p^{\prime})\Big{[}\gamma^{\mu}F_{1}^{N}(t)+\frac{i\sigma^{\mu\nu}\Delta_{\nu}}{2M}F_{2}^{N}(t)\Big{]}u(p)\,\equiv\int\\!\differential[4]{k}\,\widetilde{\Gamma}^{\mu}(k),$
(9)

where $J^{\mu}$ is given in terms of its individual quark flavor contributions
as in Ref. He  _et al._ (2022). One can also verify that the Dirac and Pauli
form factors can be obtained by integrating the splitting functions over $y$.

To compute the hadronic splitting functions $f$ and $g$, we use the standard
chiral SU(2)${}_{L}\times$SU(2)R effective Lagrangian, the details of which
can be found, for example, in Refs. Jenkins and Manohar (1991); Bernard
(2008); Hacker  _et al._ (2005). In the one-loop approximation used in our
analysis, the chiral Lagrangian can be expanded to reveal the interaction
terms, given by

$\displaystyle{\cal L}_{\rm int}$ $\displaystyle=$
$\displaystyle\frac{g_{A}}{2f_{\pi}}\left(\bar{p}\,\gamma^{\mu}\gamma_{5}p\,\partial_{\mu}\pi^{0}+\sqrt{2}\,\bar{p}\,\gamma^{\mu}\gamma_{5}n\,\partial_{\mu}\pi^{+}\right)$
(10) $\displaystyle+$ $\displaystyle\frac{\cal
C}{\sqrt{12}f_{\pi}}\left(-2\,\bar{p}\,\Theta^{\nu\mu}\Delta_{\mu}^{+}\,\partial_{\nu}\pi^{0}-\sqrt{2}\,\bar{p}\,\Theta^{\nu\mu}\Delta_{\mu}^{0}\,\partial_{\nu}\pi^{+}+\sqrt{6}\,\bar{p}\,\Theta^{\nu\mu}\Delta_{\mu}^{++}\,\partial_{\nu}\pi^{-}\right)$
$\displaystyle+$
$\displaystyle\frac{i}{4f_{\pi}^{2}}\,\bar{p}\,\gamma^{\mu}p\,(\pi^{+}\partial_{\mu}\pi^{-}-\pi^{-}\partial_{\mu}\pi^{+})+\frac{2i(b_{10}+b_{11})}{f_{\pi}^{2}}\,\bar{p}\,\sigma^{\mu\nu}p\,\partial_{\mu}\pi^{+}\,\partial_{\nu}\pi^{-}+{\rm
H.c.},$

where H.c. refers to the Hermitian conjugate, and $f_{\pi}=92.4$ MeV is the
pion decay constant. The axial vector coupling constant for the octet baryon
is set to the standard value $g_{A}=1.26$, and the octet-decuplet transition
coupling ${\cal C}$ is chosen to be $-\frac{6}{5}g_{A}$ from SU(6) spin-flavor
symmetry. The octet–decuplet transition operator $\Theta^{\mu\nu}$ is defined
as Nath  _et al._ (1971)

$\Theta^{\mu\nu}=g^{\mu\nu}-\gamma^{\mu}\gamma^{\nu}.$ (11)

The coefficients $b_{10}$ and $b_{11}$ in Eq. (10) for the nucleon-pion
contact interaction have previously been determined to be $b_{10}=1.24$ GeV-1,
and $b_{11}=0.46$ GeV-1 Kubis and Meißner (2001). In our previous work He  _et
al._ (2022) we discussed in detail the construction of the nonlocal chiral
Lagrangian; here we simply present the final nonlocal interaction ${\cal
L}^{\rm(nonloc)}_{\rm int}$ we used in this calculation,

$\displaystyle{\cal L}^{\rm(nonloc)}_{\rm int}(x)$ $\displaystyle=$
$\displaystyle\bar{p}(x)\bigg{(}\frac{g_{A}}{2f_{\pi}}\,\gamma^{\mu}\gamma_{5}p(x)-\frac{\cal
C}{\sqrt{3}f_{\pi}}\,\Theta^{\mu\nu}\Delta^{+}_{\nu}(x)\bigg{)}\\!\int\\!\differential[4]{a}\,F(a)\,\partial_{\mu}\pi^{0}(x+a)$
(12) $\displaystyle+\
\bar{p}(x)\bigg{(}\frac{g_{A}}{\sqrt{2}f_{\pi}}\,\gamma^{\mu}\gamma_{5}n(x)-\frac{\cal
C}{\sqrt{6}f_{\pi}}\,\Theta^{\mu\nu}\Delta^{0}_{\nu}(x)\bigg{)}\\!\int\\!\differential[4]{a}\,F(a)\,\partial_{\mu}\pi^{+}(x+a)$
$\displaystyle+\ \frac{\cal
C}{\sqrt{2}f_{\pi}}\,\bar{p}(x)\Theta^{\mu\nu}\Delta^{++}_{\nu}(x)\\!\int\\!\differential[4]{a}\,F(a)\,\partial_{\mu}\pi^{-}(x+a)+{\rm
H.c.}$ $\displaystyle+\
\frac{i}{4f_{\pi}^{2}}\,\bar{p}(x)\gamma^{\mu}p(x)\int\\!\differential[4]{a}\\!\int\\!\differential[4]{b}\
F(a)\,F(b)\,\left[\pi^{+}(x+a)\,\partial_{\mu}\pi^{-}(x+b)-\partial_{\mu}\pi^{+}(x+a)\pi^{-}(x+b)\right]$
$\displaystyle+\
\frac{2i(b_{10}+b_{11})}{f_{\pi}^{2}}\,\bar{p}(x)\sigma^{\mu\nu}p(x)\int\\!\differential[4]{a}\\!\int\\!\differential[4]{b}\
F(a)\,F(b)\,\partial_{\mu}\pi^{+}(x+a)\,\partial_{\nu}\pi^{-}(x+b).$

The function $F(a)$ here describes the strength of the correlation between the
baryon and the meson, and will become the regulator function in the momentum
space. In the limit when $F(a)\to\delta(a)$, one can show that Eq. (12)
reduces to the local version of the Lagrangian given in Eq. (10).

###  II.3 Nonzero skewness splitting functions

In this work we will be particularly concerned with the contributions to
proton GPDs from the pseudoscalar meson coupling diagrams in Fig. 1, as a
first step towards a complete one-loop calculation of the GPDs at nonzero
skewness. Contributions to GPDs from addition couplings to the octet and
decuplet baryons, including baryon rainbow and Kroll-Ruderman diagrams He  _et
al._ (2022), will be investigated in future work.

The contribution to the nonzero skewness splitting functions from the rainbow
diagram in Fig. 1(a) with the pseudoscalar meson and octet baryon is given by

$\displaystyle\bar{u}(p^{\prime})\,\Gamma_{(\rm a)}^{+}\,u(p)$
$\displaystyle=$
$\displaystyle\bar{u}(p^{\prime})\,\frac{C^{2}_{B\phi}}{f^{2}}\int\\!\frac{\differential[4]{k}}{(2\pi)^{4}}(\not{k}+\not{\Delta})\gamma_{5}\,\widetilde{F}(k+\Delta)\frac{i}{D_{\phi}(k+\Delta)}\,(2k^{+}+\Delta^{+})$
(13) $\displaystyle\hskip
65.44142pt\times\frac{i}{D_{\phi}(k)}\frac{i(\not{p}-\not{k}+M_{B})}{D_{B}(p-k)}\gamma_{5}\not{k}\,\widetilde{F}(k)\,\delta\Big{(}y+\xi-\frac{k^{+}}{P^{+}}\Big{)}\,u(p)$
$\displaystyle\equiv$
$\displaystyle\bar{u}(p^{\prime})\bigg{[}\gamma^{+}f_{\phi B}^{({\rm
rbw})}(y,\xi,t)+\frac{i\sigma^{+\nu}\Delta_{\nu}}{2M}g_{\phi B}^{({\rm
rbw})}(y,\xi,t)\bigg{]}u(p),$

where the propagator factors $D_{\phi}(k)$ and $D_{B}(p)$ are defined as

$D_{\phi}(k)=k^{2}-m_{\phi}^{2}+i\epsilon,\qquad\qquad
D_{B}(p)=p^{2}-M_{B}^{2}+i\epsilon,$ (14)

and $m_{\phi}$ and $M_{B}$ represent the pseudoscalar meson and baryon masses,
respectively. The function $\widetilde{F}$ is a regulator function, obtained
by performing a Fourier transformation on the factor $F(a)$ in Eq. (12), and
is used to regularize the ultraviolet divergence in the loop integration. For
simplicity we choose this to be a function of the meson momentum only (see
Sec. IV below). The explicit expressions for the splitting functions $f_{\phi
B}^{({\rm rbw})}$ and $g_{\phi B}^{({\rm rbw})}$ and the couplings
$C^{2}_{B\phi}$ are given in Appendix B.

The bubble diagrams are illustrated in Figs. 1(b) and 1(c). For the regular
bubble diagram in Fig. 1(b), the corresponding splitting function is given by

$\displaystyle\bar{u}(p^{\prime})\,\Gamma_{(\rm b)}^{+}\,u(p)$
$\displaystyle=$
$\displaystyle-\bar{u}(p^{\prime})\,\frac{C_{\phi\phi}}{2f^{2}}\int\\!\frac{\differential[4]{k}}{(2\pi)^{4}}\,i(2\not{k}+\not{\Delta})\,\widetilde{F}(k+\Delta)\widetilde{F}(k)\,\frac{i}{D_{\phi}(k+\Delta)}(2k^{+}+\Delta^{+})$
(15) $\displaystyle\hskip
65.44142pt\times\frac{i}{D_{\phi}(k)}\,\delta\Big{(}y+\xi-\frac{k^{+}}{P^{+}}\Big{)}\,u(p),$
$\displaystyle\equiv$
$\displaystyle\bar{u}(p^{\prime})\gamma^{+}u(p)\,f_{\phi}^{\rm(bub)}(y,\xi,t).$

The additional bubble diagram in Fig. 1(c) is derived from the last
interaction term in Eq. (10), and its contribution to the matrix element can
be written as

$\displaystyle\bar{u}(p^{\prime})\,\Gamma_{(\rm c)}^{+}\,u(p)$
$\displaystyle=$
$\displaystyle-\bar{u}(p^{\prime})\,\frac{C^{\prime}_{\phi\phi}}{f^{2}}\int\\!\frac{\differential[4]{k}}{(2\pi)^{4}}\,\sigma^{\lambda\nu}\Delta_{\lambda}k_{\nu}\,\widetilde{F}(k+\Delta)\widetilde{F}(k)\frac{i}{D_{\phi}(k+\Delta)}(2k^{+}+\Delta^{+})$
(16) $\displaystyle\hskip
65.44142pt\times\frac{i}{D_{\phi}(k)}\,\delta\Big{(}y+\xi-\frac{k^{+}}{P^{+}}\Big{)}\,u(p),$
$\displaystyle\equiv$ $\displaystyle\bar{u}(p^{\prime})\gamma^{+}u(p)\
g_{\phi}^{\prime\rm(bub)}(y,\xi,t),$

following the convention of Ref. He  _et al._ (2022).

Finally, for the decuplet baryon rainbow diagram in Fig. 1(d), the
contribution to $\Gamma^{+}$ is given by a form similar to that for the octet
baryon diagram in Fig. 1(a),

$\displaystyle\bar{u}(p^{\prime})\,\Gamma_{(\rm d)}^{+}\,u(p)$
$\displaystyle=$
$\displaystyle-\bar{u}(p^{\prime})\,\frac{C_{T\phi}^{2}}{f^{2}}\int\\!\frac{\differential[4]{k}}{(2\pi)^{4}}\,(k+\Delta)_{\lambda}\,\Theta^{\lambda\alpha}\widetilde{F}(k+\Delta)\frac{i}{D_{\phi}(k+\Delta)}\,(2k^{+}+\Delta^{+})$
(17)
$\displaystyle\times\frac{i}{D_{\phi}(k)}\frac{i}{\not{p}-\not{k}-M_{T}}S_{\alpha\beta}(p-k)\,\Theta^{\beta\rho}k_{\rho}\,\widetilde{F}(k)\,\delta\Big{(}y+\xi-\frac{k^{+}}{P^{+}}\Big{)}\,u(p)$
$\displaystyle\equiv$
$\displaystyle\bar{u}(p^{\prime})\bigg{[}\gamma^{+}f_{\phi T}^{({\rm
rbw})}(y,\xi,t)+\frac{i\sigma^{+\nu}\Delta_{\nu}}{2M}g_{\phi T}^{({\rm
rbw})}(y,\xi,t)\bigg{]}u(p),$

where $M_{T}$ is the mass of baryon decuplet. Again, the expressions for the
integrals of the splitting functions $f_{\phi T}^{({\rm rbw})}$ and $g_{\phi
T}^{({\rm rbw})}$, along with the coefficients $C^{2}_{T\phi}$, as well as the
corresponding bubble distributions $f_{\phi}^{\rm(bub)}$ and
$g_{\phi}^{\prime\rm(bub)}$, are given in Appendix B.

##  III GPDs from pseudoscalar meson loops

The recent analysis in Ref. He  _et al._ (2022) derived formulas for quark and
antiquark GPDs at zero skewness in terms of convolutions of the nucleon to
pseudoscalar meson plus baryon splitting functions, and GPDs associated with
the mesons and baryons in the intermediate state. In this section, we discuss
the extension of this formalism to the case of nonzero skewness. This exercise
requires paying careful attention to the different regions of parton momentum
fraction $x$, meson light-cone momentum fraction $y$, and skewness parameter
$\xi$.

For the rainbow diagrams in Figs. 1(a) and 1(d), the relevant splitting
functions are nonzero when $y$ is in the region $-\xi\leq y\leq 1$, for
positive $\xi$ values. The splitting function in this region can be convoluted
with the pion GPD to obtain the quark GPD in the physical nucleon. Figure 2
illustrates the decomposition of the convolution formula for the rainbow
diagram in Fig. 1(a) into different subprocesses. The splitting functions for
Figs. 2(a), 2(b) and 2(d) are all located in the $y>\xi$ region, while the
quark GPD in the pion belongs to the quark DGLAP, ERBL and antiquark DGLAP
regions, respectively. After convoluting with the hadronic splitting function,
these contribute to the quark GPD of the proton in the quark DGLAP, ERBL, and
antiquark DGLAP regions, respectively. For the subprocess in Fig. 2(c) the
splitting function is located at $y\in[-\xi,\xi]$, and can be considered as
the meson–meson pair annihilation process, analogous to the DA-like dynamics
at the quark level.

By combining these various processes, the contribution to the skewness nonzero
GPDs from the diagram in Fig. 1(a) can be expressed in the convolution form as

$\displaystyle\int_{x}^{1}\frac{\differential{y}}{y}f_{\phi B}^{({\rm
rbw})}(y,\xi,t)\,H_{q/\phi}\Big{(}\frac{x}{y},\frac{\xi}{y},t\Big{)},\hskip
130.88284pt[\xi<x<y]$ (18a)
$\displaystyle\int_{\xi}^{1}\frac{\differential{y}}{y}f_{\phi B}^{({\rm
rbw})}(y,\xi,t)\,H_{q/\phi}\Big{(}\frac{x}{y},\frac{\xi}{y},t\Big{)},\hskip
130.88284pt[x<\xi<y]$ (18b)
$\displaystyle\int_{-\xi}^{\xi}\frac{\differential{y}}{2y}\,f_{\phi B}^{({\rm
rbw})}(y,\xi,t)\frac{1}{\pi}\int_{s_{0}}^{\infty}\\!\\!\differential{s}\frac{{\rm
Im}\Phi_{q/\phi}{\big{(}\frac{1}{2}(1\\!+\\!\frac{x}{\xi}),\frac{1}{2}(1\\!+\\!\frac{y}{\xi}),s\big{)}}}{s-t+i\epsilon},\hskip
14.79555pt[|x|,|y|<\xi]\\!\\!$ (18c)
$\displaystyle\int_{-x}^{1}\frac{\differential{y}}{y}f_{\phi B}^{({\rm
rbw})}(y,\xi,t)\,H_{q/\phi}\Big{(}\frac{x}{y},\frac{\xi}{y},t\Big{)},\hskip
102.43008pt[\xi<-x<y<1]$ (18d)

where $H_{q/\phi}$ and $\Phi_{q/\phi}$ represent the valence GPD and
generalized distribution amplitude (GDA), respectively, in the intermediate
pseudoscalar meson, and $f_{\phi B}^{({\rm rbw})}$ is the rainbow splitting
function defined in Eq. (13). The corresponding convolution formula for the
quark magnetic GPD $E_{q}(x,\xi,t)$ can be obtained by replacing $f_{\phi
B}^{({\rm rbw})}(y,\xi,t)$ by the Pauli splitting function $g_{\phi B}^{({\rm
rbw})}(y,\xi,t)$. Similarly, the contribution from the pseudoscalar
meson–decuplet diagram in Fig. 1(d) can be obtained by replacing the octet
splitting functions $\\{f,g\\}_{\phi B}^{(\text{rbw})}(y,\xi,t)$ by the
corresponding decuplet functions $\\{f,g\\}_{\phi T}^{({\rm rbw})}(y,\xi,t)$.

The splitting functions appearing in Eqs. (18a), (18b) and (18d) are located
in the region $\xi<y$, and represent the contribution to the quark GPD in the
DGLAP region, the quark GPD in the ERBL region, and the antiquark GPD in the
DGLAP region, respectively. The convolution in Eq. (18c) also contributes to
the GPD in the ERBL region from the splitting function at $y<|\xi|$, with the
splitting function convoluted with the GDA of the meson, defined in the
timelike region. The results in the spacelike region can be obtained by
analytically continuing the intergral over $s$ in Eq. (18c), where $s_{0}\sim
2m_{\phi}$ is the threshold for the meson pair production. The total
contribution of the rainbow diagrams in Figs. 1(a) and 1(d) to the GPDs are
then obtained by summing the four subprocesses in Fig. 2.

Figure 2: Representation of the convolution formula in Eqs. (18a), (18b),
(18c), and (18d), respectively, with the {dashed, thick solid, thin solid}
lines representing the {pseudoscalar meson, proton, quark}. The processes in
diagrams (a) and (d) represent the DGLAP region for the quark and antiquark,
respectively, while the processes in (b) and (c) contribute to the ERBL
region.

For the bubble diagrams in Figs. 1(b) and 1(c), the splitting function is
nonzero only when the meson momentum fraction $y$ is in the range
$y\in[-\xi,\xi]$, so that these only contribute to GPDs in the ERBL region.
Their contributions to the quark electric GPD in the proton can be written as

$\displaystyle H_{q}^{(\rm bub)}(x,\xi,t)$ $\displaystyle=$
$\displaystyle\int_{-\xi}^{\xi}\frac{\differential{y}}{2y}\,f_{\phi}^{(\rm
bub)}(y,\xi,t)\,\frac{1}{\pi}\int_{s_{0}}^{\infty}\differential{s}\frac{\text{Im}\Phi_{q/\phi}{\big{(}\frac{1}{2}(1\\!+\\!\frac{x}{\xi}),\frac{1}{2}(1\\!+\\!\frac{y}{\xi}),s\big{)}}}{s-t+i\epsilon},$
(19)

where the splitting function for the bubble diagram, $f_{\phi}^{(\rm bub)}$,
is given in Eq. (15). The analogous expression for the magnetic GPD
$E_{q}^{(\rm bub)}(x,\xi,t)$ is obtained from Eq. (19) by replacing the bubble
splitting function $f_{\phi}^{(\rm bub)}$ by the corresponding function
$g_{\phi}^{\prime\rm(bub)}$ introduced in Eq. (16).

Integrating the GPDs over $x$ from $-1$ to $1$, we can write the lowest moment
of $H_{q}^{(\rm rbw)}$ in terms of the pseudoscalar meson elastic form factor.
Adding the contributions from Eqs. (18a), (18b) and (18d), integrated over the
intervals $[\xi,1]$, $[-\xi,\xi]$ and $[-1,-\xi]$, respectively, we have

$\displaystyle\\!\int_{\xi}^{1}\\!\differential{y}f_{\phi B}^{(\rm
rbw)}(y,\xi,t)\Bigg{[}\int_{\xi/y}^{1}\\!\differential{z}H_{q/\phi}\Big{(}z,\frac{\xi}{y},t\Big{)}+\int_{-1}^{-\xi/y}\\!\\!\\!\differential{z}H_{q/\phi}\Big{(}z,\frac{\xi}{y},t\Big{)}+\int_{-\xi/y}^{\xi/y}\\!\\!\differential{z}H_{q/\phi}\Big{(}z,\frac{\xi}{y},t\Big{)}\Bigg{]}$
$\displaystyle\hskip 99.58464pt=\int_{\xi}^{1}\differential{y}f_{\phi B}^{(\rm
rbw)}(y,\xi,t)\,F_{\phi}(t),$ (20)

where we have used the relation between the lowest moment of the quark GPD in
the pseudoscalar meson and the meson form factor. For the ERBL region in Eq.
(18c), the contribution to the integral of the GPD can be written as

$\displaystyle\int_{-\xi}^{\xi}\differential{x}H_{q}^{(\rm rbw)}(x,\xi,t)$
$\displaystyle=$ $\displaystyle\int_{-\xi}^{\xi}\differential{y}\,f_{\phi
B}^{(\rm rbw)}(y,\xi,t)\frac{\xi}{\pi
y}\int_{-1}^{1}\frac{\differential{z}}{2}\int_{s_{0}}^{\infty}\differential{s}\frac{\text{Im}\Phi_{q/\phi}\big{(}\frac{1}{2}(1+z),\frac{1}{2}(1+\frac{y}{\xi}),s\big{)}}{s-t+i\epsilon}$
$\displaystyle=$ $\displaystyle\int_{-\xi}^{\xi}\differential{y}\,f_{\phi
B}^{(\rm rbw)}(y,\xi,t)\frac{\xi}{\pi
y}\int_{s_{0}}^{\infty}\differential{s}\frac{\text{Im}\int_{0}^{1}\differential{\eta}\,\Phi_{q/\phi}({\eta,\frac{1}{2}(1+\frac{y}{\xi}),s)}}{s-t+i\epsilon}$
$\displaystyle=$ $\displaystyle\int_{-\xi}^{\xi}\differential{y}\,f_{\phi
B}^{(\rm
rbw)}(y,\xi,t)\frac{1}{\pi}\int_{s_{0}}^{\infty}\differential{s}\frac{\text{Im}F_{\phi}(s)}{s-t+i\epsilon}$
$\displaystyle=$ $\displaystyle\int_{-\xi}^{\xi}\differential{y}\,f_{\phi
B}^{(\rm rbw)}(y,\xi,t)\,F_{\phi}(t),$ (21)

where we have used the relation
$\int_{0}^{1}\differential{\eta}\Phi_{q/\phi}(\eta,\kappa,s)=(2\kappa-1)F_{\phi}(s)$.
Combining Eqs. (III) and (III), one then obtains the lowest moment of the GPD,

$\displaystyle\int_{-\xi}^{1}\differential{x}\,H_{q}^{(\rm rbw)}(x,\xi,t)$
$\displaystyle=$ $\displaystyle\int_{-\xi}^{1}\differential{y}\,f_{\phi
B}^{(\rm rbw)}(y,\xi,t)\,F_{\phi}(t),$ (22)

which corresponds to the Dirac form factor according to Eqs. (8) and (9).

In the numerical calculations in this analysis, we consider specifically the
case of the pion, $\phi=\pi$. For the pion GPDs we use the parametrization
from Amrath et al. Amrath  _et al._ (2008),

$\displaystyle H_{q/\pi}(x,\xi,t)$ $\displaystyle=$
$\displaystyle\int_{-1}^{1}\differential{\beta}\int_{-1+|\beta|}^{1-|\beta|}\differential{\alpha}\,\delta(x-\beta-\xi\alpha)\,h_{b}(\beta,\alpha)\,H_{q/\pi}(\beta,0,t)$
(23) $\displaystyle+$
$\displaystyle\frac{\xi}{|\xi|}D_{q/\pi}\Big{(}\frac{x}{\xi},t\Big{)}\,\theta(\xi-|x|),$

where the function $h_{b}$ is defined by

$\displaystyle
h_{b}(\beta,\alpha)=\frac{\Gamma(2b+2)}{2^{2b+1}\Gamma^{2}(b+1)}\frac{\big{[}(1-|\beta|)^{2}-\alpha^{2}\big{]}^{b}}{(1-|\beta|)^{2b+1}}.$
(24)

We take the parameter $b=2$, and the GPD in the zero skewness case is defined
according to the ansatz,
$H_{q/\pi}(\beta,0,t)=H_{q/\pi}(\beta,0,0)\,F_{\pi}(t)$. The second term in
Eq. (23) represents the contribution to the $D$-term of the gravitational form
factor of the pion, which is related to the pressure distribution of quarks
and gluons in the hadron. In order to obtain the correct first moment, this
term must be an odd function of $x$. In the present analysis, we use for the
$D$-term the simple form,

$\displaystyle D_{q/\pi}(z,t)=\frac{15}{4}z(1-z^{2})\,D_{q/\pi}(t),$ (25)

where $D_{q/\pi}(t)$ is the gravitational form factor of the pion, for which
we choose a monopole form,
$D_{q/\pi}(t)=D_{q/\pi}(0)/(1-t/\Lambda_{\pi}^{2})$, with normalization
$D_{q/\pi}(0)=-0.157$ and cutoff paramaeter $\Lambda_{\pi}=1.44$ GeV from the
recent lattice QCD calculation by Hackett et al. Hackett  _et al._ (2023).

In the forward limit the quark GPD in the pion can be related to the quark
PDF, $H_{q/\pi}(x,0,0)=q_{\pi}(x)$ for $x>0$, and
$H_{q/\pi}(x,0,0)=-\bar{q}_{\pi}(-x)$ for $x<0$ by crossing symmetry. Since we
only consider the lowest Fock state of the pion, composed of a valence quark
and antiquark, we set $q_{\pi}(x)=q^{v}_{\pi}(x)$, the valence quark PDF in
the pion, and $\bar{q}_{\pi}(x)=0$ for the sea quarks in the pion. Recent
analyses of pion PDFs have been performed by several groups Barry  _et al._
(2018); Novikov  _et al._ (2020); Barry  _et al._ (2021); Kotz  _et al._
(2023); for simplicity we use the older parametrization from Ref. Aicher  _et
al._ (2010).

For the contribution to the pion GPD from the ERBL region, we need in addition
information about the pion GDA, $\Phi_{q/\phi}$ [Eq. (18c)]. For this we also
use the double parametrization of the form Teryaev (2001); Kivel and Polyakov
(2002); Diehl (2003)

$\displaystyle\Phi_{q/\pi}(\eta,z,s)$ $\displaystyle=$ $\displaystyle
2(2z-1)\int_{-1}^{1}\differential{\beta}\int_{-1+|\beta|}^{1-|\beta|}\\!\\!\differential{\alpha}\delta\big{(}(2\eta-1)-(2z-1)\beta-\alpha)\big{)}\,h_{b}(\beta,\alpha)\,H_{q/\pi}(\beta,0,s)$
(26) $\displaystyle+$ $\displaystyle 2D_{q/\pi}(2\eta-1,s).$

Such a parametrization can be shown to satisfy the polynomial property in Eq.
(41).

As discussed above, in this calculation we only consider contributions to the
GPDs from direct couplings to the virtual pion loop, as in Fig. 1. Using quark
flavor symmetry, the pion GPDs and pion GDAs for the different possible charge
states can be related according to

$\displaystyle H_{u/\pi^{+}}$ $\displaystyle=$ $\displaystyle
H_{\bar{d}/\pi^{+}}=H_{d/\pi^{-}}=H_{\bar{u}/\pi^{-}}=2\,H_{q/\pi^{0}},$ (27)

and

$\displaystyle\Phi_{u/\pi^{+}}$ $\displaystyle=$
$\displaystyle\Phi_{\bar{d}/\pi^{+}}=\Phi_{d/\pi^{-}}=\Phi_{\bar{u}/\pi^{-}}=2\,\Phi_{q/\pi^{0}},$
(28)

respectively, where $q=u,\bar{u},d$ or $\bar{d}$ for the quark flavors in the
neutral $\pi^{0}$. Combining the contributions from the different diagrams in
Fig. 1, including rainbow and bubble, and the different isospin channels, the
$u$-quark GPD in the proton can be written as

$\displaystyle H_{u}(x,\xi,t)$ $\displaystyle=$ $\displaystyle H^{({\rm
rbw})\pi^{+}n}_{u}(x,\xi,t)+H^{({\rm
rbw})\pi^{+}\Delta^{0}}_{u}(x,\xi,t)+H^{({\rm
rbw})\pi^{0}p}_{u}(x,\xi,t)+H^{({\rm rbw})\pi^{0}\Delta^{+}}_{u}(x,\xi,t)$
(29) $\displaystyle+$ $\displaystyle H^{({\rm
bub})\pi^{+}\pi^{-}}_{u}(x,\xi,t)$ $\displaystyle-$ $\displaystyle H^{({\rm
rbw})\pi^{0}p}_{\bar{u}}(-x,\xi,t)-H^{({\rm
rbw})\pi^{0}\Delta^{+}}_{\bar{u}}(-x,\xi,t)-H^{({\rm
rbw})\pi^{-}\Delta^{++}}_{\bar{u}}(-x,\xi,t),$

where we have explicitly written out the specific pion-baryon charge
contributions, and similarly for the $u$-quark magnetic GPD $E_{u}$. The
results of $d$-quark GPDs can be obtained using the isospin symmetry
$\\{H,E\\}_{d}(x,\xi,t)=-\\{H,E\\}_{u}(-x,\xi,t)$, since we only consider the
pion loop contribution in this work. Note that Eq. (29) includes contributions
from the quark DGLAP, ERBL, and antiquark DGLAP regions. As mentioned above,
contributions from the baryon coupling rainbow diagrams, as well as Kroll-
Ruderman diagrams He  _et al._ (2022), will be investigated in future work Gao
_et al._ (2024).

##  IV Phenomenology of meson loop contributions to GPDs

In this section we will present the phenomenological results for the
generalized splitting functions derived in the preceeding sections, along with
the meson loop contributions to the light quark GPDs obtained using the
convolution formula in Eq. (18). We also apply the formalism to calculation of
the Dirac, Pauli, and gravitational form factors, obtained from the lowest two
moments of the GPDs. In the numerical calculations the loop integrals are
regularized using a covariant regulator of a dipole form,

$\widetilde{F}(k)=\left(\frac{\Lambda^{2}-m^{2}_{\phi}}{\Lambda^{2}-k^{2}}\right)^{2},$
(30)

with $\Lambda$ a mass parameter. In practice, we use the value
$\Lambda=1.0(1)$ GeV obtained from previous analyses of meson loop
contributions to PDFs and GPDs He and Wang (2018a, b); He  _et al._ (2022).

###  IV.1 Splitting functions

The 3-dimensional splitting functions $y\,f_{i}(y,\xi,t)$ for the rainbow and
bubble diagrams are shown in Fig. 3, where for the purpose of comparison we
focus on the case of a charged $\pi^{+}$ loop. For illustration, we fix
$\xi=0.1$, and display the splitting functions for $y$ between $-\xi$ and 1,
and for $-t<1$ GeV2. As expected, the absolute values of these splitting
functions decrease with increasing $-t$. For the octet baryon rainbow diagram,
the splitting function $y\,f_{\pi^{+}n}^{\rm(rbw)}$ has a similar peak value
within the regions $y>\xi$ and $|y|<\xi$. However, the absolute value of the
Pauli splitting function $g_{\pi^{+}n}^{\rm(rbw)}$ in the $|y|<\xi$ region is
smaller than that in the $y>\xi$ region. Both $f_{\pi^{+}n}^{\rm(rbw)}$ and
$g_{\pi^{+}n}^{\rm(rbw)}$ are continuous at the boundary $y=\xi$. The shapes
of these are similar in the $y>\xi$ region; however, the Dirac splitting
function’s shape is sharper when $y$ ranges from $-\xi$ to $\xi$.

For the decuplet baryon rainbow diagram, the shape of the splitting function
$y\,g_{\pi^{+}\Delta^{0}}^{\rm(rbw)}$ is similar as that of
$y\,f_{\pi^{+}n}^{\rm(rbw)}$, although their signs in the $y>\xi$ region are
opposite. The Dirac-like splitting function for the decuplet rainbow
$f_{\pi^{+}\Delta^{0}}^{\rm(rbw)}$ decreases rapidly near the endpoints
$y=\xi$ and $y=-\xi$, and is not continuous at the boundary $y=\xi$. Such a
discontinuity arises from terms with higher powers of the meson momentum in
the loop integral for the splitting function, for example, the $(k\cdot
p)^{2}k\cdot p^{\prime}$ terms in Eq. (55a). For the bubble diagrams in Figs.
1(b) and 1(c), the relevant splitting functions are nonzero only when $y$
satisfies $|y|<\xi$, and $y\,f_{\pi^{+}}^{\rm(bub)}$ and
$y\,g_{\pi^{+}}^{\prime{\rm(bub)}}$ are antisymmetric about the $y=0$ axis.
The absolute value of the Pauli-like $g_{\pi^{+}}^{\prime\rm(bub)}$ splitting
function is approximately an order of magnitude larger than that of the Dirac-
like function $f_{\pi^{+}}^{\rm(bub)}(y)$.

Figure 3: Three-dimensional splitting functions $y\,f_{i}(y,\xi,t)$ as a
function of $y$ and $t$, with fixed $\xi=0.1$, for the octet baryon rainbow
(upper panels), decuplet baryon rainbow (middle panels), bubble and magnetic
bubble (lower panels) diagrams.

To more clearly illustrate the shapes of the splitting functions, we show the
two-dimensional projections of the functions in Fig. 4 for $-t=1$ GeV2 and
$\xi=0.1$. The splitting function $y\,f_{\pi^{+}n}^{\rm(rbw)}$ is
approximately antisymmetric for $-\xi<y<\xi$. Two sharp humps are seen when
$y$ is close to the $\pm\xi$ points. For the splitting function
$g_{\pi^{+}n}^{\rm(rbw)}$, two humps are also seen when $-\xi<y<\xi$, although
the absolute values of the maxima of the humps are smaller than those in the
$\xi<y<1$ region.

For the decuplet baryon rainbow diagram, the two-dimensional splitting
function $g_{\pi^{+}\Delta^{0}}^{\rm(rbw)}$ has a similar shape to that of
$f_{\pi^{+}n}^{\rm(rbw)}$. The shape of $f_{\pi^{+}\Delta^{0}}^{\rm(rbw)}$ is
different from that of the other functions when $-\xi<y<\xi$, with values near
the $y=-\xi$ and $y=\xi$ boundary decreasing rapidly, and the shape is neither
symmetric nor antisymmetric about the $y=0$ axis. Such a behavior is in fact
due to the effect from specific terms, which behave as $\delta(y)/y$ when
$\xi$ tends to zero. The splitting functions for the bubble diagrams,
$y\,f_{\pi^{+}}^{\rm(bub)}$ and $y\,g_{\pi^{+}}^{\prime\rm(bub)}$, are
antisymmetric about the $y=0$ axis, so that by dividing by $y$ the resulting
functions $f_{\pi^{+}}^{\rm(bub)}$ and $g_{\pi^{+}}^{\prime\rm(bub)}$ will be
symmetric.

Figure 4: Two-dimensional projections of the Dirac-like $yf_{i}(y,\xi,t)$ and
Pauli-like $yg_{i}(y,\xi,t)$ splitting functions for the octet baryon rainbow
(top panels), decuplet baryon rainbow (middle panels), and bubble and magnetic
bubble diagrams (lower panels), for $\xi=0.1$ and $-t=1$ GeV2. The dashed
vertical lines in the rainbow diagram plots correspond to the point $y=\xi$.

In Fig. 5, we show the results for the splitting functions
$f_{\pi^{+}n}^{\rm(rbw)}$, $f_{\pi^{+}\Delta^{0}}^{\rm(rbw)}$,
$f_{\pi^{+}}^{\rm(bub)}$, and $g_{\pi^{+}}^{\prime\rm(bub)}$ when $y$ is in
the region $[-\xi,\xi]$, for $\xi$ values between 0.001 and 0.1. To more
conveniently compare the results with different $\xi$, we scale the horizontal
axis by dividing $y$ by $\xi$, and the vertical axis by multiplying the
splitting functions by $\xi$. The scaled results with different $\xi$ values
are very similar and converge as $\xi$ becomes small, with the results for
$\xi<0.001$ almost indistinguishable from those at $\xi=0.001$. This behavior
is reminiscent of the behavior of a $\delta$-function, as has been discussed
in previous chiral loop calculations Salamu  _et al._ (2015, 2019a). In
contrast, for the decuplet baryon splitting function one observes convergence
at small $\xi$ after multiplying $f_{\pi^{+}\Delta^{0}}$ by $\xi^{2}$. Since
it is antisymmetric in $y$, this indicates that the splitting function will
approach $\delta(y)/y$ rather than $\delta(y)$ in the zero-skewness limit.

Figure 5: Scaled splitting functions $f_{i}(y,\xi,t)$ as a function of $y/\xi$
in the $-\xi<y<\xi$ region for $\xi=0.001$ (red lines) and $\xi=0.1$ (blue
lines), for $-t=1$ GeV2. All the functions are scaled by a factor of $\xi$,
except for the decuplet function $f_{\pi^{+}\Delta^{0}}^{\rm(rbw)}(y,\xi,t)$,
which is scaled by a factor $\xi^{2}$.

###  IV.2 Pion loop contributions to quark GPDs

Performing the convolutions of the splitting functions with the corresponding
quark GPDs of the virtual pions, we show the contributions of the chiral loop
to the $H$ and $E$ GPDs for $u$ quarks in Fig. 6 as a function of $x$ and $t$
at a fixed $\xi=0.1$. For the contributions from the couplings to the virtual
pion loops, as in Fig. 1, the results for the $d$-quark GPDs can be obtained
from the $u$-quark distributions using isospin symmetry,
$\\{H,E\\}_{d}(x,\xi,t)=-\\{H,E\\}_{u}(-x,\xi,t)$. The most striking structure
is seen in the region of low $x$, $|x|\lesssim 0.2$. For the case of the
electric $u$-quark GPD, $xH_{u}$ is positive in the DGLAP region, and has two
valleys in the ERBL region. For the magnetic GPD, $xE_{u}$ is negative in the
$x<0$ region and positive when $x>0$. The distributions fall rapidly as
$|x|\to 1$, and for increasing values of the momentum transfer squared, $-t$.

Figure 6: Total three-dimensional $u$-quark GPDs $xH_{u}$ and $xE_{u}$ from
Eq. (29) as a function of $x$ and $t$, for fixed $\xi=0.1$. The corresponding
$d$-quark distributions can be obtained from the $u$-quark GPDs using the
isospin symmetery relation, $\\{H,E\\}_{d}(x,\xi,t)=-\\{H,E\\}_{u}(-x,\xi,t)$,
which holds for the contributions from virtual pion loops as in Fig. 1.

Figure 7: Meson loop contributions to the GPDs $H_{u}$ and $E_{u}$ for
$\xi=0.1$ and $-t=1$ GeV2. The corresponding $d$-quark GPDs can be obtained
from the relation $\\{H,E\\}_{d}(x,\xi,t)=-\\{H,E\\}_{u}(-x,\xi,t)$.

To more clearly illustrate the dependence of the GPDs on the parton momentum
fraction $x$, in Fig. 7 we plot the two-dimensional projections of the GPDs
$xH_{u}$ and $xE_{u}$ for $-t=0.25$ GeV2 and $-t=1$ GeV2. Although the
splitting function $f_{\pi^{+}\Delta^{0}}^{\rm(rbw)}$ is discontinuous at
$y=\xi$, as evident from Fig. 4, the quark GPDs remain continuous at
$x=\pm\xi$. The reason is that the convolution formula in Eqs. (18a), (18b)
and (18d) are only related to the splitting function in the $y>\xi$ region,
and the results obtained using these should be continuous at $\xi$ if the
input pion GPD is continuous at $x=\xi$. On the other hand, the contribution
to the quark GPD from Eq. (18c) vanishes $x=\xi$ because of the endpoint
property of the distribution amplitude, $\Phi_{q/\pi}(1,\kappa,s)=0$.

Comparing the results with different $-t$, the absolute values of the GPDs at
$-t=0.25$ GeV2 are around 4–5 times larger than those at $-t=1$ GeV2. For the
$u$-quark distribution, one finds that the quark GPD $H_{u}$ in the DGLAP
region for $x>0$ has a larger magnitude than $H_{u}$ at $x<0$ (which by
crossing symmetry is equivalent to the $\bar{u}$ distribution at $x>0$).
According to the above isospin relation for the contributions from the pion
coupling diagrams, $xH_{u}(x>0)$ is identical to $xH_{d}(x<0)$, the latter
which is equivalent to the $x\bar{d}$ distribution at $x>0$. Therefore the
contributions from the loop diagrams in Fig. 1 naturally give an enhancement
of the $\bar{d}$ distribution compared with the $\bar{u}$, reminiscent of the
empirical result for the $\bar{d}$ and $\bar{u}$ PDF asymmetry in the
collinear region Salamu  _et al._ (2015).

Moreover, the result in the $0<x<\xi$ region is positive for $-t=0.25$ GeV2,
but negative for $-t=1$ GeV2, which can be understood from the fact that the
pion GPD or GDA includes the $D_{q/\pi}(t)$ form factor in the ERBL region
that has a different $t$ dependence compared with the pion form factor
$F_{\pi}(t)$. The magnitude of the Pauli GPD $E_{u}$ is larger than that of
the Dirac GPD $H_{u}$, and the quark and antiquark distributions in the DGLAP
region have opposite signs. This implies that the GPD for the $\bar{u}$ flavor
has a different sign to that of the Dirac GPD for $\bar{d}$, and the absolute
value of the Pauli GPD $E$ for $\bar{d}$ is larger than that for $\bar{u}$.

###  IV.3 Form factors

From the calculated $H_{q}$ and $E_{q}$ GPDs one can compute the Dirac, Pauli,
and gravitational form factors from the lowest two moments of the GPDs. The
Dirac and Pauli form factors can be obtained from Eq. (6a) and (6b),
respectively. For the gravitational form factors, the pion loop contributions
to $A^{q}$, $B^{q}$ and $C^{q}$ can be written as

$\displaystyle A^{q}(t)$ $\displaystyle=$
$\displaystyle\sum_{i}A^{2(0)}_{q/\pi}(t)\,{\cal A}^{2(0)}_{i}(t),$ (31a)
$\displaystyle B^{q}(t)$ $\displaystyle=$
$\displaystyle\sum_{i}A^{2(0)}_{q/\pi}(t)\,{\cal B}^{2(0)}_{i}(t),$ (31b)
$\displaystyle C^{q}(t)$ $\displaystyle=$
$\displaystyle\sum_{i}\Big{[}A^{2(0)}_{q/\pi}(t)\,{\cal
C}^{(2)}_{i}(t)+C^{(2)}_{q/\pi}(t)\,{\cal C}^{\prime}_{i}(t)\Big{]},$ (31c)

where $A^{2(0)}_{q/\pi}$ and $C^{(2)}_{q/\pi}$ are the quark gravitational
form factors in the pion, obtained from the pion GPD in Eq. (23) through

$\int_{-1}^{1}\differential{x}x\,H_{q/\pi}(x,\xi,t)=A^{2(0)}_{q/\pi}(t)+(2\xi)^{2}\,C^{(2)}_{q/\pi}(t).$
(32)

Note that the $C^{(2)}_{q/\pi}$ form factor is proportional to the pion
gravitational form factor $D_{q/\pi}(t)$ in Eq. (25),
$C^{(2)}_{q/\pi}(t)=\frac{1}{4}D_{q/\pi}(t)$. The sum over the index $i$ in
Eqs. (31) includes the various hadronic configurations in Fig. 1, and the
corresponding hadronic form factors ${\cal A}^{2(0)}_{i}$, ${\cal
B}^{2(0)}_{i}$ and ${\cal C}_{i}^{(2)}$ are obtained from the first moments of
splitting functions,

$\displaystyle\int_{-\xi}^{1}\differential{y}\,y\,f_{i}(y,\xi,t)$
$\displaystyle=$ $\displaystyle{\cal A}^{2(0)}_{i}(t)+(2\xi)^{2}\,{\cal
C}^{(2)}_{i}(t),$ (33a)
$\displaystyle\int_{-\xi}^{1}\differential{y}\,y\,g_{i}(y,\xi,t)$
$\displaystyle=$ $\displaystyle{\cal B}^{2(0)}_{i}(t)-(2\xi)^{2}\,{\cal
C}^{(2)}_{i}(t),$ (33b)

with the form factor ${\cal C}^{\prime}_{i}(t)$ for the rainbow diagrams in
Figs. 1(a) and 1(d) obtained from

$\displaystyle{\cal C}^{\prime}_{i}(t)$ $\displaystyle=$
$\displaystyle\int_{-\xi}^{1}\frac{\differential{y}}{y}\,\,f_{i}(y,\xi,t)=-\int_{-\xi}^{1}\frac{\differential{y}}{y}\,\,g_{i}(y,\xi,t).$
(34)

The ${\cal C}^{\prime}_{i}(t)$ form factor for Figs. 1(b) and 1(c) can also be
obtained from Eq. (34) by changing the integral region to $y\in[-\xi,\xi]$.
However, since for these diagrams the functions
$\frac{1}{y}f_{\pi^{+}}^{\rm(bub)}$ and
$\frac{1}{y}g_{\pi^{+}}^{\prime\rm(bub)}$ are antisymmetric with respect to
$y=0$, the integrals in Eq. (34), and hence ${\cal C}^{\prime}_{i}(t)$,
vanish. (Note that the functions $yf_{\pi^{+}}^{\rm(bub)}$ and
$yg_{\pi^{+}}^{\prime\rm(bub)}$ shown in Fig. 4 are also antisymmetric with
respect to $y$.)

Figure 8: Gravitational form factors $A^{q}$ (red band), $B^{q}$ (blue band),
and $C^{q}$ (green band) for quark flavor $q=u$ or $d$ as a function of $t$
arising from the pseudoscalar meson loop diagrams in Fig. 1. The results for
the $u$ and $d$ flavors for these diagrams are identical due to isospin
symmetry.

The numerical results for the gravitational $A^{q}$, $B^{q}$ and $C^{q}$ form
factors are shown in Fig. 8 as a function of $t$. Because of isospin symmetry
the results for the $q=u$ flavor are identical to those for the $q=d$ flavor.
Despite the absolute value of the GPD $E_{q}$ being larger than that of the
GPD $H_{q}$, the first moment of $E_{q}$ largely cancels between the
positive-$x$ and negative-$x$ regions, leading to similar values for the
$A^{q}$ and $B^{q}$ form factors. In contrast, the $C^{q}$ form factor is
negative and has a smaller magnitude compared with the other two form factors.

As pointed out in the preceeding discussions, in this analysis we only
consider contributions to GPDs arising from the meson loop diagrams in Fig. 1.
Generally, comparison with experimental measurements or other phenomenological
calculations of the total quark gravitational form factors will require
inclusion of contributions from coupling to intermediate state baryons in the
loops or Kroll-Ruderman type contributions, which will be considered in future
analyses. On the other hand, in a valence quark approximation the antiquark
distributions arise predominantly from couplings to the meson loop, which
allows their contributions to form factors and sea quark flavor asymmetries to
be studied directly in terms of the diagrams in Fig. 1 alone.

Figure 9: Meson loop contributions to the Dirac $F_{1}^{q}$ and Pauli
$F_{2}^{q}$ (upper panels) and gravitational form factors $A^{q}$ and $B^{q}$
(lower panel) versus $-t$ for $q=\bar{d}-\bar{u}$.

In Fig. 9 we illustrate the $\bar{d}-\bar{u}$ contributions to the Dirac
$F_{1}^{q}$, Pauli $F_{2}^{q}$, and gravitational $A^{q}$ and $B^{q}$ form
factors. For the pion loop coupling diagrams in Fig. 1, these are given in
terms of the moments of the $H_{u}$ and $E_{u}$ GPDs at zero skewness by

$\displaystyle F_{1}^{\bar{d}-\bar{u}}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\differential{x}\big{(}H_{u}(x,0,t)+H_{u}(-x,0,t)\big{)},$
(35a) $\displaystyle F_{2}^{\bar{d}-\bar{u}}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\differential{x}\big{(}E_{u}(x,0,t)+E_{u}(-x,0,t)\big{)},$
(35b) $\displaystyle A^{\bar{d}-\bar{u}}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\differential{x}\,x\big{(}H_{u}(x,0,t)+H_{u}(-x,0,t)\big{)},$
(35c) $\displaystyle B^{\bar{d}-\bar{u}}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\differential{x}\,x\big{(}E_{u}(x,0,t)+E_{u}(-x,0,t)\big{)},$
(35d)

where we have used the crossing symmetry and isospin relations discussed in
Sec. IV.2 to relate the GPDs for the $u$ and $d$ flavors. From the calculated
$H^{u}$ GPD, we find the $\bar{d}-\bar{u}$ contribution to the Dirac form
factor at $t=0$ is $F_{1}^{\bar{d}-\bar{u}}(0)=0.11(2)$. Note that
$F_{1}^{\bar{d}-\bar{u}}(0)$ corresponds to the light antiquark PDF asymmetry,

$\displaystyle F_{1}^{\bar{d}-\bar{u}}(0)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\differential{x}\big{(}\bar{d}(x)-\bar{u}(x)\big{)},$
(36)

which has been the subject of considerable interest in the literature
Ackerstaff  _et al._ (1998); Towell  _et al._ (2001); Dove  _et al._ (2021);
Salamu  _et al._ (2015). Our result is consistent with phenomenological
extractions from experimental Drell-Yan lepton-pair production data in $pp$
and $pd$ collisions,
$\int_{0.015}^{0.35}\differential{x}(\bar{d}-\bar{u})=0.0803(11)$ Towell  _et
al._ (2001), as well as from global QCD analyses of all high-energy scattering
data Cocuzza  _et al._ (2021); Accardi  _et al._ (2023). The Pauli form
factor, on the other hand, is an order of magnitude larger than the Dirac form
factor, $F_{2}^{\bar{d}-\bar{u}}(0)=1.15(25)$, which is consistent with
previous calculations of meson loop effects in chiral effective theory He  _et
al._ (2022). Although challenging to obtain directly from experimental
measurements, this contribution may be accessible in future lattice QCD
calculations.

For the $\bar{d}-\bar{u}$ contributions to the gravitational form factors, we
find in the forward limit $A^{\bar{d}-\bar{u}}(0)=0.01_{-0.002}^{+0.003}$ and
$B^{\bar{d}-\bar{u}}(0)=0.035(7)$, which is consistent with the previous
results from Ref. He  _et al._ (2022). Note that the $A^{\bar{d}-\bar{u}}$
form factor represents the $\bar{d}-\bar{u}$ asymmetry for the light antiquark
momentum fractions,

$\displaystyle A^{\bar{d}-\bar{u}}(0)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\differential{x}x\big{(}\bar{d}(x)-\bar{u}(x)\big{)}.$
(37)

Our calculated result can be compared with the truncated moment obtained by
the SeaQuest Collaboration from $pp$ and $pd$ Drell-Yan data at Fermilab Dove
_et al._ (2023),
$\int_{0.13}^{0.45}\differential{x}\,x\big{(}\bar{d}(x)-\bar{u}(x)\big{)}\approx
0.00318$. The numerical value for the form factor $B^{\bar{d}-\bar{u}}(0)$ is
about 3 times larger than the value of $A^{\bar{d}-\bar{u}}(0)$, and is also
consistent with the prediction in the large-$N_{c}$ limit Goeke  _et al._
(2001).

##  V Summary

Within the framework of nonlocal chiral effective theory, we have calculated
the one-loop contributions to the spin-averaged GPDs of the proton at nonzero
skewness $\xi$ from coupling to pseudoscalar meson loops with intermediate
octet and decuplet baryon states. We derived the generalization of convolution
formulas for nonzero skewness, which, unlike the zero skewness case, now
involve both the DGLAP and ERBL regions. We verified that the moments of the
generalized convolution formulas for the GPDs satisfy the correct polynominal
properties.

To regularize the ultraviolet divergences in the loop integrations for the
GPDs, we introduced a relativistic regulator derived self-consistently from
the nonlocal Lagrangian. Using constraints for this regulator obtained from
previous analyses of meson loop contributions to PDFs and GPDs He and Wang
(2018a, b); He  _et al._ (2022), in addition to phenomenological input for the
pion GPD Amrath  _et al._ (2008) and distribution amplitude Teryaev (2001);
Kivel and Polyakov (2002); Diehl (2003), we studied the detailed dependence of
the $H_{q}$ and $E_{q}$ GPDs for the $u$ and $d$ quarks on the parton momentum
fraction $x$, skewness $\xi$, and four-momentum transfer squared $t$. At the
kinematics considered in this analysis, we generally find that the absolute
values of the Pauli GPD $E_{q}$ are significantly larger than those of the
Dirac GPD $H_{q}$, with both suppressed for increasing $-t$.

For the $\bar{d}-\bar{u}$ flavor combination, we computed the contribution
from pion loops to the Dirac and Pauli electromagnetic and the gravitational
form factors as a function of $t$. At $t=0$ the normalization of the
$F_{1}^{\bar{d}-\bar{u}}$ form factor is equivalent to the $x$-integrated
value of the $\bar{d}-\bar{u}$ PDF asymmetry, and for the gravitational
$A^{\bar{d}-\bar{u}}(0)$ form factor the normalization corresponds to the
$x$-weighted moment of $\bar{d}-\bar{u}$. Our results are consistent with
phenomenological extractions of both asymmetries from experimental data Towell
_et al._ (2001); Cocuzza  _et al._ (2021), and with previous zero skewness
meson loop GPD calculations He  _et al._ (2022).

An obvious extension of the present work is the inclusion of contributions
from the coupling to intermediate state baryons, similar to the zero skewness
GPD analysis in Ref. He  _et al._ (2022). This will allow comparison with
quark as well as antiquark GPDs, as would be needed for studies of strange
quark contributions to proton GPDs, for instance Salamu  _et al._ (2019a, b),
or gravitational form factors. Such a strategy would also allow the
calculation of chiral loop contributions to spin-dependent GPDs, and even
generalized transverse momentum distributions, which will provide further
insights into the three-dimensional structure of the proton.

###### Acknowledgements.

We thank Yuxun Guo and Shuzhe Shi for helpful discussions. F.H. is supported
by the National Science Foundation under Award PHY-1847893. This work was
supported by NSFC under Grant No. 11975241, the DOE Contract No. DE-
AC05-06OR23177, under which Jefferson Science Associates, LLC operates
Jefferson Lab, DOE Contract No. DE-FG02-03ER41260, and in part within the
framework of the Quark-Gluon Tomography (QGT) Topical Collaboration, under
Contract No. DE-SC0023646.

##  Appendix A Proof of polynomiality of GPDs

Polynomiality requires that the $n$-th moments of GPDs can be expanded as a
series in $\xi$,

$\displaystyle\int_{-1}^{1}\differential{x}x^{n-1}\,H_{q}(x,\xi,t)=\sum_{i=0,\rm
even}^{n-1}\left(2\xi\right)^{i}A^{n(i)}_{q}(t)+(2\xi)^{n}C^{(n)}_{q}(t)\Big{|}_{n\,\rm
even},$ (38)

and similarly for the $E_{q}$. In analogy with the proof for the form factor
in Sec. III, the contributions from Eqs. (18a), (18b) and (18d) can be written
as

$\displaystyle\int_{-\xi}^{1}\differential{x}x^{n-1}\,H_{q}(x,\xi,t)$
$\displaystyle=$
$\displaystyle\int_{\xi}^{1}\differential{y}y^{n-1}\,f(y,\xi,t)\int_{-1}^{1}\differential{\tilde{z}}\,\tilde{z}^{n-1}\,H_{q/\pi}(\tilde{z},\xi/y,t)$
(39)
$\displaystyle=\int_{\xi}^{1}\differential{y}\,y^{n-1}\,f(y,\xi,t)\left(\sum_{i=0,\rm
even}^{n-1}\left(\frac{2\xi}{y}\right)^{i}A^{n(i)}_{q/\pi}(t)+\left(\frac{2\xi}{y}\right)^{n}C^{(n)}_{q/\pi}(t)\Big{|}_{n\,\rm
even}\right).$

The $n$-th moments of GPDs from Eq. (18c) can be expressed as

$\displaystyle\int_{-\xi}^{\xi}\differential{x}x^{n-1}\,H_{q}(x,\xi,t)$
$\displaystyle=$
$\displaystyle\int_{-\xi}^{\xi}\differential{y}\,\xi^{n-1}\,f(y,\xi,t)\frac{1}{\pi}\frac{\xi}{y}\int_{s_{0}}^{\infty}\differential{s}\frac{\text{Im}\int_{0}^{1}\differential{\eta}\,(2\eta-1)^{n-1}\,\Phi_{q/\pi}({\eta,\frac{\frac{y}{\xi}+1}{2},s)}}{s-t+i\epsilon}$
(40)
$\displaystyle=\int_{-\xi}^{\xi}\differential{y}\,\xi^{n-1}\,f(y,\xi,t)\left(\sum_{i=0,\rm
even}^{n-1}2^{i}\left(\frac{y}{\xi}\right)^{n-i-1}A^{n(i)}_{q/\pi}(t)+2^{n}\frac{\xi}{y}C^{(n)}_{q/\pi}(t)\Big{|}_{n\,\rm
even}\right),$

where we have used the property of GDAs,

$\displaystyle\int_{0}^{1}\differential{\eta}\,(2\eta-1)^{n-1}\,\Phi_{q/\pi}(\eta,z,s)$
$\displaystyle=$ $\displaystyle\frac{2^{n-1}}{(p^{+}+p^{\prime+})^{n}}\langle
0|\bar{q}\gamma^{+}(\mathop{i\partial^{+}}\limits^{\leftrightarrow})^{n-1}q(0)|\pi(p)\pi(p^{\prime})\rangle$
(41)
$\displaystyle=\frac{2^{n-1}}{(p^{\prime+}+p^{+})^{n}}\langle\pi(-p^{\prime})|\bar{q}\gamma^{+}(\mathop{i\partial^{+}}\limits^{\leftrightarrow})^{n}q(0)|\pi(p)\rangle$
$\displaystyle=\frac{2^{n-1}}{(p^{\prime+}+p^{+})^{n}}\bigg{(}\sum_{i=0,\rm
even}^{n-1}(p^{+}-p^{\prime+})(p^{+}+p^{\prime+})^{i}\left(\frac{p^{+}-p^{\prime+}}{2}\right)^{n-1-i}A_{q/\pi}^{n(i)}(s)$
$\displaystyle\hskip
56.9055pt+2\,(p^{+}+p^{\prime+})^{n}C^{(n)}_{q/\pi}(s)\Big{|}_{n\,\rm
even}\bigg{)}$ $\displaystyle=\sum_{i=0,\rm
even}^{n-1}2^{i}\left(\frac{p^{\prime+}-p^{+}}{p^{\prime+}+p^{+}}\right)^{n-i}A_{q/\pi}^{n(i)}(s)+2^{n}C^{(n)}_{q/\pi}(s)\Big{|}_{n\,\rm
even}$ $\displaystyle=\sum_{i=0,\rm
even}^{n-1}2^{i}\left(2z-1\right)^{n-i}A_{q/\pi}^{n(i)}(s)+2^{n}C^{(n)}_{q/\pi}(s)\Big{|}_{n\,\rm
even}.$

The parameter $z$ here is defined as $z=p^{\prime+}/(p^{+}+p^{\prime+})$ in
the GDA. Combining the results in Eq. (39) and (40), the $n$-th moments of
GPDs can be written

$\displaystyle\int_{\xi}^{1}\differential{x}x^{n-1}\,H_{q}(x,\xi,t)$
$\displaystyle=$
$\displaystyle\int_{-\xi}^{1}\differential{y}\,y^{n-1}\,f(y,\xi,t)\left(\sum_{i=0,\rm
even}^{n-1}\left(\frac{2\xi}{y}\right)^{i}A^{n(i)}_{q/\pi}(t)+\left(\frac{2\xi}{y}\right)^{n}C^{(n)}_{q/\pi}(t)\Big{|}_{n\,\rm
even}\right)$ (42) $\displaystyle=\sum_{i=0,\rm
even}^{n-1}(2\xi)^{i}A^{n(i)}_{q/\pi}(t)\int_{-\xi}^{1}\differential{y}y^{n-i-1}\,f(y,\xi,t)+(2\xi)^{n}C^{(n)}_{q/\pi}(t)\int_{-\xi}^{1}d\,y\,\,\frac{f(y,\xi,t)}{y}$
$\displaystyle=\sum_{i=0,\rm
even}^{n-1}(2\xi)^{i}A^{n(i)}_{q/\pi}(t)\left(\sum_{j=0,\rm
even}^{n-i-1}(2\xi)^{j}{\cal A}^{n(j)}(t)+\left(2\xi\right)^{n-i}{\cal
C}^{(n)}(t)\Big{|}_{n\,\rm even}\right)$
$\displaystyle+\,(2\xi)^{n}\,C^{(n)}_{q/\pi}(t\,){\cal C^{\prime}}(t),$

where ${\cal A}^{n(j)}$, ${\cal C}^{(n)}$, ${\cal C^{\prime}}$ are the moments
of splitting function, they are defined as

$\displaystyle\int_{-\xi}^{1}\differential{y}y^{n-1}\,f(y,\xi,t)$
$\displaystyle=$ $\displaystyle\sum_{j=0,\rm even}^{n-1}(2\xi)^{j}{\cal
A}^{n(j)}(t)+\left(2\xi\right)^{n}{\cal C}^{(n)}(t)\Big{|}_{n\,\rm even},$
(43a) $\displaystyle\int_{-\xi}^{1}\differential{y}y^{-1}\,f(y,\xi,t)$
$\displaystyle=$ $\displaystyle{\cal C^{\prime}}(t).$ (43b)

Equation (42) can be rewritten as

$\displaystyle\int_{\xi}^{1}\differential{x}x^{n-1}\,H_{q}(x,\xi,t)$
$\displaystyle=$ $\displaystyle\sum_{i=0,\rm
even}^{n-1}(2\xi)^{i}\sum_{j=0,\rm even}^{i}A^{n(j)}_{q/\pi}(t){\cal
A}^{n(i-j)}(t)$ (44) $\displaystyle+$
$\displaystyle(2\xi)^{n}\left(\sum_{i=0,\rm
even}^{n-1}A^{n(i)}_{q/\pi}(t)\,{\cal C}^{(n)}(t)+C^{(n)}_{q/\pi}(t)\,{\cal
C^{\prime}}(t)\right).$

Now we obtain the correct polynomiality of GPD $H_{q}$. The polynomiality of
the magnetic GPD $E_{q}$ can be obtained by replacing splitting function $f$
by $g$.

##  Appendix B Splitting function integrals

In this appendix, we present the explicit expressions for the integrals of the
splitting functions appearing in Sec. II.3, along with the meson-baryon
couplings. For Fig 1(a), the splitting functions $f_{\phi B}^{({\rm rbw})}$
and $g_{\phi B}^{({\rm rbw})}$ can be written as

$\displaystyle f_{\phi B}^{({\rm rbw})}(y,\xi,t)$ $\displaystyle=$
$\displaystyle\frac{C_{B\phi}^{2}}{f^{2}}\int\\!\frac{\differential[4]{k}}{(2\pi)^{4}}\frac{-iF_{\phi
B}^{({\rm
rbw})}}{D_{B}(p-k)D_{\phi}(k+\Delta)D_{\phi}(k)}\widetilde{F}(k)\widetilde{F}(k+\Delta)\,\delta\Big{(}y+\xi-\frac{k^{+}}{P^{+}}\Big{)},$
$\displaystyle g_{\phi B}^{({\rm rbw})}(y,\xi,t)$ $\displaystyle=$
$\displaystyle\frac{C_{B\phi}^{2}}{f^{2}}\int\\!\frac{\differential[4]{k}}{(2\pi)^{4}}\frac{-iG_{\phi
B}^{({\rm
rbw})}}{D_{B}(p-k)D_{\phi}(k+\Delta)D_{\phi}(k)}\widetilde{F}(k)\widetilde{F}(k+\Delta)\,\delta\Big{(}y+\xi-\frac{k^{+}}{P^{+}}\Big{)},$

where the factors in the numerator of the integrand are given by

$\displaystyle F_{\phi B}^{({\rm rbw})}$ $\displaystyle=$
$\displaystyle\frac{P^{+}y}{2\left(4M^{2}\xi^{2}-\xi^{2}t+t\right)}\bigg{\\{}4(k\cdot
p)^{2}\Big{[}(\xi-1)t-4M^{2}\xi^{2}\Big{]}-4k\cdot p\,k\cdot
p^{\prime}\Big{[}4M^{2}\xi^{2}+(1+\xi)\Big{]}$ (46a) $\displaystyle+$
$\displaystyle 2k\cdot
p^{\prime}\left(2M\overline{M}_{\\!B}+k^{2}\right)\left(4M^{2}\xi^{2}+\xi
t+t\right)-2k\cdot p\Big{[}2MM^{B}\left(4M^{2}\xi^{2}-2\xi^{2}t+\xi
t+t\right)$ $\displaystyle+$ $\displaystyle
k^{2}\left((\xi-1)t-4M^{2}\xi^{2}\right)+8M^{4}\xi^{2}-4M^{2}\xi^{2}t+2M^{2}\xi
t+2M^{2}t+\xi t^{2}+t^{2}y\Big{]}\bigg{\\}}$ $\displaystyle G_{\phi B}^{({\rm
rbw})}$ $\displaystyle=$
$\displaystyle\frac{2MP^{+}y}{\left(\xi^{2}-1\right)t-4M^{2}\xi^{2}}\bigg{\\{}4M(1-\xi)\xi(k\cdot
p)^{2}-4M\xi(\xi+1)k\cdot p\,k\cdot p^{\prime}$ $\displaystyle+$
$\displaystyle 2M\xi(\xi+1)k\cdot
p^{\prime}\left(2M\overline{M}_{\\!B}+k^{2}\right)+2k\cdot
p\Big{[}\left((\xi+1)M^{B}\left((\xi-1)t-2M^{2}\xi\right)\right)$
$\displaystyle+$ $\displaystyle k^{2}M(\xi-1)\xi-2M^{3}(\xi+1)\xi-
Mt\left(-\xi^{2}+\xi+y+1\right)\Big{]}$ $\displaystyle+$ $\displaystyle
k^{2}\Big{[}M^{B}\left(4M^{2}\xi^{2}-\xi^{2}t+t\right)+4M^{3}\xi^{2}+Mt\left(-\xi^{2}+\xi+y+1\right)\Big{]}+2M^{2}t\overline{M}_{\\!B}(\xi+y)\bigg{\\}}$

with $\overline{M}_{\\!B}\equiv M+M_{B}$. The coefficients $C^{2}_{B\phi}$ in
Eqs. (45) are given from the Lagrangian density in Eq. (10),

$\displaystyle C_{p\pi^{0}}$ $\displaystyle=$
$\displaystyle\frac{g_{A}^{2}}{4},$ (47a) $\displaystyle C_{n\pi^{+}}$
$\displaystyle=$ $\displaystyle\frac{g_{A}^{2}}{2}.$ (47b)

For the bubble diagram in Fig. 1(b), the splitting function $f^{({\rm
bub})}_{\phi}(y,t)$ is given by

$\displaystyle f^{({\rm
bub})}_{\phi}(y,\xi,t)=\frac{iC_{\phi\phi}}{2f^{2}}\int\\!\frac{\differential[4]{k}}{(2\pi)^{4}}\frac{F_{\phi}^{({\rm
bub})}}{D_{\phi}(k+\Delta)D_{\phi}(k)}\,\widetilde{F}(k+\Delta)\,\widetilde{F}(k)\,\delta\Big{(}y+\xi-\frac{k^{+}}{P^{+}}\Big{)},$
(48)

with

$\displaystyle F_{\phi}^{({\rm bub})}=\frac{y\,P^{+}\left[4k\cdot
P\left(4M^{2}\xi^{2}+t\right)+2\xi
t\,k\cdot\Delta+t^{2}(\xi+y)\right]}{4M^{2}\xi^{2}-\xi^{2}t+t},$ (49)

where the coefficient for the pion loop is given by

$\displaystyle C_{\pi\pi}$ $\displaystyle=$ $\displaystyle\frac{1}{2}.$ (50)

Similarly, for the additional bubble diagram in Fig. 1(c), the splitting
function $g^{\prime({\rm bub})}_{\phi}(y,\xi,t)$ can be expressed as

$\displaystyle g^{\prime({\rm
bub})}_{\phi}(y,\xi,t)=\frac{iC^{\prime}_{\phi\phi}}{f^{2}}\int\\!\frac{\differential[4]{k}}{(2\pi)^{4}}\frac{G_{\phi}^{({\rm
bub})}}{D_{\phi}(k+\Delta)D_{\phi}(k)}\,\widetilde{F}(k+\Delta)\,\widetilde{F}(k)\,\delta\Big{(}y+\xi-\frac{k^{+}}{P^{+}}\Big{)}.$
(51)

with

$\displaystyle G_{\phi}^{({\rm bub})}=\frac{2MP^{+}y\left[2k\cdot
P\left(\xi^{2}-1\right)t+4k\cdot\Delta
M^{2}\xi+2M^{2}t(\xi+y)\right]}{\left(\xi^{2}-1\right)t-4M^{2}\xi^{2}},$ (52)

and the coefficients $C^{\prime}_{\phi\phi}$ is given by

$\displaystyle C^{\prime}_{\phi\phi}$ $\displaystyle=$ $\displaystyle
2(b_{10}+b_{11}).$ (53)

For Fig 1(d), the splitting functions $f_{\phi T}^{({\rm rbw})}(y,\xi,t)$ and
$g_{\phi T}^{({\rm rbw})}(y,\xi,t)$ can be written as

$\displaystyle f^{({\rm rbw})}_{\phi T}(y,\xi,t)$ $\displaystyle=$
$\displaystyle\frac{C_{T\phi}^{2}}{f^{2}}\int\\!\frac{\differential[4]{k}}{(2\pi)^{4}}\frac{iF_{\phi
T}^{({\rm
rbw})}}{D_{T}(p-k)D_{\phi}(k+\Delta)D_{\phi}(k)}\,\widetilde{F}(k+\Delta)\,\widetilde{F}(k)\,\delta\Big{(}y+\xi-\frac{k^{+}}{P^{+}}\Big{)},$
$\displaystyle g^{({\rm rbw})}_{\phi T}(y,\xi,t)$ $\displaystyle=$
$\displaystyle\frac{C_{T\phi}^{2}}{f^{2}}\int\\!\frac{\differential[4]{k}}{(2\pi)^{4}}\frac{iG_{\phi
T}^{({\rm
rbw})}}{D_{T}(p-k)D_{\phi}(k+\Delta)D_{\phi}(k)}\,\widetilde{F}(k+\Delta)\,\widetilde{F}(k)\,\delta\Big{(}y+\xi-\frac{k^{+}}{P^{+}}\Big{)},$

where the factors in the numerator of the integrand are given by

$\displaystyle F_{\phi T}^{({\rm rbw})}$ $\displaystyle=$
$\displaystyle-\frac{P^{+}y}{12M_{T}^{2}\left(4M^{2}\xi^{2}-\xi^{2}t+t\right)}\bigg{\\{}8(k\cdot
p)^{2}k\cdot p^{\prime}\left(4M^{2}\xi^{2}+(1-\xi)t\right)$ (55a)
$\displaystyle+$ $\displaystyle 8(k\cdot p^{\prime})^{2}k\cdot
p\left(4M^{2}\xi^{2}+(1-\xi)t\right)-4k\cdot p\,k\cdot
p^{\prime}\Big{[}24M^{3}M_{T}\xi^{2}+16M^{4}\xi^{2}+4M^{2}(\xi+1)t$
$\displaystyle+$ $\displaystyle
2MM_{T}\left(2\xi^{2}+6\xi+3\right)t-t^{2}(2\xi+y+1)\Big{]}+4(k\cdot
p)^{2}\Big{[}4M^{3}M_{T}\xi^{2}+8M^{4}\xi^{2}$ $\displaystyle+$ $\displaystyle
2M^{2}\left(-2\xi^{2}+3\xi+1\right)t+MM_{T}\left(-12\xi^{2}+11\xi+1\right)t-(\xi-1)t^{2}\Big{]}$
$\displaystyle-$ $\displaystyle 4M(k\cdot
p^{\prime})^{2}(2M-M_{T})\left(4M^{2}\xi^{2}+\xi t+t\right)-2k\cdot
p\Big{[}2M^{2}t\left(t\left(2\xi^{2}+2\xi+3y-1\right)-6M_{T}^{2}\xi\right)$
$\displaystyle+$ $\displaystyle
4k^{2}M\left(4M^{3}\xi^{2}+M\left(-2\xi^{2}+\xi+1\right)t-3M_{T}(\xi-1)\xi
t\right)+32M^{5}M_{T}\xi^{2}$ $\displaystyle+$ $\displaystyle
4M^{3}M_{T}\left(-9\xi^{2}+2\xi+2\right)t+32M^{6}\xi^{2}+8M^{4}\left(-3\xi^{2}+\xi+1\right)t$
$\displaystyle+$ $\displaystyle
MM_{T}t\left(12M_{T}^{2}(\xi-1)\xi+t\left(4\xi^{2}+12\xi+11y-5\right)\right)-t^{3}(\xi+y)$
$\displaystyle+$ $\displaystyle 3t^{2}M_{T}^{2}(\xi-1)\Big{]}+2k\cdot
p^{\prime}\Big{[}32M^{5}M_{T}\xi^{2}+4M^{3}M_{T}\left(\xi^{2}+2\xi+2\right)t+32M^{6}\xi^{2}$
$\displaystyle-$ $\displaystyle
4k^{2}M\left(4M^{3}\xi^{2}-M\left(2\xi^{2}+\xi-1\right)t-3M_{T}\xi(\xi+1)t\right)-2M^{2}t\left(6M_{T}^{2}\xi+t(2\xi+y+1)\right)$
$\displaystyle+$ $\displaystyle
8M^{4}\left(-\xi^{2}+\xi+1\right)t+MM_{T}t\left(t(2\xi+y+1)-12M_{T}^{2}\xi(\xi+1)\right)+3M_{T}^{2}(\xi+1)t^{2}\Big{]}$
$\displaystyle+$ $\displaystyle
4k^{2}M\Big{[}16M^{4}M_{T}\xi^{2}+4M^{2}M_{T}\left(1-4\xi^{2}\right)t+16M^{5}\xi^{2}+4M^{3}\left(1-3\xi^{2}\right)t$
$\displaystyle+$ $\displaystyle
Mt^{2}\left(2\xi^{2}+\xi+y-2\right)+3M_{T}t^{2}\left(\xi^{2}+\xi+y-1\right)\Big{]}$
$\displaystyle+$ $\displaystyle
t^{2}\overline{M}_{T}(\xi+y)\left(8M^{3}-2M\left(6M_{T}^{2}+t\right)+3M_{T}t\right)\bigg{\\}},$
$\displaystyle G_{\phi T}^{({\rm rbw})}$ $\displaystyle=$
$\displaystyle\frac{MP^{+}y}{3M_{T}^{2}\left(4M^{2}\xi^{2}-\xi^{2}t+t\right)}\bigg{\\{}8M(\xi-1)\xi(k\cdot
p)^{2}k\cdot p^{\prime}+8M\xi(\xi+1)k\cdot p(k\cdot p^{\prime})^{2}$ (55b)
$\displaystyle-$ $\displaystyle 4k\cdot p\,k\cdot
p^{\prime}\Big{[}6M^{2}M_{T}(\xi+2)\xi+4M^{3}(\xi+1)\xi-
Mt(2\xi+y+1)+M_{T}\left(\xi^{2}-1\right)t\Big{]}$ $\displaystyle+$
$\displaystyle 4(k\cdot
p)^{2}\Big{[}M^{2}M_{T}(\xi+11)\xi+2M^{3}(\xi+3)\xi-M\left(\xi^{2}+\xi-2\right)t-3M_{T}\left(\xi^{2}-1\right)t\Big{]}$
$\displaystyle-$ $\displaystyle 4(k\cdot
p^{\prime})^{2}M^{2}\xi(\xi+1)(2M-M_{T})-2k\cdot
p\Big{[}2M^{3}\left(t\left(-3\xi^{2}+2\xi+3y+2\right)-6M_{T}^{2}\xi\right)$
$\displaystyle+$ $\displaystyle
k^{2}\left(12M^{2}M_{T}\xi+4M^{3}(\xi+1)\xi-2M\left(\xi^{2}-1\right)t-3M_{T}\left(\xi^{2}-1\right)t\right)$
$\displaystyle+$ $\displaystyle
M^{2}M_{T}\left(t\left(-9\xi^{2}+12\xi+11y+4\right)-12M_{T}^{2}\xi\right)+8M^{4}M_{T}(\xi+1)\xi+8M^{5}(\xi+1)\xi$
$\displaystyle-$ $\displaystyle
Mt\left(t\left(-\xi^{2}+\xi+y+1\right)-3M_{T}^{2}(\xi-1)\right)+M_{T}\left(\xi^{2}-1\right)t\left(3M_{T}^{2}+t\right)\Big{]}$
$\displaystyle+$ $\displaystyle 2k\cdot
p^{\prime}\Big{[}k^{2}\left(12M^{2}M_{T}\xi-4M^{3}(\xi-1)\xi+2M\left(\xi^{2}-1\right)t+3M_{T}\left(\xi^{2}-1\right)t\right)$
$\displaystyle-$ $\displaystyle
2M^{3}\left(6M_{T}^{2}\xi+t(\xi+2)\xi+yt\right)+M^{2}M_{T}\left(t(\xi+2)\xi+yt-12M_{T}^{2}\xi\right)$
$\displaystyle+$ $\displaystyle
8M^{4}M_{T}(\xi+1)\xi+8M^{5}(\xi+1)\xi+3MM_{T}^{2}(\xi+1)t-3M_{T}^{3}\left(\xi^{2}-1\right)t\Big{]}$
$\displaystyle+$ $\displaystyle
k^{2}\Big{[}16M^{4}M_{T}\xi^{2}+4M^{2}M_{T}t\left(-4\xi^{2}+3\xi+3y+1\right)+16M^{5}\xi^{2}$
$\displaystyle+$ $\displaystyle
4M^{3}t\left(-3\xi^{2}+\xi+y+1\right)+2M\left(\xi^{2}-1\right)t^{2}+3M_{T}\left(\xi^{2}-1\right)t^{2}\Big{]}$
$\displaystyle+$ $\displaystyle
Mt\overline{M}_{T}(\xi+y)\left(8M^{3}-2M\left(6M_{T}^{2}+t\right)+3M_{T}t\right)\bigg{\\}},$

with $\overline{M}_{T}\equiv M+M_{T}$ the average of the nucleon and decuplet
baryon masses. The decuplet baryon couplings for the pion case are give by

$\displaystyle C^{2}_{\Delta^{+}\pi^{0}}$ $\displaystyle=$
$\displaystyle\frac{{\cal C}^{2}}{3},$ (56a) $\displaystyle
C^{2}_{\Delta^{0}\pi^{+}}$ $\displaystyle=$ $\displaystyle\frac{{\cal
C}^{2}}{6},$ (56b) $\displaystyle C^{2}_{\Delta^{++}\pi^{-}}$ $\displaystyle=$
$\displaystyle\frac{{\cal C}^{2}}{2}$ (56c)

## References

  * Ji (1998) X.-D. Ji, J. Phys. G 24, 1181 (1998), arXiv:hep-ph/9807358 .
  * Goeke  _et al._ (2001) K. Goeke, M. V. Polyakov, and M. Vanderhaeghen, Prog. Part. Nucl. Phys. 47, 401 (2001), arXiv:hep-ph/0106012 .
  * Belitsky and Radyushkin (2005) A. V. Belitsky and A. V. Radyushkin, Phys. Rep. 418, 1 (2005), arXiv:hep-ph/0504030 .
  * Altarelli and Parisi (1977) G. Altarelli and G. Parisi, Nucl. Phys. B126, 298 (1977).
  * Gribov and Lipatov (1972) V. N. Gribov and L. N. Lipatov, Sov. J. Nucl. Phys. 15, 438 (1972).
  * Dokshitzer (1977) Y. L. Dokshitzer, Sov. Phys. JETP 46, 641 (1977).
  * Efremov and Radyushkin (1980) A. V. Efremov and A. V. Radyushkin, Phys. Lett. B 94, 245 (1980).
  * Lepage and Brodsky (1980) G. P. Lepage and S. J. Brodsky, Phys. Rev. D 22, 2157 (1980).
  * Lepage and Brodsky (1979) G. P. Lepage and S. J. Brodsky, Phys. Lett. B 87, 359 (1979).
  * Ji (1997) X. Ji, Phys. Rev. D 55, 7114 (1997), arXiv:hep-ph/9609381 .
  * Radyushkin (1996) A. V. Radyushkin, Phys. Lett. B 385, 333 (1996), arXiv:hep-ph/9605431 .
  * Mankiewicz  _et al._ (1998) L. Mankiewicz, G. Piller, and T. Weigl, Eur. Phys. J. C 5, 119 (1998), arXiv:hep-ph/9711227 .
  * Collins and Freund (1999) J. C. Collins and A. Freund, Phys. Rev. D 59, 074009 (1999), arXiv:hep-ph/9801262 .
  * Adloff  _et al._ (2001) C. Adloff _et al._ , Phys. Lett. B 517, 47 (2001), arXiv:hep-ex/0107005 .
  * Aktas  _et al._ (2005) A. Aktas _et al._ , Eur. Phys. J. C 44, 1 (2005), arXiv:hep-ex/0505061 .
  * Aaron  _et al._ (2009) F. D. Aaron _et al._ , Phys. Lett. B 681, 391 (2009), arXiv:0907.5289 [hep-ex] .
  * Breitweg  _et al._ (1999) J. Breitweg _et al._ , Eur. Phys. J. C 6, 603 (1999), arXiv:hep-ex/9808020 .
  * Chekanov  _et al._ (2003) S. Chekanov _et al._ , Phys. Lett. B 573, 46 (2003), arXiv:hep-ex/0305028 .
  * Airapetian  _et al._ (2001) A. Airapetian _et al._ , Phys. Rev. Lett. 87, 182001 (2001), arXiv:hep-ex/0106068 .
  * Airapetian  _et al._ (2011) A. Airapetian _et al._ , Phys. Lett. B 704, 15 (2011), arXiv:1106.2990 [hep-ex] .
  * Airapetian  _et al._ (2012) A. Airapetian _et al._ , J. High Energy Phys. 07, 032 (2012), arXiv:1203.6287 [hep-ex] .
  * d’Hose  _et al._ (2004) N. d’Hose, E. Burtin, P. A. M. Guichon, and J. Marroncle, Eur. Phys. J. A 19S1, 47 (2004).
  * Fuchey (2015) E. Fuchey, PoS QCDEV2015, 048 (2015).
  * Stepanyan  _et al._ (2001) S. Stepanyan _et al._ , Phys. Rev. Lett. 87, 182002 (2001), arXiv:hep-ex/0107043 .
  * Bedlinskiy  _et al._ (2012) I. Bedlinskiy _et al._ , Phys. Rev. Lett. 109, 112001 (2012), arXiv:1206.6355 [hep-ex] .
  * Seder  _et al._ (2015) E. Seder _et al._ , Phys. Rev. Lett. 114, 032001 (2015), [Addendum: Phys. Rev. Lett. 114, 089901 (2015)], arXiv:1410.6615 [hep-ex] .
  * Jo  _et al._ (2015) H. S. Jo _et al._ , Phys. Rev. Lett. 115, 212003 (2015), arXiv:1504.02009 [hep-ex] .
  * Hattawy  _et al._ (2019) M. Hattawy _et al._ , Phys. Rev. Lett. 123, 032502 (2019), arXiv:1812.07628 [nucl-ex] .
  * Burkert  _et al._ (2021) V. Burkert _et al._ , Eur. Phys. J. A 57, 186 (2021), arXiv:2103.12651 [nucl-ex] .
  * Defurne  _et al._ (2015) M. Defurne _et al._ , Phys. Rev. C 92, 055202 (2015), arXiv:1504.05453 [nucl-ex] .
  * Georges  _et al._ (2022) F. Georges _et al._ , Phys. Rev. Lett. 128, 252002 (2022), arXiv:2201.03714 [hep-ph] .
  * Belitsky  _et al._ (2002) A. V. Belitsky, D. Müller, and A. Kirchner, Nucl. Phys. B629, 323 (2002), arXiv:hep-ph/0112108 .
  * Accardi  _et al._ (2016) A. Accardi _et al._ , Eur. Phys. J. A 52, 268 (2016), arXiv:1212.1701 [nucl-ex] .
  * Avakian (2019) H. Avakian, PoS SPIN2018, 171 (2019).
  * Diehl  _et al._ (2005) M. Diehl, T. Feldmann, R. Jakob, and P. Kroll, Eur. Phys. J. C 39, 1 (2005), arXiv:hep-ph/0408173 .
  * Kumerički and Müller (2010) K. Kumerički and D. Müller, Nucl. Phys. B841, 1 (2010), arXiv:0904.0458 [hep-ph] .
  * Goldstein  _et al._ (2011) G. R. Goldstein, J. O. Hernandez, and S. Liuti, Phys. Rev. D 84, 034007 (2011), arXiv:1012.3776 [hep-ph] .
  * Guo  _et al._ (2022) Y. Guo, X. Ji, and K. Shiells, J. High Energy Phys. 09, 215 (2022), arXiv:2207.05768 [hep-ph] .
  * Guo  _et al._ (2023) Y. Guo, X. Ji, M. G. Santiago, K. Shiells, and J. Yang, (2023), arXiv:2302.07279 [hep-ph] .
  * Müller and Schäfer (2006) D. Müller and A. Schäfer, Nucl. Phys. B739, 1 (2006), arXiv:hep-ph/0509204 .
  * Alexandrou  _et al._ (2020) C. Alexandrou, K. Cichy, M. Constantinou, K. Hadjiyiannakou, K. Jansen, A. Scapellato, and F. Steffens, Phys. Rev. Lett. 125, 262001 (2020), arXiv:2008.10573 [hep-lat] .
  * Alexandrou  _et al._ (2022) C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, and A. Vaquero, PoS LATTICE2021, 250 (2022), arXiv:2112.06750 [hep-lat] .
  * Bertone  _et al._ (2021) V. Bertone, H. Dutrieux, C. Mezrag, H. Moutarde, and P. Sznajder, Phys. Rev. D 103, 114019 (2021).
  * Moffat  _et al._ (2023) E. Moffat, A. Freese, I. Cloët, T. Donohoe, L. Gamberg, W. Melnitchouk, A. Metz, A. Prokudin, and N. Sato, Phys. Rev. D 108, 036027 (2023).
  * Qiu and Yu (2023) J.-W. Qiu and Z. Yu, Phys. Rev. D 107, 014007 (2023), arXiv:2210.07995 [hep-ph] .
  * Qiu and Yu (2022) J.-W. Qiu and Z. Yu, J. High Energy Phys. 08, 103 (2022), arXiv:2205.07846 [hep-ph] .
  * Lin (2021) H.-W. Lin, Phys. Rev. Lett. 127, 182001 (2021), arXiv:2008.12474 [hep-ph] .
  * Bhattacharya  _et al._ (2022) S. Bhattacharya, K. Cichy, M. Constantinou, J. Dodson, X. Gao, A. Metz, S. Mukherjee, A. Scapellato, F. Steffens, and Y. Zhao, Phys. Rev. D 106, 114512 (2022), arXiv:2209.05373 [hep-lat] .
  * Alexandrou  _et al._ (2006) C. Alexandrou, G. Koutsou, J. W. Negele, and A. Tsapalis, Phys. Rev D. 74, 034508 (2006), arXiv:hep-lat/0605017 .
  * Ji (2013) X. Ji, Phys. Rev. Lett. 110, 262002 (2013), arXiv:1305.1539 [hep-ph] .
  * Orginos  _et al._ (2017) K. Orginos, A. Radyushkin, J. Karpie, and S. Zafeiropoulos, Phys. Rev. D 96, 094503 (2017), arXiv:1706.05373 [hep-ph] .
  * Ma and Qiu (2018a) Y.-Q. Ma and J.-W. Qiu, Phys. Rev. D 98, 074021 (2018a), arXiv:1404.6860 [hep-ph] .
  * Ma and Qiu (2018b) Y.-Q. Ma and J.-W. Qiu, Phys. Rev. Lett. 120, 022003 (2018b), arXiv:1709.03018 [hep-ph] .
  * Ji  _et al._ (1997) X. Ji, W. Melnitchouk, and X. Song, Phys. Rev. D 56, 5511 (1997), arXiv:hep-ph/9702379 .
  * Petrov  _et al._ (1998) V. Y. Petrov, P. V. Pobylitsa, M. V. Polyakov, I. Bornig, K. Goeke, and C. Weiss, Phys. Rev. D 57, 4325 (1998), arXiv:hep-ph/9710270 .
  * Penttinen  _et al._ (2000) M. Penttinen, M. V. Polyakov, and K. Goeke, Phys. Rev. D 62, 014024 (2000), arXiv:hep-ph/9909489 .
  * Choi  _et al._ (2001) H.-M. Choi, C.-R. Ji, and L. S. Kisslinger, Phys. Rev. D 64, 093006 (2001), arXiv:hep-ph/0104117 .
  * Choi  _et al._ (2002) H.-M. Choi, C.-R. Ji, and L. S. Kisslinger, Phys. Rev. D 66, 053011 (2002), arXiv:hep-ph/0204321 .
  * Tiburzi and Miller (2002) B. C. Tiburzi and G. A. Miller, Phys. Rev. D 65, 074009 (2002), arXiv:hep-ph/0109174 .
  * Theussl  _et al._ (2004) L. Theussl, S. Noguera, and V. Vento, Eur. Phys. J. A 20, 483 (2004), arXiv:nucl-th/0211036 .
  * Boffi  _et al._ (2003) S. Boffi, B. Pasquini, and M. Traini, Nucl. Phys. B649, 243 (2003), arXiv:hep-ph/0207340 .
  * Schweitzer  _et al._ (2002) P. Schweitzer, S. Boffi, and M. Radici, Phys. Rev. D 66, 114004 (2002), arXiv:hep-ph/0207230 .
  * Scopetta and Vento (2004) S. Scopetta and V. Vento, Phys. Rev. D 69, 094004 (2004), arXiv:hep-ph/0307150 .
  * Ossmann  _et al._ (2005) J. Ossmann, M. V. Polyakov, P. Schweitzer, D. Urbano, and K. Goeke, Phys. Rev. D 71, 034011 (2005), arXiv:hep-ph/0411172 .
  * Wakamatsu and Tsujimoto (2005) M. Wakamatsu and H. Tsujimoto, Phys. Rev. D 71, 074001 (2005), arXiv:hep-ph/0502030 .
  * Mineo  _et al._ (2005) H. Mineo, S. N. Yang, C.-Y. Cheung, and W. Bentz, Phys. Rev. C 72, 025202 (2005).
  * Pasquini and Boffi (2007) B. Pasquini and S. Boffi, Nucl. Phys. A782 (2007), 10.1016/j.nuclphysa.2006.10.006, arXiv:hep-ph/0607213 .
  * Goeke  _et al._ (2008) K. Goeke, V. Guzey, and M. Siddikov, Eur. Phys. J. C 56, 203 (2008), arXiv:0804.4424 [hep-ph] .
  * Thomas (1983) A. W. Thomas, Phys. Lett. B 126, 97 (1983).
  * Thomas  _et al._ (2000) A. W. Thomas, W. Melnitchouk, and F. M. Steffens, Phys. Rev. Lett. 85, 2892 (2000), arXiv:hep-ph/0005043 .
  * Chen and Ji (2001) J.-W. Chen and X. Ji, Phys. Rev. Lett. 87, 152002 (2001), [Erratum: Phys. Rev. Lett. 88, 249901 (2002)], arXiv:hep-ph/0107158 .
  * Salamu  _et al._ (2015) Y. Salamu, C.-R. Ji, W. Melnitchouk, and P. Wang, Phys. Rev. Lett. 114, 122001 (2015), arXiv:1409.5885 [hep-ph] .
  * Luan and Lu (2023) X. Luan and Z. Lu, (2023), arXiv:2302.11278 [hep-ph] .
  * Young  _et al._ (2003) R. D. Young, D. B. Leinweber, and A. W. Thomas, Prog. Part. Nucl. Phys. 50, 399 (2003), arXiv:hep-lat/0212031 .
  * Leinweber  _et al._ (2004) D. B. Leinweber, A. W. Thomas, and R. D. Young, Phys. Rev. Lett. 92, 242002 (2004), arXiv:hep-lat/0302020 .
  * Donoghue  _et al._ (1999) J. F. Donoghue, B. R. Holstein, and B. Borasoy, Phys. Rev. D 59, 036002 (1999), arXiv:hep-ph/9804281 .
  * Wang  _et al._ (2007) P. Wang, D. B. Leinweber, A. W. Thomas, and R. D. Young, Phys. Rev. D 75, 073012 (2007), arXiv:hep-ph/0701082 .
  * Allton  _et al._ (2005) C. R. Allton, W. Armour, D. B. Leinweber, A. W. Thomas, and R. D. Young, Phys. Lett. B 628, 125 (2005), arXiv:hep-lat/0504022 .
  * Wang  _et al._ (2009) P. Wang, D. B. Leinweber, A. W. Thomas, and R. D. Young, Phys. Rev. C 79, 065202 (2009), arXiv:0807.0944 [hep-ph] .
  * Wang  _et al._ (2012) P. Wang, D. B. Leinweber, A. W. Thomas, and R. D. Young, Phys. Rev. D 86, 094038 (2012), arXiv:1210.5072 [hep-ph] .
  * Hall  _et al._ (2014) J. M. M. Hall, D. B. Leinweber, and R. D. Young, Phys. Rev. D 89, 054511 (2014), arXiv:1312.5781 [hep-lat] .
  * Shanahan  _et al._ (2013) P. E. Shanahan, A. W. Thomas, and R. D. Young, Phys. Rev. D 87, 074503 (2013), arXiv:1205.5365 [nucl-th] .
  * Shanahan  _et al._ (2014) P. E. Shanahan, A. W. Thomas, R. D. Young, J. M. Zanotti, R. Horsley, Y. Nakamura, D. Pleiter, P. E. L. Rakow, G. Schierholz, and H. Stüben, Phys. Rev. D 89, 074511 (2014), arXiv:1401.5862 [hep-lat] .
  * Wang (2011) P. Wang, Chin. Phys. C 35, 223 (2011), arXiv:1006.0842 [hep-th] .
  * He and Wang (2018a) F. He and P. Wang, Phys. Rev. D 97, 036007 (2018a), arXiv:1711.05896 [nucl-th] .
  * He and Wang (2018b) F. He and P. Wang, Phys. Rev. D 98, 036007 (2018b), arXiv:1805.11986 [hep-ph] .
  * Salamu  _et al._ (2019a) Y. Salamu, C.-R. Ji, W. Melnitchouk, A. W. Thomas, and P. Wang, Phys. Rev. D 99, 014041 (2019a), arXiv:1806.07551 [hep-ph] .
  * Salamu  _et al._ (2019b) Y. Salamu, C.-R. Ji, W. Melnitchouk, A. W. Thomas, P. Wang, and X. G. Wang, Phys. Rev. D 100, 094026 (2019b), arXiv:1907.08551 [hep-ph] .
  * He and Wang (2019) F. He and P. Wang, Phys. Rev. D 100, 074032 (2019), arXiv:1904.06815 [hep-ph] .
  * He  _et al._ (2022) F. He, C.-R. Ji, W. Melnitchouk, A. W. Thomas, and P. Wang, Phys. Rev. D 106, 054006 (2022), arXiv:2202.00266 [hep-ph] .
  * Wang  _et al._ (2023) P. Wang, F. He, C.-R. Ji, and W. Melnitchouk, Prog. Part. Nucl. Phys. 129, 104017 (2023), arXiv:2206.09361 [hep-ph] .
  * Copeland and Mehen (2024) M. Copeland and T. Mehen, (2024), arXiv:2405.14965 [hep-ph] .
  * Jenkins and Manohar (1991) E. E. Jenkins and A. V. Manohar, Phys. Lett. B 255, 558 (1991).
  * Bernard (2008) V. Bernard, Prog. Part. Nucl. Phys. 60, 82 (2008), arXiv:0706.0312 [hep-ph] .
  * Hacker  _et al._ (2005) C. Hacker, N. Wies, J. Gegelia, and S. Scherer, Phys. Rev. C 72, 055203 (2005), arXiv:hep-ph/0505043 .
  * Nath  _et al._ (1971) L. M. Nath, B. Etemadi, and J. D. Kimel, Phys. Rev. D 3, 2153 (1971).
  * Kubis and Meißner (2001) B. Kubis and U.-G. Meißner, Eur. Phys. J. C 18, 747 (2001), arXiv:hep-ph/0010283 .
  * Amrath  _et al._ (2008) D. Amrath, M. Diehl, and J.-P. Lansberg, Eur. Phys. J. C 58, 179 (2008), arXiv:0807.4474 [hep-ph] .
  * Hackett  _et al._ (2023) D. C. Hackett, P. R. Oare, D. A. Pefkou, and P. E. Shanahan, Phys. Rev. D 108, 114504 (2023), arXiv:2307.11707 [hep-lat] .
  * Barry  _et al._ (2018) P. C. Barry, N. Sato, W. Melnitchouk, and C.-R. Ji, Phys. Rev. Lett. 121, 152001 (2018).
  * Novikov  _et al._ (2020) I. Novikov _et al._ , Phys. Rev. D 102, 014040 (2020), arXiv:2002.02902 [hep-ph] .
  * Barry  _et al._ (2021) P. C. Barry, C.-R. Ji, N. Sato, and W. Melnitchouk, Phys. Rev. Lett. 127, 232001 (2021).
  * Kotz  _et al._ (2023) L. Kotz, A. Courtoy, P. Nadolsky, F. Olness, and M. Ponce-Chavez, (2023), arXiv:2311.08447 [hep-ph] .
  * Aicher  _et al._ (2010) M. Aicher, A. Schafer, and W. Vogelsang, Phys. Rev. Lett. 105, 252003 (2010).
  * Teryaev (2001) O. V. Teryaev, Phys. Lett. B 510, 125 (2001), arXiv:hep-ph/0102303 .
  * Kivel and Polyakov (2002) N. Kivel and M. V. Polyakov, (2002), arXiv:hep-ph/0203264 .
  * Diehl (2003) M. Diehl, Phys. Rep. 388, 41 (2003), arXiv:hep-ph/0307382 .
  * Gao  _et al._ (2024) Z. Gao _et al._ , in preparation (2024).
  * Ackerstaff  _et al._ (1998) K. Ackerstaff _et al._ , Phys. Rev. Lett. 81, 5519 (1998), arXiv:hep-ex/9807013 .
  * Towell  _et al._ (2001) R. S. Towell _et al._ , Phys. Rev. D 64, 052002 (2001), arXiv:hep-ex/0103030 .
  * Dove  _et al._ (2021) J. Dove _et al._ , Nature 590, 561 (2021), arXiv:2103.04024 [hep-ph] .
  * Cocuzza  _et al._ (2021) C. Cocuzza, W. Melnitchouk, A. Metz, and N. Sato, Phys. Rev. D 104, 074031 (2021), arXiv:2109.00677 [hep-ph] .
  * Accardi  _et al._ (2023) A. Accardi, X. Jing, J. F. Owens, and S. Park, Phys. Rev. D 107, 113005 (2023), arXiv:2303.11509 [hep-ph] .
  * Dove  _et al._ (2023) J. Dove _et al._ , Phys. Rev. C 108, 035202 (2023), arXiv:2212.12160 [hep-ph] .