# Three-dimensional Imaging of Pion using Lattice QCD: Generalized Parton
Distributions

Heng-Tong Ding  Key Laboratory of Quark & Lepton Physics (MOE) and Institute
of Particle Physics, Central China Normal University, Wuhan 430079, China
Xiang Gao  Physics Division, Argonne National Laboratory, Lemont, IL 60439,
USA  Swagato Mukherjee  Physics Department, Brookhaven National Laboratory,
Upton, New York 11973, USA  Peter Petreczky  Physics Department, Brookhaven
National Laboratory, Upton, New York 11973, USA  Qi Shi<EMAIL_ADDRESS>Key
Laboratory of Quark & Lepton Physics (MOE) and Institute of Particle Physics,
Central China Normal University, Wuhan 430079, China  Physics Department,
Brookhaven National Laboratory, Upton, New York 11973, USA  Sergey Syritsyn
Department of Physics and Astronomy, Stony Brook University, Stony Brook, New
York 11790  Yong Zhao  Physics Division, Argonne National Laboratory, Lemont,
IL 60439, USA

###### Abstract

In this work, we report a lattice calculation of $x$-dependent valence pion
generalized parton distributions (GPDs) at zero skewness with multiple values
of the momentum transfer $-t$. The calculations are based on an $N_{f}=2+1$
gauge ensemble of highly improved staggered quarks with Wilson-Clover valence
fermion. The lattice spacing is 0.04 fm, and the pion valence mass is tuned to
be 300 MeV. We determine the Lorentz-invariant amplitudes of the quasi-GPD
matrix elements for both symmetric and asymmetric momenta transfers with
similar values and show the equivalence of both frames. Then, focusing on the
asymmetric frame, we utilize a hybrid scheme to renormalize the quasi-GPD
matrix elements obtained from the lattice calculations. After the Fourier
transforms, the quasi-GPDs are then matched to the light-cone GPDs within the
framework of large momentum effective theory with improved matching, including
the next-to-next-to-leading order perturbative corrections, and leading
renormalon and renormalization group resummations. We also present the
3-dimensional image of the pion in impact-parameter space through the Fourier
transform of the momentum transfer $-t$.

##  I Introduction

Since its discovery in 1947, the pion has been the subject of intense
research, recognized for its dual identity as both a Goldstone boson linked to
chiral symmetry breaking and a QCD-bound state. Over the years, substantial
efforts have been dedicated to studying its internal structure by analyzing
experimental data. Key methodologies have involved extracting form factors
(FFs) through the pion-electron scattering [1] and discerning parton
distribution functions (PDFs) via the Drell-Yan process [2]. However, these
approaches are limited to revealing only the one-dimensional structure of the
hadron. For a more comprehensive, three-dimensional perspective, the focus has
shifted to generalized parton distributions (GPDs), a concept introduced in
the 1990s [3, 4, 5, 6]. Access to GPDs is typically gained through exclusive
reactions, notably deeply virtual Compton scattering [6] and deeply virtual
meson production [7, 8].

Presently, a series of experiments are being conducted or planned, which can
enrich our understanding of the pion internal structure. Prominent among these
are the JLab 12 GeV program [9], the Apparatus for Meson and Baryon
Experimental Research (AMBER) at CERN SPS [10], the forthcoming Electron-Ion
Collider (EIC) at Brookhaven National Laboratory [11], and the Electron-Ion
Collider in China (EicC) [12]. These initiatives are designed to probe the
pion structure in various kinematic regimes. However, extracting GPDs from
experimental data is fraught with challenges, including the chiral-odd nature
of certain distributions, such as the transversity GPDs, and the complexity
involved in pion production. On the other hand, lattice QCD results, offering
complementary insights and potentially guiding experimental efforts, are
highly sought after. The advent of Large Momentum Effective Theory (LaMET)
[13, 14, 15] in 2013 made it possible to directly calculate the Bjorken-$x$
dependence of the parton distributions, extending the scope of GPD studies
beyond just the first few Mellin moments [16, 17, 18, 19, 20, 21, 22, 23, 24,
25]. This spurred a flurry of lattice research on nucleon and meson GPDs [26,
27, 28, 29, 30, 31, 32, 33, 34] based on the quasi-distribution (see reviews
in Refs. [15, 35]), which is defined in terms of matrix elements of non-local
operators involving quark and antiquark fields separated by a spatial
distance, rather than a light-cone distance, and thus is calculable on the
lattice.

While there are several lattice QCD studies of the nucleon GPDs [27, 28, 29,
30, 31, 33, 34], including a few based on alternate methods [36, 37, 38, 39],
there have been limited LaMET calculations of the pion GPDs [26, 32].
Traditionally, GPD calculations are conducted in the Breit (symmetric) frame,
requiring the momentum transfer to be symmetrically distributed between the
initial and final states. However, this approach incurs substantial
computational costs. Recently, a frame-independent method was proposed to
extract GPDs from lattice calculations in an arbitrary frame [30]. This
innovative approach holds the potential to significantly reduce computational
expenses by working in a non-Breit (asymmetric) frame where the initial or
final state momentum is fixed.

Two crucial aspects of obtaining the light-cone parton distributions from the
lattice calculations are the renormalization and matching process. The GPDs
and quasi-GPDs (qGPDs) are typically defined in different renormalization
schemes. For the qGPDs extracted from the lattice analyses, renormalization is
required to eliminate the ultraviolet (UV) divergences stemming from the
Wilson line. Commonly employed renormalization methods include the
regularization-independent momentum-subtraction (RI/MOM) [40, 41, 42, 43] and
various ratio schemes [44, 45, 46, 47, 48]. However, these methods adhere to a
factorization relation with the light-cone correlation only at short
distances. At long distances, they introduce nonperturbative effects [49],
which will impact the qGPDs through the Fourier transform of the matrix
elements, thereby affecting the LaMET matching results in Bjorken-$x$ space.
To overcome this issue, in this work, we use a hybrid-scheme renormalization
[49]. The key point of this scheme is to renormalize the matrix elements at
short and long distances separately, which can remove the linear divergence at
long distances without introducing additional nonperturbative effects.

Furthermore, we match the qGPDs in the lattice renormalization scheme to the
light-cone GPDs in the $\overline{\rm MS}$ scheme through LaMET [50, 51, 52,
53]. During this process, the accuracy of the perturbative calculation plays a
significant role in the precision of the final GPD results. In the zero-
skewness limit, the matching kernel of GPDs is the same as the one for PDFs
[51], which has been derived up to the next-to-next-to-leading order (NNLO)
[47, 54]. We also utilize renormalization group resummation (RGR) [55, 56] to
resum the small-$x$ logarithmic terms. Additionally, we consider leading-
renormalon resummation (LRR) [57, 58], which can remove the renormalon
ambiguity in the Wilson-line mass matching, to eliminate the linear power
corrections [58]. Therefore, we achieve a state-of-the-art calculation of the
valence pion GPDs using an adapted hybrid-scheme renormalization, along with
the implementation of the combined NNLO+RGR+LRR matching.

In this study, we present our lattice calculations of the valence pion GPDs
using the LaMET approach, featuring a very small lattice spacing of 0.04 fm
within the non-Breit frame. The organization of this paper is as follows: In
Section II, we review the definitions of pion GPDs and qGPDs in an arbitrary
frame, outline our lattice setup, and describe the methodologies employed to
extract the matrix elements from the combined analyses of the two-point and
three-point correlation functions. In this section, we also detail the
renormalization procedure of the qGPD matrix elements in the coordinate space
using the hybrid scheme and discuss the challenges of performing the Fourier
transform to get the qGPDs. Subsequently, we apply the LaMET matching approach
to derive the valence light-cone GPDs in Section III. This section contains
our main results, including the sensitivity of our results to the perturbative
accuracy of the matching and the dependence on the renormalization scales.
Finally, Section IV provides a summary of our findings.

##  II  Valence pion quasi-GPDs in asymmetric frame on the lattice

###  II.1 General considerations

Our goal is to calculate the valence pion GPDs in an asymmetric frame using
the frame-independent approach laid out in Ref. [30]. For this, we consider
the following non-local iso-vector matrix element

$M^{\mu}(z,\bar{P},\Delta)=\left\langle\pi(P^{f})|O^{\gamma_{\mu}}(z)|\pi(P^{i})\right\rangle,$
(1)

where $\mu=(t,x,y,z)$ represents the Lorentz indices, $P^{i}$ and $P^{f}$ are
the initial and final momenta of the pion, respectively,
$\bar{P}=(P^{i}+P^{f})/2$, $\Delta=P^{f}-P^{i}$, and

$O^{\gamma_{\mu}}(z)=\frac{1}{2}\left[\bar{u}(-\frac{z}{2})\gamma^{\mu}{\cal
W}_{-\frac{z}{2},\frac{z}{2}}u(\frac{z}{2})-\bar{d}(-\frac{z}{2})\gamma^{\mu}{\cal
W}_{-\frac{z}{2},\frac{z}{2}}d(\frac{z}{2})\right],$ (2)

where ${\cal W}_{-\frac{z}{2},\frac{z}{2}}$ is the Wilson line that connects
the quark and antiquark fields, ensuring the gauge invariant of the matrix
element. Here, we use the following normalization of the pion states
$\left\langle\pi(P)|\pi(P^{\prime})\right\rangle=E_{p}(2\pi)^{3}\delta^{3}(\mathbf{P}-\mathbf{P^{\prime}})$.
Following Ref. [30], we write the matrix element in terms of the Lorentz-
invariant amplitudes $A_{i}(z\cdot\bar{P},z\cdot\Delta,\Delta^{2},z^{2})$ for
$i=1,2,3$, as follows:

$M^{\mu}(z,\bar{P},\Delta)=\bar{P}^{\mu}A_{1}+m_{\pi}^{2}z^{\mu}A_{2}+\Delta^{\mu}A_{3}.$
(3)

The coordinate space light-cone valence pion GPDs $H$ is usually defined in
terms of $M^{+}(z,\bar{P},\Delta)$ as
$M^{+}(z,\bar{P},\Delta)=P^{+}H(z,\bar{P},\Delta)$, i.e. the momentum space
light-cone GPDs is written as

$H(x,\xi,t)=P^{+}\int\frac{dz^{-}}{2\pi}e^{iz^{-}P^{+}x}H(z,\bar{P},\Delta),$
(4)

where the skewness parameter $\xi=-\Delta^{+}/(2P^{+})$ and $t=\Delta^{2}$.
However, it is possible to define the light-cone GPDs in a frame-independent,
i.e., Lorentz-invariant way [30] as

$H_{\rm
LI}(z\cdot\bar{P},z\cdot\Delta,\Delta^{2},0)=A_{1}(z\cdot\bar{P},z\cdot\Delta,\Delta^{2},0)+\frac{z\cdot\Delta}{z\cdot\bar{P}}A_{3}(z\cdot\bar{P},z\cdot\Delta,\Delta^{2},0).$
(5)

Motivated by this, the Lorentz-invariant definition of qGPDs can be written as

$\widetilde{H}_{\rm LI}(z\cdot\bar{P},z\cdot\Delta,\Delta^{2},z^{2})\equiv
A_{1}(z\cdot\bar{P},z\cdot\Delta,\Delta^{2},z^{2})+\frac{z\cdot\Delta}{z\cdot\bar{P}}A_{3}(z\cdot\bar{P},z\cdot\Delta,\Delta^{2},z^{2}).$
(6)

This is a natural choice for the qGPDs because in the light-cone limit
$z^{2}\rightarrow 0$, it should be equal to the light-cone GPDs up to the
leading order in $\alpha_{s}$

$\lim_{z^{2}\rightarrow 0}\widetilde{H}_{\rm
LI}(z\cdot\bar{P},z\cdot\Delta,\Delta^{2},z^{2})=H_{\rm
LI}(z\cdot\bar{P},z\cdot\Delta,\Delta^{2},0)+{\cal O}(\alpha_{s}).$ (7)

In the Lorentz-invariant framework, the momentum-space GPD is defined as [30]

$H(x,\xi,t)=\int\frac{d(z\cdot\bar{P})}{2\pi}e^{ixz\cdot\bar{P}}H_{\rm
LI}(z\cdot\bar{P},-2\xi(z\cdot\bar{P}),t,0),$ (8)

and similarly for qGPD

$\widetilde{H}(x,\xi,t)=\int\frac{d(z\cdot\bar{P})}{2\pi}e^{ixz\cdot\bar{P}}\widetilde{H}_{\rm
LI}(z\cdot\bar{P},-2\xi(z\cdot\bar{P}),t,z^{2}).$ (9)

By studying the behavior under hermiticity and time-reversal transformations
simultaneously, it was found that the amplitude $A_{3}$ is an odd function of
$z\cdot\Delta=-2\xi(z\cdot\bar{P})$ [30, 33]. In the forward limit, GPDs
should smoothly approach PDFs, implying that $A_{3}$ should be equal to 0 in
this work at zero skewness. In the following subsections, we will discuss our
lattice setup and lattice QCD results on $A_{i}$.

###  II.2 Lattice setup

Our lattice calculations use the gauge ensembles provided by the HotQCD
collaboration [59], utilizing a 2+1 flavor setup with Highly Improved
Staggered Quark (HISQ) action [60]. The lattice configuration has dimensions
of $N_{s}\times N_{t}=64^{3}\times 64$ and a lattice spacing of $a=0.04$ fm.
The sea quark masses are adjusted to yield a pion mass of 160 MeV. In the
valence sector, we use the Wilson-Clover action with one level of hypercubic
(HYP) smearing [61]. For the coefficient of the clover term, we utilize the
tree-level tadpole-improved value computed with the fourth root of the
plaquette, yielding 1.02868 [62]. The valence quark masses in the Wilson-
Clover action are tuned to $am_{q}=-0.033$, resulting in a pion mass of 300
MeV.

Central to our computational approach is the use of momentum-boosted smeared
sources [63]. The quark propagators are obtained through the application of
the multigrid algorithm [64] to invert the Wilson-Dirac operator using the
QUDA software suite [65, 66, 67] on GPUs. In this work, we consider the pion
boosted along the $z$-direction with momenta $P_{z}=2\pi n_{z}/(aN_{s})$ with
$n_{z}$ being an integer. To obtain a good signal for the boosted pions,
momentum-boosted sources and sinks are constructed using a Gaussian profile,
with boost momenta $k_{z}=2\pi m_{z}/(aN_{s})$ in the $z$-direction. Source
construction is done in the Coulomb gauge, and the Gaussian profile is created
with the radius of 0.208 fm [68, 62]. Employing these quark propagators, we
have computed both the two-point and three-point hadron correlation functions.
To increase statistics per configuration, we combined multiple exact (high-
precision) and sloppy (low-precision) samples and implemented the all-mode
averaging (AMA) technique [69]. The lattice parameters used in this study are
summarized in Table 1.

We can get the matrix elements by analyzing the two-point and three-point
correlation functions. The two-point correlation function is defined as

$C_{{\rm
2pt}}(\mathbf{P},t_{s})=\sum_{\mathbf{x}}e^{-i\mathbf{P}\cdot\mathbf{x}}\left\langle\pi_{\rm
s}(\mathbf{x},t_{s})\pi_{\rm s}^{\dagger}(\mathbf{0},0)\right\rangle,$ (10)

where $\mathbf{P}$ denotes the spatial momentum, $\mathbf{x}$ and $\mathbf{0}$
represent the spatial coordinates, and $t_{s}$ and $0$ correspond to the time
coordinates. Here, $\pi^{\dagger}_{s}$ and $\pi_{s}$ stand for the pion
creation and annihilation operators, respectively, with the subscript “${\rm
s}$” indicating “smeared”. The three-point correlation function is defined as

$C_{{\rm
3pt}}^{\mu}(\mathbf{P}^{f},\mathbf{q};t_{s},\tau;z)=\sum_{\mathbf{x,y}}e^{-i\mathbf{P}^{f}\cdot\mathbf{x}}e^{i\mathbf{q}\cdot\mathbf{y}}\left\langle\pi_{\rm
s}(\mathbf{x},t_{s})O^{\gamma_{\mu}}(\mathbf{y},\tau;z)\pi_{\rm
s}^{\dagger}(\mathbf{0},0)\right\rangle,$ (11)

where $\mathbf{P}^{f}=2\pi\mathbf{n}^{f}/(aN_{s})$ and $\mathbf{q}$ denote the
spatial final momentum and momentum transfer, respectively. The spatial
initial momentum is given by $\mathbf{P}^{i}=\mathbf{P}^{f}-\mathbf{q}$. The
quark bilinear operator, $O^{\gamma_{\mu}}(z)$, defined in Eq. (2), is also
characterized by its spacetime insertion position ($\mathbf{y},\tau$). In this
work, we only consider the case of zero skewness, that is, $n_{z}^{\Delta}=0$
and $P_{z}=P_{z}^{f}=P_{z}^{i}$. For each value of $P_{z}$, we consider
several values of $\mathbf{q}=2\pi\mathbf{n}^{\Delta}/(aN_{s})$. The different
choices of momenta used in our study are summarized in Table 1, including the
momentum transfer $t$ in physical units. As one can see from the table, for
the smallest non-zero value of $P_{z}$, we have performed calculations in both
Breit and non-Breit frames at very similar values of $t$.

The approach used to obtain the pion matrix element here closely follows our
previous works [62, 70, 71, 72, 73, 74]. We fit the large $t_{s}$ behavior of
the two-point correlation functions with two or three states and obtain the
energies of these states for different momenta. The corresponding energy
levels are then used in the analysis of the three-point function and the
extraction of the pion matrix elements $M^{\mu}$. We consider the following
ratio of

$R^{\mu}(\mathbf{P}^{f},\mathbf{P}^{i};t_{s},\tau;z)\equiv\sqrt{E_{0}^{f}E_{0}^{i}}\frac{C_{{\rm
3pt}}^{\mu}(\mathbf{P}^{f},\mathbf{P^{i}};t_{s},\tau;z)}{C_{{\rm
2pt}}(\mathbf{P}^{f},t_{s})}\times\left[\frac{C_{{\rm
2pt}}(\mathbf{P}^{i},t_{s}-\tau)C_{{\rm 2pt}}(\mathbf{P}^{f},\tau)C_{{\rm
2pt}}(\mathbf{P}^{f},t_{s})}{C_{{\rm 2pt}}(\mathbf{P}^{f},t_{s}-\tau)C_{{\rm
2pt}}(\mathbf{P}^{i},\tau)C_{{\rm 2pt}}(\mathbf{P}^{i},t_{s})}\right]^{1/2},$
(12)

where $E_{0}^{i}$ and $E_{0}^{f}$ correspond to the ground-state pion energies
in the initial and final states, respectively, and the kinematic factor
$\sqrt{E_{0}^{f}E_{0}^{i}}$ is used for the normalization. For large $t_{s}$
and $\tau$, this ratio gives the matrix element $M^{\mu}$. We assume that for
large $t_{s}$ and $\tau$, two or three states mainly contribute to
$R^{\mu}(\mathbf{P}^{f},\mathbf{P}^{i};t_{s},\tau;z)$ and perform the
corresponding fits to obtain the matrix element. Some details of the fits are
provided in Appendix A.

Table 1:  Details of the measurements and statistics for the Breit and non-Breit frames are shown. The symbol $t_{s}$ represents the source-sink separation. We present the momentum in units of $2\pi/(aN_{s})$, including the final momentum ($\mathbf{n}^{f}$) and the boost momentum along the $z$-direction ($m_{z}$). $P_{z}$ denotes the physical values of the momentum in the $z$-direction. For the momentum transfer, we show both the three-dimensional lattice unit momentum transfer ($\mathbf{n}^{\Delta}=\mathbf{n}^{f}-\mathbf{n}^{i}$), where $\mathbf{n}^{i}$ denotes the initial momentum, and the 4-dimensional physical unit momentum transfer ($t\equiv\Delta^{2}$). Additionally, we provide the numbers of gauge configurations (#cfgs) as well as the counts for both exact (#ex) and sloppy (#sl) inversion samples per configuration. Frame | $t_{s}/a$ | $\mathbf{n}^{f}=(n^{f}_{x},n^{f}_{y},n^{f}_{z})$ | $m_{z}$ |  $P_{z}$[GeV] | $\mathbf{n}^{\Delta}=(n^{\Delta}_{x},n^{\Delta}_{y},n^{\Delta}_{z})$ |  $-t$[GeV2] | #cfgs | (#ex, #sl)
---|---|---|---|---|---|---|---|---
Breit | 9,12,15,18 | (1, 0, 2) | 2 | 0.968 | (2, 0, 0) | 0.938 | 115 | (1, 32)
non-Breit | 9,12,15,18 | (0,0,0) | 0 | 0 | (0,0,0) | 0 | 314 | (3, 96)
9,12,15,18 | (0,0,2) | 2 | 0.968 | (1,2,0) | 0.952 | 314 | (4, 128)
9,12,15 | (0,0,3) | 2 | 1.453 |  [(0,0,0), (1,0,0) (1,1,0), (2,0,0) (2,1,0), (2,2,0)] |  [0, 0.229, 0.446, 0.855, 1.048, 1.589] | 314 | (4, 128)
9,12,15 | (0,0,4) | 3 | 1.937 |  [0, 0.231, 0.455, 0.887, 1.095, 1.690] | 564 | (4, 128)

###  II.3 Lattice results on the Lorentz-invariant amplitudes

Having determined the matrix elements $M^{\mu}$ for different values of the
kinematic variables in Eq. (3), we can deduce the Lorentz-invariant amplitudes
$A_{i}$. The results of $A_{i}$ obtained from the Breit and non-Breit frames
with comparable values of $-t$ should be consistent. Therefore, we select two
sets of lattice data from different frames, as shown in the first and third
rows of Table 1, which have identical values of $P_{z}$ and very close $-t$
values. The dataset from the Breit frame has $-t=0.938$ GeV2, while the
dataset from the non-Breit frame has $-t=0.952$ GeV2. According to Eq. (3) and
$z^{\mu}=(0,0,0,z)$ in our case, for the non-Breit frame, we can
simultaneously solve $A_{1}$ and $A_{3}$ by combining the data of $\gamma_{t}$
and $\gamma_{\perp}$ components. For the Breit frame with $\Delta^{t}=0$ and
$\bar{P}_{\perp}=0$, we can directly get $A_{1}$ and $A_{3}$ from $\gamma_{t}$
and $\gamma_{\perp}$ components, respectively.

The amplitudes obtained from these two datasets with $P_{z}=0.968$ GeV are
shown in Fig. 1. The left two panels show the results obtained using both
2-state and 3-state fits within the Breit and non-Breit frames, respectively.
The 2-state fit results are represented by the open symbols, while the 3-state
fit results are denoted as the filled symbols, labeled as “2s” and “3s” in the
figure, respectively. One can see that the contamination from the excited
states is apparent in the Breit frame, so the 3-state fit is needed. For the
non-Breit frame, it can be seen that these two kinds of fit results are
comparable within 1-$\sigma$ error, and the 2-state fit results are more
stable. In both frames, we find that $A_{3}$ is zero within error, as
expected. In the right panel of Fig. 1, we investigate the frame independence
of $A_{1}$ by selecting the 3-state fit results for the Breit frame and the
2-state fit results for the non-Breit frame. Furthermore, since we observe
that $A_{3}$ is compatible with zero within errors also for the non-Breit
frame, we perform the extraction of $A_{1}$ from the non-Breit frame based on
setting $A_{3}=0$ and using the $\gamma_{t}$ component. The results from this
analysis are shown as the green triangles in the right panel of Fig. 1. As
depicted in the figure, the non-Breit frame results, both with and without
setting $A_{3}=0$, exhibit good agreement and also are consistent with the
Breit frame results.

Figure 1:  The results for the amplitudes $A_{i}(i=1,3)$, obtained from the
Breit frame (left panel) and non-Breit frame (middle panel) with identical
$P_{z}=0.968$ GeV and very similar transfer momenta, are shown, including the
outcomes of both the 2-state and 3-state fits. In the right panel, we compare
the results of $A_{1}$ obtained from the Breit frame using the 3-state fit
with two different determinations of $A_{1}$ obtained from the non-Breit frame
using the 2-state fits with and without setting $A_{3}$ to zero, see text.

Therefore, we conclude that it is possible to obtain $A_{1}$ from the
calculations in the non-Breit frame by initially setting $A_{3}$ to zero.
Specifically, since in our lattice QCD calculations, both the initial and
final pion states are boosted along the $z$ direction by an equal amount
$P_{z}$, the Lorentz-invariant definition for the coordinate-space qGPDs for
$\xi=0$ can be directly expressed as

$\widetilde{H}_{\rm
LI}(zP_{z},\Delta^{2},z^{2})=A_{1}(zP_{z},\Delta^{2},z^{2})=M^{t}(zP_{z},\Delta^{2},z^{2})/\bar{P}^{t}.$
(13)

For convenience, we denote the dimensionless matrix element as $M\equiv
M/\bar{P}$ hereafter.

Since we aim to get the valence pion GPDs within the LaMET framework, we
selected the two largest available momenta in the non-Breit frame for our
calculations, as detailed in the last two rows of Table 1. Additionally, these
two datasets can be used to study the $P_{z}$ dependence of the final light-
cone GPDs, which can reflect the effectiveness of our perturbative matching,
as discussed later. For each momentum, we analyzed six values of the momentum
transfer $-t$. Using the matrix elements we have obtained, several steps are
required to derive the qGPDs, including renormalization, extrapolation for the
large region of $z$, and the Fourier transform. Subsequently, the light-cone
GPDs can be obtained from the qGPDs using the LaMET.

###  II.4 Renormalization

As we mentioned before, the renormalization of the qGPD matrix elements is
crucial for removing the UV divergences from the Wilson line as well as for
matching to the light-cone GPDs, which are usually defined in the
$\overline{\rm MS}$ scheme. Since the non-local quark bilinear operator is
multiplicatively renormalizable, if we call the matrix elements extracted
directly from the lattice calculations as the bare matrix elements, the
relation between the bare and renormalized matrix elements can be expressed as
[75, 76, 77, 49]

$M^{\rm B}(z,a)=Z(a)e^{-\delta m(a)|z|}e^{-\bar{m}_{0}|z|}M^{\rm R}(z,a),$
(14)

where the superscripts “B” and “R” denote the “bare” and “renormalized”
quantities, respectively, $Z(a)$ includes the $z$-independent logarithmic
divergence, the term with $\delta m(a)$ accounts for the linear divergence,
and the term with $\bar{m}_{0}$ is introduced to address the scheme dependence
of $\delta m$ and to match the lattice scheme to the $\overline{\rm MS}$
scheme [78, 70, 57, 58]. Theoretically, $\bar{m}_{0}$ should be a constant and
independent of $z$.

The hybrid scheme renormalization is defined by merging the ratio scheme for
short distances with the explicit subtraction of self-energy divergences in
the Wilson line for long distances [49]

$M^{\rm R}(z,z_{s};P_{z},t)=\left\\{\begin{aligned} &\frac{M^{\rm
B}(z,P_{z},t)}{M^{\rm B}(z,0,0)},&|z|\leq|z_{s}|;\\\ &\frac{M^{\rm
B}(z,P_{z},t)}{M^{\rm B}(z_{s},0,0)}e^{(\delta
m+\bar{m}_{0})|z-z_{s}|},&|z|>|z_{s}|.\end{aligned}\right.$ (15)

Here, $z_{s}$ represents the position where the ratio scheme is matched onto
the explicit subtraction of the divergence in the Wilson line, which is part
of the renormalization scheme. To reduce artifacts and improve the signal, we
can further divide $M^{\rm R}$ by the renormalized electromagnetic form factor
$F(P_{z},t)\equiv M^{\rm B}(0,P_{z},t)/M^{\rm B}(0,0,0)$ such as

$M^{\bar{\rm R}}(z,z_{s};P_{z},t)=M^{\rm R}(z,z_{s};P_{z},t)/F(P_{z},t).$ (16)

Notably, we need to multiply $F(P_{z},t)$ back into our results when
calculating the pion qGPDs, as presented in the next subsection.

At short distances, all the divergences can be directly canceled by such a
ratio. However, for long distances, in order to perform the renormalization,
we first need to determine the values of $\delta m$ and $\bar{m}_{0}$. In line
with our previous studies, we determine $\delta m$ using lattice QCD results
on the static quark-antiquark potential and the free energy of a static quark
at non-zero temperatures [79]: $a\delta m=0.1508(12)$ for $a=0.04$ fm lattice
[70]. The value of $\bar{m}_{0}$ can be obtained using the bare matrix element
results at zero momentum and zero momentum transfer. Namely, by comparing the
lattice computations at $P_{z}=t=0$ with their corresponding $\overline{\rm
MS}$ Operator Product Expansion (OPE) expressions for such a ratio

$e^{(\delta m+\bar{m}_{0})\delta z}\frac{M^{\rm B}(z+\delta z)}{M^{\rm
B}(z)}=\frac{C_{0}(\alpha_{s}(\mu_{0}(z+\delta z)),\mu_{0}^{2}(z+\delta
z)^{2})}{C_{0}(\alpha_{s}(\mu_{0}(z)),\mu_{0}^{2}z^{2})}\times\exp\left[\int_{\alpha_{s}(\mu_{0}(z+\delta
z))}^{\alpha_{s}(\mu_{0}(z))}{d\alpha_{s}(\mu^{\prime})\over\beta[\alpha_{s}(\mu^{\prime})]}\gamma_{O}[\alpha_{s}(\mu^{\prime})]\right],$
(17)

where $\gamma_{O}$ is the anomalous dimension of the quark bilinear operator,
known up to three loops [80], and the expressions for $\alpha_{s}$ and $\beta$
can be found in Appendix B, we can extract the values of $\bar{m}_{0}$. Here,
we set $|\delta z|/a=1$. The use of the zero momentum matrix element results
in only the 0-th Wilson coefficient $C_{0}$ remaining, which is calculated up
to the NNLO [47, 54]. Furthermore, we have incorporated RGR [55] at the next-
to-next-to-leading-logarithm (NNLL) accuracy. The RGR allows us to effectively
sum the large logarithmic terms arising in the perturbative expansions and to
evolve the running coupling from the physical scale $\mu_{0}=2\kappa
e^{-\gamma_{E}}/|z|$ to the factorization scale $\mu$ = 2 GeV. Here, $\kappa$
is a proportionality constant, and $\gamma_{E}$ is the usual Euler constant.
In principle, we should choose $\kappa\sim 1$ to match the natural physical
scale $1/|z|$, but we found that Eq. (17) with $\kappa<1$ cannot describe the
lattice data with a constant $\bar{m}_{0}$ in $z$. Therefore, to estimate the
scale variation uncertainty, we choose our central value of $\kappa$ as 1.414
and vary it between 1 and 2, where the strong coupling constant remains
reasonably small ($\alpha_{s}<$ 0.3) at the short distances we considered.
Additionally, we include LRR [58, 57] in the Wilson coefficient that captures
the dominant contributions from renormalons and improves the convergence
properties of the perturbative series [58].

In Fig. 2, we show the results of $\bar{m}_{0}$ obtained with NNLO(NLO)+LRR
and NNLO(NLO)+LRR+RGR Wilson coefficients. As mentioned before, $\bar{m}_{0}$
is excepted to be a constant and independent of $z$. However, in practice,
this is not always the case for all $z$ values, as illustrated in Fig. 2. For
the two smallest values of $z$, we observe lattice artifacts in the
determination of $\bar{m}_{0}$, as expected. However, the dependence of
$\bar{m}_{0}$ on $z$ becomes much milder for larger $z$ values. The only
exception is the case $\kappa=1$, where the running coupling becomes very
large for $z>$ 0.2 fm, causing the perturbative expansions for the Wilson
coefficients to break down. The values of $\bar{m}_{0}$ determined using
NLO+LRR also show some dependence on $z$ due to missing higher-order
corrections. However, thanks to LRR, the dependence of $\bar{m}_{0}$ on $z$ is
significantly reduced already at this level compared to our prior analyses of
the pion PDFs [70]. All other determinations of $\bar{m}_{0}$ obtained with
$\kappa\geq 1.414$ are consistent with each other and show very small
$z$-dependence for $z>0.08$ fm. In the following, we will take the
$\bar{m}_{0}$ results at $z=0.12$ fm to perform the hybrid-scheme
renormalization. For example, $\bar{m}_{0}=0.1922(3)$ is taken in the case of
NNLO+LRR+RGR ($\kappa=1$). We also choose $z_{s}=0.12$ fm as our default
choice for the hybrid scheme.

Figure 2: The additional Wilson line renormalization mass parameter
$\bar{m}_{0}$, obtained using NNLO(NLO)+LRR and NNLO(NLO)+LRR+RGR expressions
for $C_{0}$ at different values of $z$. The filled symbols correspond to NNLO
results for $C_{0}$, while the open symbols correspond to NLO results for
$C_{0}$.

###  II.5 Large $z$ extrapolation and Fourier transform

In this work, our aim is to get the valence light-cone GPDs in momentum ($x$)
space. In the isospin symmetric limit, the valence light-cone GPDs, as well as
qGPDs, are equal to the isovector GPDs, as discussed in Ref. [68]. To obtain
the momentum-space isovector (valence) qGPDs, we need to perform a Fourier
transform on $\lambda=zP_{z}$ of the coordinate-space isovector qGPDs

$\widetilde{H}(x,P_{z},t)=2\int_{-\infty}^{\infty}\frac{{\rm
d}\lambda}{2\pi}e^{ix\lambda}\widetilde{H}_{\rm
LI}(zP_{z},t,z^{2})=F(P_{z},t)\int_{-\infty}^{\infty}\frac{{\rm
d}\lambda}{\pi}e^{ix\lambda}M^{\bar{\rm R}}(z,P_{z},t),$ (18)

where $\widetilde{H}(x,P_{z},t)$ is the valence pion qGPDs, and the factor of
2 at the beginning of the middle formula comes from the definition of
$O^{\gamma_{\mu}}$ in Eq. (2) and the fact that the integral of the valence
$u$ quark distribution in the momentum space from $0$ to $1$ should be one.

In Fig. 3, we present our results of the renormalized matrix elements, denoted
by the red square symbols, as a function of $\lambda$. Due to the large
statistical errors at large $z$ and the finite volume effects, we can reliably
calculate the renormalized matrix elements only up to a certain value of
$\lambda$. This prevents us from performing the Fourier transform directly on
$M^{\bar{\rm R}}$. To address this problem, we perform an extrapolation of the
renormalized matrix elements from a truncation position $z_{t}$ and replace
the lattice results with the extrapolated values for $z\geq z_{t}$. We expect
that $M^{\bar{\rm R}}$ should vanish exponentially at large $z$ [49, 70],
which means that the extrapolation mainly affect the small-$x$ region of GPDs
[70]. Therefore, we fit the lattice results of $M^{\bar{\rm R}}$ at large $z$
using the following ansatz

$M^{\bar{\rm R}}=A\frac{e^{-mz}}{(zP_{z})^{d}},$ (19)

where $\\{A,m,d\\}$ are the fit parameters. In practice, we choose 6 data
points above a critical position $z_{c}$ where the matrix element results
become unreliable, such as showing negativity or lacking a decaying trend, to
do the fit. Subsequently, employing this decay model with the fitted
parameters, we simulate results starting from the truncation position
$z_{t}=z_{c}-3$, which is the midpoint of the fit range, and extend towards
long distances where the matrix element approaches zero. We constrain the fit
parameters with the priors $\\{A,m,d\\}>0$ and extrapolate the data to
$z/a=64$, which can lead to extrapolated results decaying to zero in most
cases. However, in the cases of $-t=[0,0.231]$ GeV2 for $P_{z}=1.937$ GeV, due
to the slower decay, we need to impose a stricter constraint $m>0.2$ GeV [70]
as a fit prior and extrapolate the matrix elements to a longer distance
$z/a=200$. The extrapolated results, indicated by the purple circle symbols in
Fig. 3, exhibit desired decay behavior at long distances and start to approach
zero as the distances become sufficiently large. Therefore, the extrapolation
can correct the artificial behavior of the lattice data at large $\lambda$.

Figure 3:  The results of the extrapolated together with the original
renormalized matrix elements obtained with NNLO+LRR+RGR ($\kappa=1$)
coefficients are shown as a function of $\lambda=zP_{z}$.

To obtain the $x$ dependence of the qGPDs, we need to integrate the
extrapolated renormalized matrix elements over all values of $\lambda=zP_{z}$
in the Fourier transform. In Fig. 4, we show the valence pion qGPD results at
$P_{z}=1.453$ and 1.937 GeV, which are obtained using NNLO+LRR+RGR
($\kappa=1.414$) renormalization. We present results for six $-t$ values
corresponding to each $P_{z}$ value. We see that as $-t$ increases, the
distributions start at smaller values and decay slower with increasing $x$. If
we compare the data with similar $-t$ sizes between these two panels, we can
find that there are considerable differences between $P_{z}=1.453$ GeV and
$1.937$ GeV results. We expect that the perturbative matching to light-cone
GPDs, discussed in the next section, can correct for a significant portion of
these differences.

Figure 4:  The pion valence quark qGPD results for all $-t$ values obtained
with NNLO+LRR+RGR ($\kappa=1.414$) coefficient are shown as a function of the
Bjorken-$x$. Left: $P_{z}=1.453$ GeV, right: $P_{z}=1.937$ GeV.

##  III  Numerical results on the valence pion light-cone GPDs

We utilize the LaMET approach for the perturbative matching of qGPDs to the
light-cone GPDs

$\displaystyle H(x,\mu,t)$ $\displaystyle=\int_{-\infty}^{\infty}\frac{{\rm
d}k}{|k|}\int_{-\infty}^{\infty}\frac{{\rm d}y}{|y|}{\cal C}_{\rm
evo}^{-1}\left(\frac{x}{k},\frac{\mu}{\mu_{0}}\right)\times{\cal
C}^{-1}\left(\frac{k}{y},\frac{\mu_{0}}{yP_{z}},|y|\lambda_{s}\right)\widetilde{H}(y,P_{z},t,z_{s},\mu_{0})$
(20) $\displaystyle\approx\sum_{i,j}\frac{\Delta k}{|k_{i}|}\frac{\Delta
y}{|y_{j}|}{\cal C}_{\rm
evo}^{-1}\left(\frac{x}{k_{i}},\frac{\mu}{\mu_{0}}\right)\times{\cal
C}^{-1}\left(\frac{k_{i}}{y_{j}},\frac{\mu_{0}}{y_{j}P_{z}},|y_{j}|\lambda_{s}\right)\widetilde{H}(y_{j},P_{z},t,z_{s},\mu_{0}),$

where $\lambda_{s}=z_{s}P_{z}$, ${\cal C}^{-1}_{\rm evo}$ and ${\cal C}^{-1}$
are the inverses of the DGLAP evolution kernel and the hybrid-scheme matching
kernel, respectively. The detailed expressions for these kernels up to the
NNLL accuracy are given in Appendix B.

Since the qGPDs vanish at the positive and negative infinities, the
integration range is finite in practice. To obtain numerical results on the
light-cone GPDs, we can discretize the integration variables and replace the
integrals with sums, as shown in the second line of Eq. (20). This procedure
can reproduce the exact numerical integration well [70]. In Fig. 5, we
demonstrate both qGPD and light-cone GPD results obtained from NNLO+LRR+RGR
($\kappa=1.414$) perturbative matching with three-loop DGLAP evolution [81,
82]. It shows the results from two datasets with different values of $P_{z}$
but similar values of $-t$. As seen from this figure, there is a significant
difference in the $x$ dependence between the qGPDs and light-cone GPDs,
especially at lower $P_{z}$. The effects of the matching are most significant
at small and large values of $x$. Additionally, the $P_{z}$ dependence of the
distributions is significantly reduced after matching, as expected. The slight
dip in the solid-filled green band is likely a numerical artifact, which
should be smooth judging by the smoothness of results nearby. In this figure,
we only show the results for the light-cone GPDs for $x>0.15$ for the reasons
that will be explained below.

Since a larger $P_{z}$ is favored for better control of power accuracy, we
will focus on the case with the largest momentum $P_{z}=1.937$ GeV in the
following discussion. We will explore the dependence of the light-cone GPDs on
the matching accuracy, renormalization scale, and physical variables.

Figure 5:  The $x$ dependence of the qGPDs and light-cone (LC) GPDs, solved
with NNLO+LRR+RGR ($\kappa=1.414$) perturbative matching, is shown. It
contains the results from two datasets with different values of $P_{z}$ but
similar values of $-t$.

###  III.1 The dependence of the light-cone GPDs on the perturbative accuracy
of the matching and the renormalization scale

In this subsection, we examine the sensitivity of the pion light-cone GPDs to
the perturbative accuracy in the matching and the renormalization scale,
aiming to determine the reliability of our results.

In Fig. 6, we present the light-cone GPD results for the largest momentum
$P_{z}=1.937$ GeV, obtained using NLO and NNLO perturbative matching with LRR,
both with or without RGR. The line-filled bands represent the results obtained
without RGR at $\mu=2$ GeV. The solid-filled bands show the results obtained
with RGR, considering the variation of $\kappa$. The darker-solid-filled bands
display the statistical errors with $\kappa=1.414$ only, while the lighter-
solid-filled bands also include additional systematic errors associated with
scale variation from $\kappa=1$ to 2. The scale variation of the results is
more significant at small and large $x$. We have selected the data with the
smallest and largest momentum transfers for a clearer comparison. The NNLO
results exhibit smaller scale variation, particularly for data with larger
$-t$ values. Moreover, the results with and without RGR are almost consistent
in the intermediate $x$ region. However, they have a significant difference
for $x<0.2$. Since at low scales, the strong coupling constant and the
logarithms diverge in the case of RGR, causing the perturbation theory to
break down. Consequently, our light-cone GPD results obtained with RGR are
only reliable for $x>0.2$. The RGR also plays a role in the vicinity of $x=1$.
In this region, the threshold logs may also be important [55, 83]. However,
threshold resummation is beyond the scope of this work. Finally, we note that
the results between the NLO and NNLO display convergence, especially for the
case with larger $-t$ value.

Figure 6:  The $x$ dependence of the valence light-cone GPDs for $P_{z}=1.937$
GeV is presented, showing results obtained with various perturbative matching
choices and renormalization scales. The line-filled bands correspond to the
results without RGR, while the solid-filled bands correspond to the results
with RGR. For the latter, the darker-solid-filled bands represent the
statistical uncertainty with $\kappa=1.414$, while the lighter-solid-filled
bands contain additional systematic uncertainty with varying $\kappa$ from 1
to 2.

If one uses the hybrid scheme to perform the renormalization, the resulting
light-cone GPDs may depend on the choice of $z_{s}$. Therefore, it is
important to check the dependence of our results on $z_{s}$. We compare the
valence light-cone GPD results with three different values of $z_{s}$ for
$P_{z}=1.937$ GeV in Fig. 7. As a demonstration, we show our results for
$-t=0$ GeV2, obtained using NNLO+LRR+RGR. The lines represent the central
values obtained with $\kappa=1.414$, while the bands encompass both
statistical and systematic errors. This convention for the error bands is used
in the subsequent figures as well. We can see that the final valence light-
cone GPD results have little dependence on $z_{s}$. The situations are similar
for other values of $-t$. Thus, our findings are not sensitive to the choice
of $z_{s}$. Since the value of $\bar{m}_{0}$ at $z/a=3$ shows the smallest
dependence on the perturbative order and the renormalization scale, we adopt
$z_{s}/a=3$ for the following analyses.

Furthermore, the zero momentum transfer case corresponds to the PDF results,
allowing for comparison with previous PDF calculations. Specifically, we
compare our findings with the global analysis at NLO accuracy (JAM21) [84] and
our previous results using the NNLO matching (without LRR or RGR) within the
same lattice framework (BNL-ANL21) [70]. Our results agree well with the
global fit results for $0.2\lesssim x\lesssim 0.8$. Regarding the comparison
with BNL-ANL21, differences for $x<0.2$ are not surprising, as discussed
earlier. Additionally, we observe that the results obtained with LRR and RGR
in this work show a slight downward shift compared to previous findings in the
intermediate $x$ region, while they slightly exceed previous results for
$x\lesssim 0.25$.

Figure 7:  The dependence of the valence light-cone GPDs on $z_{s}$ is
illustrated at $P_{z}=1.937$ GeV and $-t=0$ GeV2. The black band displays the
NLO global analysis results [84]. The red band represents the valence pion
PDFs obtained with NNLO matching in a hybrid scheme, excluding LRR and RGR
[70].

###  III.2 The $t$-dependence of the light-cone GPDs

Having analyzed the systematic effects in the determination of the valence
pion light-cone GPDs at some representative values of $-t$, we can now study
the entire $t$ dependence of the light-cone GPDs in more detail. In Fig. 8, we
present the valence light-cone GPDs as a function of the Bjorken-$x$ across
all the available values of $-t$. As expected from the previous discussion,
for $x<0.2$, there are significant differences in all results between the
cases with (solid-filled bands) and without (line-filled bands) RGR. The
distributions exhibit significant dependence on $-t$ and $x$. Firstly, the
light-cone GPDs $H(x,\xi=0,t)$ decrease as $-t$ increases for a fixed $x$.
Secondly, the fall-off of the GPDs along $x$ is notably reduced at larger
values of $-t$. Similar behavior of the pion GPDs was also found in Ref. [32],
where a lattice spacing of $a=0.09$ fm and a symmetric frame were utilized.

Figure 8:  The valence pion light-cone GPDs as a function of $x$ for all
available values of $-t$. The solid-filled bands represent results obtained
with NNLO+LRR+RGR, and the width of the bands indicates the total uncertainty
of the results, including scale variation. The line-filled bands represent
NNLO+LRR results for $\mu=2$ GeV.

Figure 9:  Left: The $t$-dependence of the valence light-cone GPDs for various
values of $x$. Right: The pion effective radius as a function of the
Bjorken-$x$.

To further explore the $t$ and $x$ dependence of the pion valence light-cone
GPDs, we would like to parametrize our results with a suitable ansatz. Based
on past experience with the pion form factor [71], we employ the monopole form

$H(x,t)=\frac{H(x,0)}{1-t/M^{2}(x)},$ (21)

where $H(x,0)$ represents the GPDs at zero momentum transfer, equivalent to
the PDFs, and $M(x)$ is a $x$-dependent monopole mass parameter, which is
related to the pion effective radius as $\langle r^{2}(x)\rangle=6/M^{2}(x)$,
in analogy to the usual pion charge radius. This fit ansatz works very well,
as demonstrated by the bands in the left panel of Fig. 9. The left panel
displays the $t$ dependence of the light-cone GPDs for four selected values of
$x$, namely $x=0.2,~{}0.3,~{}0.5$, and 0.8. The lattice results are shown as
the data points, and the corresponding monopole fit results are displayed as
the bands. We can observe that the $t$ dependence of the valence light-cone
GPDs becomes milder as $x$ increases. For $x=0.8$, almost no $t$ dependence
can be discerned within the estimated errors. The results of the pion
effective radius are shown in the right panel of Fig. 9. The effective radius
clearly decreases with increasing $x$. This means that when quarks have higher
momentum fractions, they are likely to be confined to a smaller spatial region
in the transverse plane. This confinement results in narrower distributions in
position space, potentially reducing the effective radius. Interestingly, the
effective radius for $x=0.2$ is comparable to the pion charge radius obtained
from the pion electromagnetic form factor, specifically $\langle
r_{\pi}^{2}\rangle=0.313(27)$ fm2, as calculated in Ref. [71] using the same
pion mass.

Moreover, the Fourier transform of the valence light-cone GPDs at zero
skewness, with respect to the transverse components of the momentum transfer,
defines the valence impact-parameter-space parton distributions (IPDs):

$q(x,\mathbf{b}_{\perp})=\int\frac{{\rm
d}^{2}\mathbf{\Delta}_{\perp}}{(2\pi)^{2}}H(x,\mathbf{\Delta}^{2}_{\perp})e^{i\mathbf{b}_{\perp}\cdot\mathbf{\Delta}_{\perp}},$
(22)

where $\mathbf{b}_{\perp}$ is the impact parameter. The IPDs describe the
probability density of finding a parton with momentum fraction $x$ at a
specific transverse distance $\mathbf{b}_{\perp}$ from the center of the
transverse momentum (CoTM). They can offer a comprehensive view of the parton
distributions in both momentum and coordinate spaces within the hadron. Using
our monopole fit results of $H(x,t)$, we can perform a Fourier transform in
$\mathbf{\Delta}_{\perp}$ to obtain the valence IPDs. Fig. 10 shows the
results with varying $x$ from 0.2 to 0.9 and $\mathbf{b}_{\perp}$ from 0 to 1
fm. We can find that when the partons carry a larger momentum fraction $x$,
the distributions are more localized, which results in a more concentrated
distribution in the impact-parameter space. Conversely, as $x$ decreases, the
distributions become broader in $\mathbf{b}_{\perp}$, indicating a more
diffuse spatial distribution of lower-momentum quarks.

Figure 10:  The valence IPDs are shown as functions of the Bjorken-$x$ and the
impact parameter $\mathbf{b}_{\perp}$.

It is interesting to compare the pion valence IPDs with those of the proton,
as shown in Fig. 11. We display the comparison in a two-dimensional
distribution at $x=0.3$ and 0.5. The light-cone proton GPDs data are from Ref.
[85]. We perform a dipole fit on the proton GPDs and then get the
corresponding IPDs through the Fourier transform. For both the pion and
proton, the distributions become more concentrated at larger $x$,
corresponding to smaller effective radii, which is consistent with our
findings from the right panel of Fig. 9. Comparison of the proton and pion
results at fixed $x$ reveals that the proton exhibits a broader distribution
of valence quarks than the pion at both $x=0.3$ and 0.5. To some extent, this
observation relates to the fact that the charge radius of the proton is larger
than that of the pion. Fig. 11 also indicates that the valence quark
distributions within the pion are somewhat less peaked at the very center
compared to those within the proton.

Figure 11:  The two-dimensional distributions show a comparison of the IPDs
between the pion and proton [85] at $x$=0.3 and 0.5.

##  IV  Conclusions

We study the pion valence light-cone GPDs using the LaMET approach at a
lattice spacing of 0.04 fm, which is more than two times smaller than the
lattice spacing used in other studies of pion GPDs ($a=0.09$ fm) [32]. In this
work, we utilize the Lorentz-invariant definition of the GPDs, enabling us to
perform calculations in the non-Breit frame. This approach enables
simultaneous calculation of multiple momentum transfers and thus significantly
reduces the computational costs. In this work, we perform calculations within
both the Breit and non-Breit frames. By comparing the amplitude results of
these two frames, we demonstrated once again the validity and efficacy of the
Lorentz-invariant approach. To obtain the pion valence light-cone GPDs from
the qGPDs, we use NNLO hybrid-scheme matching along with leading renormalon
resummation and renormalization group resummation. The use of leading
renormalon resummation significantly reduces the perturbative uncertainty in
the hybrid renormalization scheme. Furthermore, the use of renormalization
group resummation allows us to determine the range of Bjorken-$x$ for which
the LaMET approach is reliable for a specific momentum value $P_{z}$. Our
results at zero momentum transfer, specifically the valence light-cone PDF
results, are in good agreement with the global analyses [84] within the range
$0.2\lesssim x\lesssim 0.8$. The $P_{z}$ dependence of the final light-cone
GPD results is significantly reduced compared to the qGPDs, which indicates
the effectiveness of the perturbative matching framework. Our lattice
calculations provide a comprehensive three-dimensional imaging of the pion
structure, which can be better visualized in terms of impact-parameter-
dependent distribution. From these findings, we observe that the effective
transverse size of the pion reduces as $x$ increases, a pattern also observed
in the proton [85]. Additionally, the effective size of the pion is
consistently smaller than that of the proton for the considered $x$ region.

###### Acknowledgements.

We thank Shohini Bhattacharya, Yushan Su, and Rui Zhang for their valuable
communications. This material is based upon work supported by the U.S.
Department of Energy, Office of Science, Office of Nuclear Physics through
Contract Nos. DE-SC0012704, DE-AC02-06CH11357, and within the frameworks of
Scientific Discovery through Advanced Computing (SciDAC) award Fundamental
Nuclear Physics at the Exascale and Beyond, and under the umbrella of the
Quark-Gluon Tomography (QGT) Topical Collaboration with Award DE-SC0023646.
HTD is supported by the NSFC under grant No. 12325508 and the National Key
Research and Development Program of China under Grant No. 2022YFA1604900. SS
is supported by the National Science Foundation under CAREER Award
PHY-1847893. YZ is partially supported by the 2023 Physical Sciences and
Engineering (PSE) Early Investigator Named Award program at Argonne National
Laboratory. This research used awards of computer time provided by the U.S.
Department of Energy’s INCITE and ALCC programs at the Argonne and the Oak
Ridge Leadership Computing Facilities. The Argonne Leadership Computing
Facility at Argonne National Laboratory is supported by the Office of Science
of the U.S. DOE under Contract No. DE-AC02-06CH11357. The Oak Ridge Leadership
Computing Facility at Oak Ridge National Laboratory is supported by the Office
of Science of the U.S. DOE under Contract No. DE-AC05-00OR22725. This research
also used the Delta advanced computing and data resource which is supported by
the National Science Foundation (award OAC 2005572) and the State of Illinois.
Delta is a joint effort of the University of Illinois Urbana-Champaign and its
National Center for Supercomputing Applications. Computations for this work
were carried out in part on facilities of the USQCD Collaboration, funded by
the Office of Science of the U.S. Department of Energy.

## References

  * Amendolia  _et al._ [1986] S. R. Amendolia _et al._ (NA7), Nucl. Phys. B 277, 168 (1986).
  * Conway  _et al._ [1989] J. S. Conway _et al._ , Phys. Rev. D 39, 92 (1989).
  * Müller  _et al._ [1994] D. Müller, D. Robaschik, B. Geyer, F. M. Dittes, and J. Hořejši, Fortsch. Phys. 42, 101 (1994), arXiv:hep-ph/9812448 .
  * Ji [1997a] X.-D. Ji, Phys. Rev. Lett. 78, 610 (1997a), arXiv:hep-ph/9603249 .
  * Radyushkin [1996a] A. V. Radyushkin, Phys. Lett. B 380, 417 (1996a), arXiv:hep-ph/9604317 .
  * Ji [1997b] X.-D. Ji, Phys. Rev. D 55, 7114 (1997b), arXiv:hep-ph/9609381 .
  * Radyushkin [1996b] A. V. Radyushkin, Phys. Lett. B 385, 333 (1996b), arXiv:hep-ph/9605431 .
  * Collins  _et al._ [1997] J. C. Collins, L. Frankfurt, and M. Strikman, Phys. Rev. D 56, 2982 (1997), arXiv:hep-ph/9611433 .
  * Dudek  _et al._ [2012] J. Dudek _et al._ , Eur. Phys. J. A 48, 187 (2012), arXiv:1208.1244 [hep-ex] .
  * Adams  _et al._ [2018] B. Adams _et al._ , Letter of Intent: A New QCD facility at the M2 beam line of the CERN SPS (COMPASS++/AMBER) (2018), arXiv:1808.00848 [hep-ex] .
  * Abdul Khalek  _et al._ [2022] R. Abdul Khalek _et al._ , Nucl. Phys. A 1026, 122447 (2022), arXiv:2103.05419 [physics.ins-det] .
  * Anderle  _et al._ [2021] D. P. Anderle _et al._ , Front. Phys. (Beijing) 16, 64701 (2021), arXiv:2102.09222 [nucl-ex] .
  * Ji [2013] X. Ji, Phys. Rev. Lett. 110, 262002 (2013), arXiv:1305.1539 [hep-ph] .
  * Ji [2014] X. Ji, Sci. China Phys. Mech. Astron. 57, 1407 (2014), arXiv:1404.6680 [hep-ph] .
  * Ji  _et al._ [2021a] X. Ji, Y.-S. Liu, Y. Liu, J.-H. Zhang, and Y. Zhao, Rev. Mod. Phys. 93, 035005 (2021a), arXiv:2004.03543 [hep-ph] .
  * Hagler  _et al._ [2003] P. Hagler, J. W. Negele, D. B. Renner, W. Schroers, T. Lippert, and K. Schilling (LHPC, SESAM), Phys. Rev. D 68, 034505 (2003), arXiv:hep-lat/0304018 .
  * Gockeler  _et al._ [2004] M. Gockeler, R. Horsley, D. Pleiter, P. E. L. Rakow, A. Schafer, G. Schierholz, and W. Schroers (QCDSF), Phys. Rev. Lett. 92, 042002 (2004), arXiv:hep-ph/0304249 .
  * Gockeler  _et al._ [2005] M. Gockeler, P. Hagler, R. Horsley, D. Pleiter, P. E. L. Rakow, A. Schafer, G. Schierholz, and J. M. Zanotti (QCDSF, UKQCD), Phys. Lett. B 627, 113 (2005), arXiv:hep-lat/0507001 .
  * Brommel  _et al._ [2007] D. Brommel _et al._ (QCDSF-UKQCD), PoS LATTICE2007, 158 (2007), arXiv:0710.1534 [hep-lat] .
  * Hagler  _et al._ [2008] P. Hagler _et al._ (LHPC), Phys. Rev. D 77, 094502 (2008), arXiv:0705.4295 [hep-lat] .
  * Alexandrou  _et al._ [2011] C. Alexandrou, J. Carbonell, M. Constantinou, P. A. Harraud, P. Guichon, K. Jansen, C. Kallidonis, T. Korzec, and M. Papinutto, Phys. Rev. D 83, 114513 (2011), arXiv:1104.1600 [hep-lat] .
  * Alexandrou  _et al._ [2013] C. Alexandrou, M. Constantinou, S. Dinter, V. Drach, K. Jansen, C. Kallidonis, and G. Koutsou, Phys. Rev. D 88, 014509 (2013), arXiv:1303.5979 [hep-lat] .
  * Bali  _et al._ [2019] G. S. Bali, S. Collins, M. Göckeler, R. Rödl, A. Schäfer, and A. Sternbeck, Phys. Rev. D 100, 014507 (2019), arXiv:1812.08256 [hep-lat] .
  * Alexandrou  _et al._ [2020a] C. Alexandrou _et al._ , Phys. Rev. D 101, 034519 (2020a), arXiv:1908.10706 [hep-lat] .
  * Alexandrou  _et al._ [2023] C. Alexandrou _et al._ , Phys. Rev. D 107, 054504 (2023), arXiv:2202.09871 [hep-lat] .
  * Chen  _et al._ [2020] J.-W. Chen, H.-W. Lin, and J.-H. Zhang, Nucl. Phys. B 952, 114940 (2020), arXiv:1904.12376 [hep-lat] .
  * Alexandrou  _et al._ [2020b] C. Alexandrou, K. Cichy, M. Constantinou, K. Hadjiyiannakou, K. Jansen, A. Scapellato, and F. Steffens, Phys. Rev. Lett. 125, 262001 (2020b), arXiv:2008.10573 [hep-lat] .
  * Lin [2021] H.-W. Lin, Phys. Rev. Lett. 127, 182001 (2021), arXiv:2008.12474 [hep-ph] .
  * Lin [2022] H.-W. Lin, Phys. Lett. B 824, 136821 (2022), arXiv:2112.07519 [hep-lat] .
  * 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] .
  * Bhattacharya  _et al._ [2023] S. Bhattacharya, K. Cichy, M. Constantinou, X. Gao, A. Metz, J. Miller, S. Mukherjee, P. Petreczky, F. Steffens, and Y. Zhao, Phys. Rev. D 108, 014507 (2023), arXiv:2305.11117 [hep-lat] .
  * Lin [2023] H.-W. Lin, Phys. Lett. B 846, 138181 (2023), arXiv:2310.10579 [hep-lat] .
  * Bhattacharya  _et al._ [2024a] S. Bhattacharya _et al._ , Phys. Rev. D 109, 034508 (2024a), arXiv:2310.13114 [hep-lat] .
  * Holligan and Lin [2023] J. Holligan and H.-W. Lin, (2023), arXiv:2312.10829 [hep-lat] .
  * Constantinou  _et al._ [2021] M. Constantinou _et al._ , Prog. Part. Nucl. Phys. 121, 103908 (2021), arXiv:2006.08636 [hep-ph] .
  * Hannaford-Gunn  _et al._ [2022] A. Hannaford-Gunn, K. U. Can, R. Horsley, Y. Nakamura, H. Perlt, P. E. L. Rakow, H. Stüben, G. Schierholz, R. D. Young, and J. M. Zanotti (CSSM/QCDSF/UKQCD), Phys. Rev. D 105, 014502 (2022), arXiv:2110.11532 [hep-lat] .
  * Hannaford-Gunn  _et al._ [2024] A. Hannaford-Gunn, K. U. Can, J. A. Crawford, R. Horsley, P. E. L. Rakow, G. Schierholz, H. Stüben, R. D. Young, and J. M. Zanotti, (2024), arXiv:2405.06256 [hep-lat] .
  * Dutrieux  _et al._ [2024] H. Dutrieux, R. Edwards, C. Egerer, J. Karpie, C. Monahan, K. Orginos, A. Radyushkin, D. Richards, E. Romero, and S. Zafeiropoulos, (2024), arXiv:2405.10304 [hep-lat] .
  * Bhattacharya  _et al._ [2024b] S. Bhattacharya, K. Cichy, M. Constantinou, A. Metz, N. Nurminen, and F. Steffens, (2024b), arXiv:2405.04414 [hep-lat] .
  * Constantinou and Panagopoulos [2017] M. Constantinou and H. Panagopoulos, Phys. Rev. D 96, 054506 (2017), arXiv:1705.11193 [hep-lat] .
  * Alexandrou  _et al._ [2017] C. Alexandrou, K. Cichy, M. Constantinou, K. Hadjiyiannakou, K. Jansen, H. Panagopoulos, and F. Steffens, Nucl. Phys. B 923, 394 (2017), arXiv:1706.00265 [hep-lat] .
  * Stewart and Zhao [2018] I. W. Stewart and Y. Zhao, Phys. Rev. D 97, 054512 (2018), arXiv:1709.04933 [hep-ph] .
  * Chen  _et al._ [2018] J.-W. Chen, T. Ishikawa, L. Jin, H.-W. Lin, Y.-B. Yang, J.-H. Zhang, and Y. Zhao, Phys. Rev. D 97, 014505 (2018), arXiv:1706.01295 [hep-lat] .
  * Radyushkin [2017] A. V. Radyushkin, Phys. Rev. D 96, 034025 (2017), arXiv:1705.01488 [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] .
  * Braun  _et al._ [2019] V. M. Braun, A. Vladimirov, and J.-H. Zhang, Phys. Rev. D 99, 014013 (2019), arXiv:1810.00048 [hep-ph] .
  * Li  _et al._ [2021] Z.-Y. Li, Y.-Q. Ma, and J.-W. Qiu, Phys. Rev. Lett. 126, 072001 (2021), arXiv:2006.12370 [hep-ph] .
  * Fan  _et al._ [2020] Z. Fan, X. Gao, R. Li, H.-W. Lin, N. Karthik, S. Mukherjee, P. Petreczky, S. Syritsyn, Y.-B. Yang, and R. Zhang, Phys. Rev. D 102, 074504 (2020), arXiv:2005.12015 [hep-lat] .
  * Ji  _et al._ [2021b] X. Ji, Y. Liu, A. Schäfer, W. Wang, Y.-B. Yang, J.-H. Zhang, and Y. Zhao, Nucl. Phys. B 964, 115311 (2021b), arXiv:2008.03886 [hep-ph] .
  * Ji  _et al._ [2015] X. Ji, A. Schäfer, X. Xiong, and J.-H. Zhang, Phys. Rev. D 92, 014039 (2015), arXiv:1506.00248 [hep-ph] .
  * Liu  _et al._ [2019] Y.-S. Liu, W. Wang, J. Xu, Q.-A. Zhang, J.-H. Zhang, S. Zhao, and Y. Zhao, Phys. Rev. D 100, 034006 (2019), arXiv:1902.00307 [hep-ph] .
  * Ma  _et al._ [2022] J. P. Ma, Z. Y. Pang, and G. P. Zhang, JHEP 08, 130, arXiv:2202.07116 [hep-ph] .
  * Yao  _et al._ [2023] F. Yao, Y. Ji, and J.-H. Zhang, JHEP 11, 021, arXiv:2212.14415 [hep-ph] .
  * Chen  _et al._ [2021] L.-B. Chen, W. Wang, and R. Zhu, Phys. Rev. Lett. 126, 072002 (2021), arXiv:2006.14825 [hep-ph] .
  * Gao  _et al._ [2021a] X. Gao, K. Lee, S. Mukherjee, C. Shugert, and Y. Zhao, Phys. Rev. D 103, 094504 (2021a), arXiv:2102.01101 [hep-ph] .
  * Su  _et al._ [2023] Y. Su, J. Holligan, X. Ji, F. Yao, J.-H. Zhang, and R. Zhang, Nucl. Phys. B 991, 116201 (2023), arXiv:2209.01236 [hep-ph] .
  * Holligan  _et al._ [2023] J. Holligan, X. Ji, H.-W. Lin, Y. Su, and R. Zhang, Nucl. Phys. B 993, 116282 (2023), arXiv:2301.10372 [hep-lat] .
  * Zhang  _et al._ [2023] R. Zhang, J. Holligan, X. Ji, and Y. Su, Phys. Lett. B 844, 138081 (2023), arXiv:2305.05212 [hep-lat] .
  * Bazavov  _et al._ [2019] A. Bazavov _et al._ , Phys. Rev. D 100, 094510 (2019), arXiv:1908.09552 [hep-lat] .
  * Follana  _et al._ [2007] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong (HPQCD, UKQCD), Phys. Rev. D 75, 054502 (2007), arXiv:hep-lat/0610092 .
  * Hasenfratz and Knechtli [2001] A. Hasenfratz and F. Knechtli, Phys. Rev. D 64, 034504 (2001), arXiv:hep-lat/0103029 .
  * Gao  _et al._ [2020] X. Gao, L. Jin, C. Kallidonis, N. Karthik, S. Mukherjee, P. Petreczky, C. Shugert, S. Syritsyn, and Y. Zhao, Phys. Rev. D 102, 094513 (2020), arXiv:2007.06590 [hep-lat] .
  * Bali  _et al._ [2016] G. S. Bali, B. Lang, B. U. Musch, and A. Schäfer, Phys. Rev. D 93, 094515 (2016), arXiv:1602.05525 [hep-lat] .
  * Brannick  _et al._ [2008] J. Brannick, R. C. Brower, M. A. Clark, J. C. Osborn, and C. Rebbi, Phys. Rev. Lett. 100, 041601 (2008), arXiv:0707.4018 [hep-lat] .
  * Clark  _et al._ [2010] M. A. Clark, R. Babich, K. Barros, R. C. Brower, and C. Rebbi (QUDA), Comput. Phys. Commun. 181, 1517 (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 (QUDA) (2011) arXiv:1109.2935 [hep-lat] .
  * Clark  _et al._ [2016] M. A. Clark, B. Joó, A. Strelchenko, M. Cheng, A. Gambhir, and R. C. Brower (QUDA) (2016) arXiv:1612.07873 [hep-lat] .
  * Izubuchi  _et al._ [2019] T. Izubuchi, L. Jin, C. Kallidonis, N. Karthik, S. Mukherjee, P. Petreczky, C. Shugert, and S. Syritsyn, Phys. Rev. D 100, 034516 (2019), arXiv:1905.06349 [hep-lat] .
  * Shintani  _et al._ [2015] E. Shintani, R. Arthur, T. Blum, T. Izubuchi, C. Jung, and C. Lehner, Phys. Rev. D 91, 114511 (2015), arXiv:1402.0244 [hep-lat] .
  * Gao  _et al._ [2022a] X. Gao, A. D. Hanlon, S. Mukherjee, P. Petreczky, P. Scior, S. Syritsyn, and Y. Zhao, Phys. Rev. Lett. 128, 142003 (2022a), arXiv:2112.02208 [hep-lat] .
  * Gao  _et al._ [2021b] X. Gao, N. Karthik, S. Mukherjee, P. Petreczky, S. Syritsyn, and Y. Zhao, Phys. Rev. D 104, 114515 (2021b), arXiv:2102.06047 [hep-lat] .
  * Gao  _et al._ [2022b] X. Gao, A. D. Hanlon, N. Karthik, S. Mukherjee, P. Petreczky, P. Scior, S. Shi, S. Syritsyn, Y. Zhao, and K. Zhou, Phys. Rev. D 106, 114510 (2022b), arXiv:2208.02297 [hep-lat] .
  * Gao  _et al._ [2022c] X. Gao, A. D. Hanlon, N. Karthik, S. Mukherjee, P. Petreczky, P. Scior, S. Syritsyn, and Y. Zhao, Phys. Rev. D 106, 074505 (2022c), arXiv:2206.04084 [hep-lat] .
  * Ding  _et al._ [2024] H.-T. Ding, X. Gao, A. D. Hanlon, S. Mukherjee, P. Petreczky, Q. Shi, S. Syritsyn, R. Zhang, and Y. Zhao, (2024), arXiv:2404.04412 [hep-lat] .
  * Ji  _et al._ [2018] X. Ji, J.-H. Zhang, and Y. Zhao, Phys. Rev. Lett. 120, 112001 (2018), arXiv:1706.08962 [hep-ph] .
  * Ishikawa  _et al._ [2017] T. Ishikawa, Y.-Q. Ma, J.-W. Qiu, and S. Yoshida, Phys. Rev. D 96, 094019 (2017), arXiv:1707.03107 [hep-ph] .
  * Green  _et al._ [2018] J. Green, K. Jansen, and F. Steffens, Phys. Rev. Lett. 121, 022004 (2018), arXiv:1707.07152 [hep-lat] .
  * Huo  _et al._ [2021] Y.-K. Huo _et al._ (Lattice Parton Collaboration (LPC)), Nucl. Phys. B 969, 115443 (2021), arXiv:2103.02965 [hep-lat] .
  * Bazavov  _et al._ [2018] A. Bazavov, N. Brambilla, P. Petreczky, A. Vairo, and J. H. Weber (TUMQCD), Phys. Rev. D 98, 054511 (2018), arXiv:1804.10600 [hep-lat] .
  * Braun  _et al._ [2020] V. M. Braun, K. G. Chetyrkin, and B. A. Kniehl, JHEP 07, 161, arXiv:2004.01043 [hep-ph] .
  * Curci  _et al._ [1980] G. Curci, W. Furmanski, and R. Petronzio, Nucl. Phys. B 175, 27 (1980).
  * Moch  _et al._ [2004] S. Moch, J. A. M. Vermaseren, and A. Vogt, Nucl. Phys. B 688, 101 (2004), arXiv:hep-ph/0403192 .
  * Ji  _et al._ [2023] X. Ji, Y. Liu, and Y. Su, JHEP 08, 037, arXiv:2305.04416 [hep-ph] .
  * Barry  _et al._ [2021] P. Barry, C.-R. Ji, N. Sato, and W. Melnitchouk, Physical Review Letters 127, 10.1103/physrevlett.127.232001 (2021).
  * Cichy  _et al._ [2023] K. Cichy _et al._ , Acta Phys. Polon. Supp. 16, 7 (2023), arXiv:2304.14970 [hep-lat] .

##  Appendix A Fit results of the two-point and three-point correlation
functions

In this appendix, we present some detailed analysis of the pion two-point and
three-point correlation functions.

As mentioned in the main text, by employing the spectral decomposition formula
along with 2-state or 3-state fits of the pion two-point correlators [62, 70,
71, 72, 73], we extract energy values of both the ground and excited states of
the pion, as well as the corresponding amplitudes. We show the results of the
ground and the first excited state energy for different momenta in Fig. 12. We
can see that the fit results for $E_{0}$ agree well with the dispersion
relation shown by the solid lines, and in most cases, the results for $E_{1}$
also agree with them within errors.

Figure 12:  The results for the ground state energy $E_{0}$ and the first
excited state energy $E_{1}$, extracted from the two-point functions, are
shown. The solid lines represent the results calculated using the dispersion
relation $E_{n}=\sqrt{m^{2}_{n}+|\textbf{P}|^{2}}$ with $m_{0}$ = 0.3 GeV and
$m_{1}$ = 1.46 GeV.

To obtain the matrix elements, we consider the ratio $R^{\mu}$ given by Eq.
(12) and perform 2-state or 3-state fits. Here, we use the energy levels and
amplitudes obtained from the corresponding fits of the two-point functions. In
Fig. 13, we select two specific cases from the data of $z=0$ and $\mu=t$ to
show the results of the ratio as $R^{t}/\bar{P}^{t}$. The left one belongs to
the case of $P_{z}$ = 0.97 GeV and $-t$ = 0.95 GeV2, while the right one
corresponds to the case of $P_{z}$ = 1.94 GeV and $-t$ = 1.69 GeV2, that is,
the largest momentum transfer in this work. The lattice results are shown as
the data points with different colors representing the different time
separations. The bands with corresponding colors denote the two-state fit
results of the lattice data. And the grey bands are the extrapolated results
of the ratio, that is, the bare matrix elements. We collect all the bare
matrix element results of the largest momentum in Fig. 14. It is evident that
the results display a decreasing trend with increasing $z$ and also as the
momentum transfer increases, as expected.

Figure 13:  Two specific examples of the ratio $R^{t}(z=0)/\bar{P}^{t}$ are
shown. Left: $P_{z}$ = 0.97 GeV and $-t$ = 0.95 GeV2, right: $P_{z}$ = 1.94
GeV and $-t$ = 1.69 GeV2. The data points represent the lattice results, and
the bands with corresponding colors are the fit results of the lattice data.
The grey bands show the results of the bare matrix elements. Figure 14:  The
bare matrix element results for the largest momentum
$|\mathbf{P}^{f}|=P_{z}=1.937$ GeV are shown as a function of $z$.

##  Appendix B  Perturbative coefficient

###  B.1  Matching coefficient

The NNLO matching kernel with LRR and RGR can be expanded to ${\cal
O}(\alpha_{s})$ as [58, 57]

$\displaystyle{\cal
C}\left(\frac{x}{y},\frac{\mu_{0}}{yP_{z}},|y|\lambda_{s}\right)$
$\displaystyle=\delta\left(\frac{x}{y}-1\right)+\frac{\alpha_{s}}{2\pi}{\cal
C}^{(1)}\left(\frac{x}{y},\frac{\mu_{0}}{yP_{z}},|y|\lambda_{s}\right)+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}{\cal
C}^{(2)}\left(\frac{x}{y},\frac{\mu_{0}}{yP_{z}},|y|\lambda_{s}\right)$ (23)
$\displaystyle\quad+\frac{\mu_{0}}{yP_{z}}\delta{\cal C}_{0,{\rm
LRR\;rest}}{\cal C}_{\rm hybrid\;LRR}+{\cal O}(\alpha_{s}^{3}),$

with

$\displaystyle\alpha_{s}(\mu_{0})=\frac{4\pi
a_{s0}}{X+a_{s0}\frac{\beta_{1}}{\beta_{0}}\ln(X)+\left(a_{s0}\right)^{2}Y},\quad
a_{s0}=\frac{0.293}{4\pi},$ (24) $\displaystyle
X=1+a_{s0}\beta_{0}\ln\left(\left(\frac{\mu_{0}}{2}\right)^{2}\right),$
$\displaystyle
Y=\frac{\beta_{2}}{\beta_{0}}\left(1-\frac{1}{X}\right)+\left(\frac{\beta_{1}}{\beta_{0}}\right)^{2}\left(\frac{\ln(X)}{X}+\frac{1}{X}-1\right).$

Here, the QCD $\beta$-function

$\beta[\alpha_{s}(\mu)]=-2\alpha_{s}(\mu)\sum_{n=0}^{\infty}a_{s0}^{n+1}(\mu)\beta_{n}$
(25)

is given by

$\displaystyle\beta_{0}={11\over 3}C_{A}-{4\over 3}T_{F}n_{f},$ (26)
$\displaystyle\beta_{1}={34\over 3}C_{A}^{2}-\left({20\over
3}C_{A}+4C_{F}\right)T_{F}n_{f},$ (27) $\displaystyle\beta_{2}={2857\over
54}C_{A}^{3}+\left(2C_{F}^{2}-{205\over 9}C_{F}C_{A}-{1415\over
27}C_{A}^{2}\right)T_{F}n_{f}+\left({44\over 9}C_{F}+{158\over
27}C_{A}\right)T_{F}^{2}n_{f}^{2},$ (28)

where $C_{F}=4/3$, $C_{A}=3$, $T_{F}=1/2$, and $n_{f}=3$ for our lattice
ensemble. The 1-loop and 2-loop matching coefficients ${\cal C}^{(1)}$ and
${\cal C}^{(2)}$ are constructed as

$\displaystyle{\cal C}^{(1)}={\cal C}^{(1)}_{\rm LL}+{\cal C}^{(1)}_{\rm
hybrid\;LL}+{\cal C}^{(1)}_{{\rm
NLL}}\ln\left(\left(\frac{\mu_{0}}{2P_{z}}\right)^{2}\right),$ (29)
$\displaystyle{\cal C}^{(2)}={\cal C}^{(2)}_{\rm LL}+{\cal C}^{(2)}_{\rm
hybrid\;LL}+\left({\cal C}^{(2)}_{{\rm NLL}}+{\cal C}^{(2)}_{{\rm
hybrid\;NLL}}\right)\ln\left(\left(\frac{\mu_{0}}{2P_{z}}\right)^{2}\right)+{\cal
C}^{(2)}_{{\rm
NNLL}}\left[\ln\left(\left(\frac{\mu_{0}}{2P_{z}}\right)^{2}\right)\right]^{2},$
(30)

where ${\cal C}_{\rm LL/NLL/NNLL}$ represent the normal LL/NLL/NNLL terms,
${\cal C}_{\rm hybrid\;LL/NLL}$ denote the LL/NLL terms in the hybrid scheme.
In the last line of Eq. (23), ${\cal C}_{\rm hybrid\;LRR}$ is the LRR term in
the hybrid scheme [58, 57]

$\displaystyle{\cal C}_{\rm hybrid\;LRR}$
$\displaystyle=0.575\mu_{0}\Bigg{\\{}\frac{e^{-A}p_{z}(1+A+A^{2})}{\epsilon_{m}^{2}\pi}+\frac{1}{\pi}\Bigg{[}\frac{e^{-A}z_{s}\sin\left[Bz_{s}p_{z}\right]}{B}+\frac{e^{-A}p_{z}}{\left[\epsilon_{m}^{2}+p_{z}^{2}B^{2}\right]^{2}}$
(31)
$\displaystyle\times\Big{(}\cos\left[Bz_{s}p_{z}\right]\left(\epsilon_{m}^{2}-B^{2}p_{z}^{2}+\epsilon_{m}^{3}z_{s}+\epsilon_{m}p_{z}^{2}B^{2}z_{s}\right)-\sin\left[Bz_{s}p_{z}\right]\left(2\epsilon+B^{2}p_{z}^{2}z_{s}+\epsilon_{m}^{2}z_{s}\right)Bp_{z}\Big{)}\Bigg{]}\Bigg{\\}}_{+},$

where $e^{-A}=e^{-\epsilon_{m}z_{s}}$ is a small exponential decaying factor,
$B=1-\frac{x}{y}$, $p_{z}=yP_{z}$, and $\\{\cdots\\}_{+}$ means a plus
function. Additionally, $\delta{\cal C}_{0,{\rm LRR\;rest}}$ represents the
rest LRR term

$\displaystyle\delta{\cal C}_{0,{\rm LRR\;rest}}=-0.575\kappa$
$\displaystyle\Bigg{\\{}\alpha_{s}+\frac{\alpha_{s}^{2}\left[2\beta_{0}^{2}+\beta_{1}(1+c_{1})\right]}{4\pi\beta_{0}}+e^{-\frac{2\pi}{\alpha_{s}\beta_{0}}}\left[4\pi\beta_{0}+\alpha_{s}c_{1}\beta_{1}+\frac{\alpha_{s}^{2}\beta_{1}c_{2}(\beta_{1}-2\beta_{0}^{2})}{4\pi\beta_{0}}\right]$
(32)
$\displaystyle\times\frac{1}{2\beta_{0}^{2}}\Re\left[E_{1+\frac{\beta_{1}}{2\beta_{0}^{2}}}\left(-\frac{2\pi}{\alpha_{s}\beta_{0}}\right)\right]\Bigg{\\}},$

where $E_{n}(z)=\int_{1}^{\infty}e^{-zt}/t^{n}{\rm d}t$ is the generalized
exponential integral, $\Re$ denotes the real part, and

$\displaystyle
c_{1}=\frac{\beta_{1}^{2}-\beta_{0}\beta_{2}}{4b\beta_{0}^{4}},$ (33)
$\displaystyle
c_{2}=\frac{\beta_{1}^{4}+4\beta_{0}^{3}\beta_{1}\beta_{2}-2\beta_{0}\beta_{1}^{2}\beta_{2}+\beta_{0}^{2}(\beta_{2}^{2}-2\beta_{1}^{3})-2\beta_{0}^{4}\beta_{3}}{32b(b-1)\beta_{0}^{8}},$

with $b=\beta_{1}/(2\beta_{0}^{2})$.

###  B.2  Evolution coefficient

The NNLL matching requires the 3-loop DGLAP evolution kernel, which can be
expressed as

$\displaystyle{\cal C}_{\rm
evo}\left(a_{s},\frac{x}{y},\frac{\mu_{0}}{\mu}\right)=$
$\displaystyle\delta\left(\frac{x}{y}-1\right)+a_{s}t{\cal
P}^{(0)}_{qq}\left(\frac{x}{y}\right)$ (34)
$\displaystyle+a_{s}^{2}\left[t{\cal P}^{V(1)}_{qq}+{t^{2}\over 2}\left({\cal
P}^{(0)}_{qq}\otimes{\cal P}^{(0)}_{qq}+\beta_{0}{\cal
P}^{(0)}_{qq}\right)\right]\left(\frac{x}{y}\right)$
$\displaystyle+a_{s}^{3}\left[t{\cal P}_{qq}^{V(2)}+t^{2}\left({\beta_{1}\over
2}{\cal P}_{qq}^{(0)}+{\beta_{0}}{\cal P}_{qq}^{V(1)}+{\cal
P}_{qq}^{(0)}\otimes{\cal P}_{qq}^{V(1)}\right)\right.$
$\displaystyle+\left.{t^{3}\over 6}\\!\left(2\beta_{0}^{2}{\cal
P}_{qq}^{(0)}\\!+\\!3\beta_{0}{\cal P}_{qq}^{(0)}\\!\otimes\\!{\cal
P}_{qq}^{(0)}\\!+\\!{\cal P}_{qq}^{(0)}\\!\otimes\\!{\cal
P}_{qq}^{(0)}\\!\otimes\\!{\cal
P}_{qq}^{(0)}\right)\\!\right]\\!\left(\frac{x}{y}\right),$

where $a_{s}=\alpha_{s}(\mu_{0})/(4\pi)$, $t=\ln(\mu^{2}_{0}/\mu^{2})$, ${\cal
P}^{(0)}_{qq}$ is the 1-loop splitting function, and ${\cal P}^{V(1)}_{qq}$
and ${\cal P}^{V(2)}_{qq}$ are the 2-loop and 3-loop splitting functions for
the valence quark [81, 82]. They are given by

$\displaystyle{\cal
P}^{(0)}_{qq}\left(z\right)=2C_{F}\left[\frac{1+z^{2}}{1-z}\right]_{+},$ (35)
$\displaystyle{\cal P}^{V(1)}_{qq}\left(z\right)=4\left[{\cal
P}^{(1)}_{qq}(z)-{\cal P}^{(1)}_{q\bar{q}}(z)\right],$ $\displaystyle{\cal
P}^{(1)}_{qq}(z)=C_{F}\left[C_{F}P_{F}(z)+\frac{1}{2}C_{A}P_{G}(z)+n_{f}T_{F}P_{n_{f}}(z)\right],$
$\displaystyle{\cal P}^{(1)}_{q\bar{q}}(z)=C_{F}\left(C_{F}-{1\over
2}C_{A}\right)P_{A}(z),$

where the definitions of $P_{F}(z)$, $P_{G}(z)$, $P_{n_{f}}(z)$ and $P_{A}(z)$
can be found in Ref. [81]. As for ${\cal P}^{V(2)}_{qq}$, we take the
approximate solution ${\cal P}^{V(2)}_{qq}(z)=P_{\rm ns}^{(2)-}(z)$ in Eq.
(4.23) of Ref. [82].

Note that after the inverse matching, one obtains the matched GPDs
$H(x,\mu_{0}=2\kappa xP_{z})$ at a varying $\overline{\rm MS}$ scale in $x$.
To implement the DGLAP evolution, we set $\mu=2$ GeV and $\mu_{0}=2\kappa
xP_{z}$ and obtain the evolution kernel as a triangular matrix down to a
minimal value $x_{\rm min}\sim 0.15$. Then we invert the triangular matrix and
apply it to the matched GPDs $H(x,\mu_{0})$, which leads to $H(x,\mu)$ at the
fixed scale $\mu$ down to $x_{\rm min}\sim 0.15$.

We could also express the evolution kernel ${\cal C}_{\rm evo}$ as a
perturbative series in $a_{s}(\mu)$, but it gives a smaller scale variation as
$a_{s}$ remains small. Instead, we use ${\cal C}_{\rm evo}$ as a perturbative
series in $a_{s}(2\kappa xP_{z})$. The latter becomes large at small $x$,
which gives us a more conservative estimate of the scale variation
uncertainty.


#  SPEECH EMOTION RECOGNITION USING QUATERNION CONVOLUTIONAL NEURAL NETWORKS

###### Abstract

Although speech recognition has become a widespread technology, inferring
emotion from speech signals still remains a challenge. To address this
problem, this paper proposes a quaternion convolutional neural network (QCNN)
based speech emotion recognition (SER) model in which Mel-spectrogram features
of speech signals are encoded in an RGB quaternion domain. We show that our
QCNN based SER model outperforms other real-valued methods in the Ryerson
Audio-Visual Database of Emotional Speech and Song (RAVDESS, 8-classes)
dataset, achieving, to the best of our knowledge, state-of-the-art results.
The QCNN also achieves comparable results with the state-of-the-art methods in
the Interactive Emotional Dyadic Motion Capture (IEMOCAP 4-classes) and Berlin
EMO-DB (7-classes) datasets. Specifically, the model achieves an accuracy of
77.87%, 70.46%, and 88.78% for the RAVDESS, IEMOCAP, and EMO-DB datasets,
respectively. In addition, our results show that the quaternion unit structure
is better able to encode internal dependencies to reduce its model size
significantly compared to other methods.

Index Terms: Speech Emotion Recognition, Signal Processing, Quaternion Deep
Learning, Convolutional Neural Networks

##  1 Introduction

Speech emotion recognition (SER) is an active, yet challenging area of
research that has many important implications in technologies such as
automated healthcare, clinical trials, voice assistants, psychological
therapy, emergency responders, call centers, video games, robot-human
interactions, and more. Nonetheless, recognizing emotion from speech signals
is difficult due to many reasons, namely: qualitative properties of
identifying emotion, background noise, variable individual-specific
accentuation, weak representation of grammatical and semantic knowledge,
temporal and spectral characteristics of the signal due to fast or slow
speech, dynamical characteristics such as soft or loud voices, and multiple
voices [1, 2, 3]. In addition, if an application can learn high-level features
from speech for emotion classification, it is possible to replicate or imitate
emotion from a text-to-speech context, further improving the aforementioned
fields of impact [4].

In general, SER is a classification task in which we extract features from a
set of emotion-labeled speech signals and use the pair of feature and label to
train a classifier which is commonly a deep neural architecture. Features such
as MFCC, chromograms, spectral contrast features, and Tonnetz representations
have been used in previous neural based SER models [5, 6, 7]. Neural based SER
models usually leverage n-dimensional convolutional neural networks (CNNs)
[8], recurrent neural networks (RNNs) [9], Long short-term memory networks
(LSTMs) [10], or as combinations and variations of these techniques [11].

Recently, the quaternion based neural networks have been shown to better
encode features than the standard real-valued based neural networks [12]. The
use of the quaternion domain has shown a unique ability in color image
processing to capture all three channels of the RGB domain without the typical
information loss of separating these channels, while reducing the size of its
model in comparison to real-valued models [13].

In this paper, to the best of our knowledge, we propose the first quaternion
based model that is able to encode an RGB representation of Mel-spectrograms
for the application of SER classification. Using this approach, we propose a
quaternion convolutional neural network (QCNN) to train on benchmark SER
datasets in the quaternion domain with quaternion converted elements of a CNN.
The results, conducted on three public SER datasets, show that the QCNN model
is able to yield state-of-the-art results on RAVDESS and outperforms all but
one competitive models in both IEMOCAP and EMO-DB in terms of accuracy. In
addition, our proposed model has the smallest size compared to the previous
competitive architectures.

The rest of this paper is organized as follows: In Section 2 we mention some
related works. In Section 3, we give a brief overview of quaternion operations
and describe the proposed model. In Section 4 we detail our experiments and
results. Finally, we draw the conclusion and give the future directions in
Section 5.

##  2 Relation to prior works

A majority of previous literature on SER methods consist of variations CNNs
such as a 1d-CNN for waveform features, CNN-LSTMs, transfer CNN models [14].
In addition, many methods have also focused on RNNs [15, 16] to discover new
features. In the context of quaternion models, [12] introduced the ability to
generate quaternion-based dense layers. Extending this work, QCNNs were
developed for color image processing [13]. The significance of quaternion
models in the context of signal processing was originally shown by [17] for
automatic speech recognition; however, this model encoded views of a time-
frame frequency as quaternions. Our work extends previous literature in that
we exploit the decibel visualization of Mel-spectrograms as RGB arrays to
encode in the quaternion domain rather than encoding the waveform features as
pure quaternions or using a real-valued model.

Figure 1:  The happy waveform converted to an RGB Mel-Spectrogram

##  3 Model

###  3.1 Preprocessing

Our feature generation pipeline consists of two parts: Mel-spectrogram array
generation and RGB conversion. After the wav files were split into emotion-
labelled folders, the Fourier transform was computed on each of the speech
waveforms to transform from the time domain to the frequency domain. While the
time axis scaled linearly by seconds, the frequency axis was non-linearly
transformed using a Mel transform, which is defined as:

$F_{\text{mel}}=2,595\log_{10}(1+\frac{f}{700})[dB]$ (1)

where $f$ denotes the frequency axis. The spectral amplitudes of the arrays
were normalized between their dataset specific Max[dB] and Min [dB]. The mel-
spectrogram was then converted to an image, with a color scale to represent
the amplitudes (see example given in Figure 1). The RGB images were then split
into training ( 80% ) and test (20%) datasets.

###  3.2 Quaternion Convolutional Neural Networks

####  3.2.1 Quaternion Algebra

Briefly, a quaternion is hypercomplex number which is an extension to complex
numbers in 3-dimensional space with additionally imaginary parts
$\boldsymbol{j}$ and $\boldsymbol{k}$. Thus, a quaternion $Q$ in domain $Q$
can be represented as
$Q=r_{0}+r_{1}\boldsymbol{i}+r_{2}\boldsymbol{j}+r_{3}\boldsymbol{k}$ where
$r\in R$ and $\boldsymbol{i},\boldsymbol{j}$, and $\boldsymbol{k}$ are
imaginary parts. Like complex number, Quaternions obey specific rules:

$i^{2}=j^{2}=k^{2}=ijk=-1$ (2)

and the adddition is defined as

$\begin{split}Qa+Qb&=(r_{a_{0}}+r_{b_{0}})\\\
&+(r_{a_{1}}+r_{b1})\boldsymbol{i}\\\ &+(r_{a_{2}}+r_{b_{2}})\boldsymbol{j}\\\
&+(r_{a_{3}}+r_{b_{3}})\boldsymbol{k}.\end{split}$ (3)

Scalar multiplication is defined as

$xQ=xr_{0}+xr_{1}\boldsymbol{i}+xr_{2}\boldsymbol{j}+xr_{3}\boldsymbol{k}.$
(4)

Figure 2: A high-level block diagram of the proposed QCNN SER system.

Element muultiplication is defined as:

$\begin{split}Q_{a}\times
Q_{b}&=(r_{a_{0}}r_{b_{0}}-r_{a_{1}}r_{b_{1}}-r_{a_{2}}r_{b_{2}}-r_{a_{3}}r_{b_{3}})\\\
&+(r_{a_{0}}r_{b_{1}}+r_{a_{1}}r_{b_{0}}+r_{a_{2}}r_{b_{3}}-r_{a_{3}}r_{b_{2}})\boldsymbol{i}\\\
&+(r_{a_{0}}r_{b_{2}}-r_{a_{1}}r_{b_{1}}+r_{a_{2}}r_{b_{0}}+r_{a_{3}}r_{b_{1}})\boldsymbol{j}\\\
&+(r_{a_{0}}r_{b_{3}}+r_{a_{1}}r_{b_{2}}-r_{a_{2}}r_{b_{1}}+r_{a_{3}}r_{b_{0}})\boldsymbol{k},\end{split}$
(5)

and rotation along axis $p$ is defined as

$\hat{R}=p\hat{Q}p^{-1}$ (6)

where

$p=\cos[\frac{\theta}{2}]+(p_{1}\times i+p_{2}\times j+p_{3}\times
k)\sin\frac{\theta}{2},$ (7) $\hat{Q}=r_{q_{0}}+r_{q_{1}}\times
i+r_{q_{2}}\times j+r_{q_{3}}\times k,$ (8)

####  3.2.2 Quaternion Convolution

We implement a quaternion color kernel similar to the one proposed in [13]. We
do this by defining the $50\times 75$ colored Mel-spectrogram image in a 3D
vector color space as a quaternion matrix $\hat{C}=[\hat{c}_{nn^{\prime}}]\in
H^{50\times 75}$. We can thus represent the color channels in $C$ as:

$\hat{C}=0+Ri+Gj+Bk$ (9)

where $R,G,$ and $B\in R^{50*75}$ are the red, green, and blue channels,
respectively. By representing a pixel in 3D vector color space, the proposed
model also deploys a quaternion kernel convolving around the input, in which
the element $\hat{w}_{ll^{\prime}}$ of the kernel $\hat{W}$ can be defined as:

$\hat{w}_{ll^{\prime}}=s_{ll^{\prime}}\Bigl{(}\cos\frac{\theta_{ll^{\prime}}}{2}+(\frac{\sqrt{3}}{3}(i+j+k))\sin\frac{\theta_{ll^{\prime}}}{2}\Bigr{)}$
(10)

where $\theta_{u^{\prime}}\in[-\pi,\pi]$ and $s_{u^{\prime}}\in R$. Thus the
quaternion convolution operation $*$ can be expressed as:

$\hat{C}\times\hat{W}=\hat{F}=[\hat{f}_{kk^{\prime}}]\in
H^{((50-L+1)\times(75-L+1))}$ (11)

where $\hat{W}$ is a $L\times L$ quaternion kernel and
$[\hat{f}_{kk^{\prime}}]$ is defined as:

$\hat{f}_{kk^{\prime}}=\sum_{l=1}\sum_{l^{\prime}=1}\frac{1}{s_{ll^{\prime}}}\hat{w}_{ll^{\prime}}\hat{c}_{(k+l)(k+l^{\prime})}\hat{w}_{ll^{\prime}}.$
(12)

By implementing quaternion convolution kernels to perform rotation and scaling
operations, the model finds better representative features in the decibel
color waveform space than real-valued convolution kernels. Additionally, the
real-valued kernels apply scaling and pixel transformation separately to the
three axes of color and create single-channel feature maps; whereas the
quaternion kernels can more intuitively capture the color space as a whole
without any type of color information loss since the channels are interrelated
as a 3D vector.

####  3.2.3 Connecting Layers

The quaternion convolution layer is connected to other typical layers to
construct the fully connected neural network. Specifically, max-pooling is
performed on the imaginary parts and ReLU is used to reset invalid quaternion
vector rotations to the nearest point in color space. In addition, the
imaginary parts of the output of a quaternion layer are represented as 3 real
numbers, such that a real-valued and vectorized output can be obtained to be
connected to real-valued fully-connected layer. A softmax layer is then
connected to this real-valued layer to train the QCNN model.

####  3.2.4 Model Architecture

Figure 2 shows a block diagram of the proposed model. The RGB-Mel spectrogram
is transformed into a 3D color vector space with rotation magnitudes for the
initial input layer. This layer is then followed by the quaternion
convolutional layer, which consists of 32 quaternion kernels (all quaternion
kernels, proposed in this paper, cover a 3x3 2D coordinate space) that perform
rotation and scaling operations. ReLU is used to reset invalid quaternion
vector rotations to the nearest point in color space. Afterward, max-pooling
is performed on the imaginary color properties. This is followed by 64
quaternion kernels. The layer is then flattened to vectorize the quaternion
layer to a real-valued layer. The real-valued dense layer is used to train the
model by computing the gradient of the categorical cross-entropy loss function
(which operates on softmax outputs for the classes of emotion).

###  3.3 Training

Weights are initialized as an imaginary quaternion with a uniform distribution
in the interval [0, 1]. The imaginary unit is then normalized, with well-known
criteria described by [19] to generate the quaternion weight $\hat{W}$ which
was already defined above. Backpropagation is used to update the quaternion
weights by applying the quaternion chain rule. Additionally, the loss function
can be computed on the real-valued layer; specifically, the categorical cross-
entropy function is used. Furthermore, the model used Adam as its stochastic
optimizer and overfitting was addressed with the use of dropout (with a
probability of 0.3). The model was trained for 50 epochs with 9 steps per
validation epoch. The full model was implemented in Python with Tensorflow
2.3.0 as the deep learning backend.

##  4 Experiments

In order to measure the accuracy and effectivity of the model, QCNN is tested
and compared against the most competitive models in SER on the RAVDESS,
IEMOCAP, and EMODB datasets. In addition, we also compared the size of the
models.

###  4.1 Speech Datasets

####  4.1.1 RAVDESS

The Ryerson Audio-Visual Database of Emotional Speech and Song (RAVDESS) [20]
consists of speech samples of 8 categorical emotions— calm, happy, sad, angry,
fearful, surprise, and disgust —spoken by 24 professional actors, vocalizing
lexically-matched statements in a neutral North American accent. The
classification accuracy of QCNN for the 8 classes of emotion on RAVDESS was
compared with results from CNN-LSTM [21], a pre-feature extracted SVM [22],
GResNets [23], a controlled human accuracy study [RAVDESS], and Deep-CNN [14].
The results are summarized in table I.

####  4.1.2 IEMOCAP

The Interactive Emotional Dyadic Motion Capture (IEMOCAP) [24] database is an
acted, multimodal and multispeaker database. It consists of 12 hours of
audiovisual data, including video, speech, motion capture of face, text
transcriptions. This paper makes use of the speech samples of 4 categorical
emotions— anger, happiness, sadness, and neutrality —spoken by 10 (5 female
and 5 male) professional actors in fluent English. The classification accuracy
of QCNN for the 4 classes of emotion on IEMOCAP was compared with results from
CNN-LSTM [21], DialogueRNN [25], BiF-AGRU [26], Deep-CNN [14], DialogueGCN
[27], BiERU-lc [28], and DSCNN [29]. The results are summarized in table II.

####  4.1.3 EMO-DB

The Berlin EMO-DB [30] consists of speech samples representing 7 categorical
emotions —anger, happiness, sadness, fear, disgust, boredom, and neutral
speech—spoken by 10 professional actors (5 female and 5 male) in German. Each
speaker acted 7 emotions for 10 different utterances (5 short and 5 long) with
emotionally neutral linguistic contents. The classification accuracy for the 7
classes of emotion on EMO-DB was compared for the Multi-task LSTM [31], RNN
[21], CNN-LSTM [21], and Deep-CNN [14] methods with QCNN. The results are
summarized in table III.

###  4.2 Model Size

The only model sizes that were reported for above models are Deep-CNN and CNN-
LSTM. However, both being competitive models, they can be used as a bench mark
to compare the QCNN model size. The results are summarized in Table IV.

Table 1: Classifier performance on RAVDESS _Model_ | _Accuracy (%)_
---|---
CNN-LSTM [21] | 53.08
Feature SVM | 60.10
GResNets [23] | 65.97
Human Accuracy [20] | 67.00
Deep-CNN [14] | 71.67
QCNN [ours] | 77.87

###  4.3 Discussion

The QCNN model presented yields state-of-the-art results on RAVDESS by 6̃%.
Although the model does not achieve state-of-the-art results on EMO-DB or
IEMOCAP, in both cases it underperforms the most competitive model by only
$<$2̃%. The balance of accuracy over model size is also considered in this
paper, and although the model slightly underperforms in the IEMPCAP and EMO-DB
datasets, the size of the model is significantly less. This is attributed to
the ability to store RGB arrays in a quaternion domain to reduce model size
[12]. However, in this paper, we did not explore encoding other features of
waveforms such as Mel-frequency cepstral coefficients (MFCC), chromograms,
spectral contrast features, or Tonnetz representations as pure quaternions. It
is possible that separate Quaternion models encoding these features,
separately or individually, could outperform the model presented in this
paper. Additionally, it is possible that an auxiliary neural network can be
used in conjunction with the QCNN to capture higher-level features. General
performance can also be improved by using effective waveform data augmentation
techniques. This research can also be extended to real-time SER for quaternion
blocks to compute and process on speech-waveforms. However, for such a model
it is expected to take advantage of other features rather than an RGB domain
of Mel-spectrograms because of delta-computational expense.

Table 2: Classifier performance on IEMOCAP _Model_ | _Unweighted Accuracy (%)_
---|---
CNN-LSTM [21] | 50.17
DialogueRNN [25] | 63.40
BiF-AGRU [26] | 63.50
Deep-CNN [14] | 64.30
DialogueGCN [27] | 65.25
BiERU-1c [28] | 66.11
QCNN [ours] | 70.46
DSCNN DSCNN[29] | 72.00

##  5 Conclusion

Current SER research is a dynamic, complex task relying on feature extraction
and computational classification. To capture this highly complex qualitative
information from speech waveforms, multiple feature extraction and
classification techniques have been introduced in literature. However, many
machine learning-based methods focus on higher-level features in the real-
valued space. This paper reports a unique approach to feature and network
encoding by using a quaternion structural model. Specifically, we encode the
RGB domain of Mel-spectrogram features in a quaternion input and use custom
quaternion convolutional layers to learn features in a quaternion space. These
layers are implemented in a standard neural network structure to train on
benchmark datasets such as RAVDESS, IEMOCAP, and EMO-DB. QCNN is reported to
yield an accuracy of 77.87%, 70.46%, and 88.78% for the RAVDESS, IEMOCAP, and
EMO-DB datasets, respectfully. In comparison to other competitive methods,
QCNN achieves state-of-the-art results on RAVDESS, and underperforms only one
method in both IEMOCAP and EMO-DB. Additionally, QCNN is able to exploit its
quaternion encoding for a reduced model size, confirming previous literature
on the topic as well as providing an opportunity for deployment on lightweight
machines such as voice-assistant devices. Given the significant results in
both accuracy and computational complexity, QCNN is expected to be a
foundational block for the continuation of SER research.

Table 3: Classifier performance on EMO-DB _Model_ | _Unweighted Accuracy (%)_
---|---
Multi-task LSTM [31] | 58.14
RNN[21] | 63.21
CNN-LSTM [21] | 69.72
QCNN [ours] | 88.78
Deep-CNN [14] | 90.01
Table 4: Model sizes (Mb) _Model_ | RAVDESS | EMO-DB | IEMOCAP
---|---|---|---
Deep CNN | 74.5 | 69.5 | 88.5
CNN-LSTM | 111.1 | 94.4 | 128.3
QCNN | 42 | 31.2 | 67.7

## References

  * [1] D. Rethage, J. Pons, and X. Serra, “A wavenet for speech denoising,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018, pp. 5069–5073.
  * [2] M. Benzeghiba, R. De Mori, O. Deroo, S. Dupont, T. Erbes, D. Jouvet, L. Fissore, P. Laface, A. Mertins, C. Ris _et al._ , “Automatic speech recognition and speech variability: A review,” _Speech communication_ , vol. 49, no. 10-11, pp. 763–786, 2007.
  * [3] S. Ramakrishnan and I. M. El Emary, “Speech emotion recognition approaches in human computer interaction,” _Telecommunication Systems_ , vol. 52, no. 3, pp. 1467–1478, 2013.
  * [4] A. Ashar, M. S. Bhatti, and U. Mushtaq, “Speaker identification using a hybrid cnn-mfcc approach,” in _2020 International Conference on Emerging Trends in Smart Technologies (ICETST)_. IEEE, 2020, pp. 1–4.
  * [5] S. An, Z. Ling, and L. Dai, “Emotional statistical parametric speech synthesis using lstm-rnns,” in _2017 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC)_. IEEE, 2017, pp. 1613–1616.
  * [6] S. Pascual, A. Bonafonte, and J. Serra, “Segan: Speech enhancement generative adversarial network,” _arXiv preprint arXiv:1703.09452_ , 2017.
  * [7] E. J. Humphrey, T. Cho, and J. P. Bello, “Learning a robust tonnetz-space transform for automatic chord recognition,” in _2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2012, pp. 453–456.
  * [8] D. Palaz, M. M. Doss, and R. Collobert, “Convolutional neural networks-based continuous speech recognition using raw speech signal,” in _2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2015, pp. 4295–4299.
  * [9] Y. Miao, M. Gowayyed, and F. Metze, “Eesen: End-to-end speech recognition using deep rnn models and wfst-based decoding,” in _2015 IEEE Workshop on Automatic Speech Recognition and Understanding (ASRU)_. IEEE, 2015, pp. 167–174.
  * [10] A. Graves, N. Jaitly, and A.-r. Mohamed, “Hybrid speech recognition with deep bidirectional lstm,” in _2013 IEEE workshop on automatic speech recognition and understanding_. IEEE, 2013, pp. 273–278.
  * [11] M. W”ollmer, B. Schuller, F. Eyben, and G. Rigoll, “Combining long short-term memory and dynamic bayesian networks for incremental emotion-sensitive artificial listening,” _IEEE Journal of selected topics in signal processing_ , vol. 4, no. 5, pp. 867–881, 2010.
  * [12] T. Isokawa, T. Kusakabe, N. Matsui, and F. Peper, “Quaternion neural network and its application,” in _International conference on knowledge-based and intelligent information and engineering systems_. Springer, 2003, pp. 318–324.
  * [13] X. Zhu, Y. Xu, H. Xu, and C. Chen, “Quaternion convolutional neural networks,” in _Proceedings of the European Conference on Computer Vision (ECCV)_ , 2018, pp. 631–647.
  * [14] D. Issa, M. F. Demirci, and A. Yazici, “Speech emotion recognition with deep convolutional neural networks,” _Biomedical Signal Processing and Control_ , vol. 59, p. 101894, 2020.
  * [15] W. Lim, D. Jang, and T. Lee, “Speech emotion recognition using convolutional and recurrent neural networks,” in _2016 Asia-Pacific signal and information processing association annual summit and conference (APSIPA)_. IEEE, 2016, pp. 1–4.
  * [16] A. Satt, S. Rozenberg, and R. Hoory, “Efficient emotion recognition from speech using deep learning on spectrograms.” in _Interspeech_ , 2017, pp. 1089–1093.
  * [17] T. Parcollet, Y. Zhang, M. Morchid, C. Trabelsi, G. Linar‘es, R. De Mori, and Y. Bengio, “Quaternion convolutional neural networks for end-to-end automatic speech recognition,” _arXiv preprint arXiv:1806.07789_ , 2018.
  * [18] E. O. Brigham and R. Morrow, “The fast fourier transform,” _IEEE spectrum_ , vol. 4, no. 12, pp. 63–70, 1967.
  * [19] X. Glorot and Y. Bengio, “Understanding the difficulty of training deep feedforward neural networks,” in _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , 2010, pp. 249–256.
  * [20] S. R. Livingstone and F. A. Russo, “The ryerson audio-visual database of emotional speech and song (ravdess): A dynamic, multimodal set of facial and vocal expressions in north american english,” _PloS one_ , vol. 13, no. 5, p. e0196391, 2018.
  * [21] J. Parry, D. Palaz, G. Clarke, P. Lecomte, R. Mead, M. Berger, and G. Hofer, “Analysis of deep learning architectures for cross-corpus speech emotion recognition.” in _INTERSPEECH_ , 2019, pp. 1656–1660.
  * [22] P. Shegokar and P. Sircar, “Continuous wavelet transform based speech emotion recognition,” in _2016 10th International Conference on Signal Processing and Communication Systems (ICSPCS)_. IEEE, 2016, pp. 1–8.
  * [23] Y. Zeng, H. Mao, D. Peng, and Z. Yi, “Spectrogram based multi-task audio classification,” _Multimedia Tools and Applications_ , vol. 78, no. 3, pp. 3705–3722, 2019.
  * [24] C. Busso, M. Bulut, C.-C. Lee, A. Kazemzadeh, E. Mower, S. Kim, J. N. Chang, S. Lee, and S. S. Narayanan, “Iemocap: Interactive emotional dyadic motion capture database,” _Language resources and evaluation_ , vol. 42, no. 4, p. 335, 2008.
  * [25] N. Majumder, S. Poria, D. Hazarika, R. Mihalcea, A. Gelbukh, and E. Cambria, “Dialoguernn: An attentive rnn for emotion detection in conversations,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, 2019, pp. 6818–6825.
  * [26] C. Li, J. Jiao, Y. Zhao, and Z. Zhao, “Combining gated convolutional networks and self-attention mechanism for speech emotion recognition,” in _2019 8th International Conference on Affective Computing and Intelligent Interaction Workshops and Demos (ACIIW)_. IEEE, 2019, pp. 105–109.
  * [27] D. Ghosal, N. Majumder, S. Poria, N. Chhaya, and A. Gelbukh, “Dialoguegcn: A graph convolutional neural network for emotion recognition in conversation,” _arXiv preprint arXiv:1908.11540_ , 2019.
  * [28] Y. Li, T. Zhao, and T. Kawahara, “Improved end-to-end speech emotion recognition using self attention mechanism and multitask learning.” in _Interspeech_ , 2019, pp. 2803–2807.
  * [29] S. Kwon _et al._ , “A cnn-assisted enhanced audio signal processing for speech emotion recognition,” _Sensors_ , vol. 20, no. 1, p. 183, 2020.
  * [30] F. Burkhardt, A. Paeschke, M. Rolfes, W. F. Sendlmeier, and B. Weiss, “A database of german emotional speech,” in _Ninth European Conference on Speech Communication and Technology_ , 2005.
  * [31] S. Goel and H. Beigi, “Cross lingual cross corpus speech emotion recognition,” _arXiv preprint arXiv:2003.07996_ , 2020.