# Transverse structure of the pion beyond leading twist with basis light-front
quantization

Zhimin Zhu<EMAIL_ADDRESS>Zhi Hu<EMAIL_ADDRESS>Jiangshan Lan
<EMAIL_ADDRESS>Chandan Mondal<EMAIL_ADDRESS>Xingbo Zhao
<EMAIL_ADDRESS>James P. Vary<EMAIL_ADDRESS>
(BLFQ Collaboration)  Institute of Modern Physics, Chinese Academy of
Sciences, Lanzhou, Gansu, 730000, China  School of Nuclear Physics, University
of Chinese Academy of Sciences, Beijing, 100049, China  CAS Key Laboratory of
High Precision Nuclear Spectroscopy, Institute of Modern Physics, Chinese
Academy of Sciences, Lanzhou 730000, China  Department of Physics and
Astronomy, Iowa State University, Ames, IA 50011, USA

###### Abstract

We investigate the twist-$3$ transverse-momentum-dependent parton distribution
functions (TMDs) of the pion with basis light-front quantization. The twist-3
TMDs are not independent and can be decomposed into twist-$2$ and genuine
twist-$3$ terms from the equations of motion (EOM). We compute the TMDs from
the resulting light-front wave functions obtained by diagonalizing the light-
front QCD Hamiltonian, determined for pion’s constituent quark-antiquark and
quark-antiquark-gluon Fock sectors, with three-dimensional confinement. We
also obtain the twist-3 parton distribution functions (PDFs) and show that
they preserve the sum rule, which affirms the robustness of our approach. This
is the first time that theoretical predictions are made for subleading twist
structures of the pion containing interference terms between two light-front
Fock sectors.

###### keywords:

Light-front quantization , Pions , Subleading twist TMDs and PDFs , Quark-
gluon correlations

##  1 Introduction

Understanding the structure of hadrons is one of the main goals of modern
particle and nuclear physics [1, 2, 3, 4, 5, 6]. However, due to our
insufficient knowledge of QCD color confinement and chiral symmetry breaking,
it is difficult to obtain information on the internal structure of hadrons
directly from the fundamental Lagrangian. On the other hand, with the
development of the factorization theorem [7, 8, 9, 6, 10, 11, 12, 13, 14], one
can decompose the cross sections of different experiments into short-distance
perturbative and large-distance nonperturbative parts. The nonperturbative
segment contains the structural information of the hadron [13, 15, 16, 17,
18]. For instance, one can extract the collinear parton distribution functions
(PDFs) from deep inelastic scattering (DIS) processes [19, 20]. The PDFs,
describing the distribution of longitudinal momentum carried by the
constituents inside hadrons [21, 7], provide us with a one-dimensional (1D)
longitudinal picture of hadrons. Therefore, the PDFs represent an initial 1D
stage for comparing theory with experiment in the very challenging
nonperturbative regime. This is our initial goal for addressing the structure
of the pion.

Beyond the 1D structure, we gain insights into the three-dimensional (3D)
structure of hadrons through the generalized parton distributions (GPDs) [22,
23, 24, 6, 25] and the transverse-momentum-dependent parton distribution
functions (TMDs) [26, 27, 28, 29]. The GPDs provide us with essential
information about the distribution in position space and the orbital motion of
partons inside hadrons [30, 31, 6, 32]. The GPDs are experimentally accessible
through the deeply virtual Compton scattering (DVCS) or the deeply virtual
meson production (DVMP) [25, 6, 33, 34, 35]. The TMDs provide information on
the distribution of the transverse momentum in addition to the longitudinal
momentum fraction [7]. The TMDs, especially of quarks, can be extracted from
the cross section associated with a specific azimuthal angle in the semi-
inclusive deep inelastic scattering (SIDIS) processes or the Drell-Yan
processes [7, 36, 37, 38]. These two classes of distribution functions allow
us to draw a 3D picture of hadrons.

Higher-twist parton distribution functions give access to a wealth of
information about the partonic structure of hadrons [21, 39, 40] such as
describing multiparton correlations inside the hadrons that correspond to the
interference between scattering from a single quark and from a coherent quark-
gluon pair [41, 42, 43, 44]. They help us to understand the quark-gluon
dynamics going beyond the probabilistic interpretation, which applies to the
parton distribution functions at the leading-twist. Twist-three PDFs and TMDs
contribute to various observables in inclusive and semi-inclusive DIS,
respectively. Due to the existence of suppression factors, the twist-three
contributions are often much smaller than the twist-two. However, in the
kinematics of fixed target experiments, the twist-three structure functions
are not small. One of the major goals of the future experimental program at
JLab$12$ is to determine the different higher-twist spin-azimuthal asymmetries
in SIDIS [45, 46, 47]. Meanwhile, the upcoming Electron-Ion-Colliders (EICs)
[48, 49, 50, 51] would access different kinematical regions [52, 1]. In this
context, twist-three TMDs have been investigated in various QCD inspired
models, for example, diquark spectator models [53, 54, 55, 56], the MIT bag
model [41, 57, 58, 59], chiral quark soliton models [60, 61, 62, 63, 64],
instanton models of QCD vacuum [65, 66], constituent parton model on the light
front [67, 68], and perturbative light-front Hamiltonian approaches with a
quark target [69, 70, 71, 72]. Recently, the unpolarized twist-2 and twist-3,
T-even quark TMDs in the pion have been evaluated by using the homogeneous
Bethe-Salpeter integral equation in Ref. [73]. In addition, phenomenological
extractions of the pion’s quark TMDs from available cross section data of
pion-nucleus induced Drell–Yan processes have been reported in Refs. [74, 75].

In this work, we present the twist-three TMDs of the pion obtained with basis
light-front quantization (BLFQ), which provides a nonperturbative framework to
solve relativistic many-body bound state problems in quantum field theories
[76, 77, 78, 79, 80, 81, 82, 83, 84, 32, 85, 86, 87, 88, 89, 90, 91, 92, 93,
94]. Previously, this approach has been successfully applied to investigate
the leading-twist TMDs of the physical spin-$\frac{1}{2}$ and spin-$1$
particles in QED, i.e., electron [89] and photon [85], respectively, as well
as for the proton [95]. We emphasize that we compute the pion TMDs at the
model scale and do not consider the TMD evolution [7, 96, 97, 98, 99]. Our
calculation of TMDs should be physically interpreted as its nonperturbative,
initial scale components. In other words, our model computation can provide a
useful starting point for the TMDs at a specific scale, but the TMD evolution
must be performed to obtain a complete distribution at a different scale.
However, this is not the focus of the present work.

With the framework of BLFQ [76], we adopt the light-front QCD Hamiltonian [86,
100] and solve for its mass eigenvalues and eigenstates. With quark ($q$),
antiquark ($\bar{q}$), and gluon ($g$) being the explicit degrees of freedom
for the strong interaction, our Hamiltonian includes light-front QCD
interactions relevant to constituent $|q\bar{q}\rangle$ and
$|q\bar{q}g\rangle$ Fock sectors of the mesons with a complementary 3D
confinement [79]. We solve this Hamiltonian in the leading two Fock components
and fix the model parameters by fitting the known mass spectra of unflavored
light mesons [86]. The resulting light-front wavefunctions (LFWFs) obtained as
eigenvectors of this Hamiltonian have been successfully employed to compute a
wide class of pion observables, e.g., the electromagnetic form factors and
associated charge radius, decay constant, quark and gluon PDFs, the
differential cross sections for charmonium production by a pion beam, etc.,
with remarkable overall success. Here, we extend our investigations of the
pion to compute its subleading-twist TMDs and predict, for the first time, the
genuine twist-3 TMDs, $\tilde{e}$, which preserves the sum rule, and
$\tilde{f}^{\perp}$ of the pion.

##  2 Meson LFWFs in the BLFQ framework

In light-front field theory, the solution of a bound state can be obtained by
solving the light-front stationary Schrödinger equation:
$P^{+}P^{-}|\Psi\rangle=M^{2}|\Psi\rangle$, where $P^{-}$ and $P^{+}$
represent the light-front Hamiltonian and the longitudinal momentum of the
system, respectively. At the fixed light-front time, $x^{+}\equiv
x^{0}+x^{3}$, the meson state with mass squared eigenvalue $M^{2}$ can be
expressed schematically in terms of various quark, antiquark, and gluon Fock
sectors,

$|\Psi\rangle=\psi_{(q\bar{q})}|q\bar{q}\rangle+\psi_{(q\bar{q}g)}|q\bar{q}g\rangle+\psi_{(q\bar{q}q\bar{q})}|q\bar{q}q\bar{q}\rangle+\cdots,$
(1)

where $\psi_{(\cdots)}$ is the LFWF associated with the Fock component
$|\cdots\rangle$. For numerical calculations, we must truncate the infinite
Fock sector expansion to a finite Fock space of Eq. (1), and here we retain
the first and the second Fock components. This implies that the meson, at the
model scale, is described by the quark-antiquark $\psi_{(q\bar{q})}$ and the
quark-antiquark-gluon $\psi_{(q\bar{q}g)}$ LFWFs.

We adopt an effective light-front Hamiltonian:
$P^{-}=P^{-}_{\text{QCD}}+P^{-}_{C}$ [86], where $P^{-}_{\text{QCD}}$ and
$P^{-}_{C}$ correspond to the light-front QCD Hamiltonian relevant to
constituent $|q\bar{q}\rangle$ and $|q\bar{q}g\rangle$ Fock components of the
mesons and a model Hamiltonian for the confining interaction, respectively.
The light-front QCD Hamiltonian with one dynamical gluon in the light-front
gauge $A^{+}=0$ is [86, 100]

$\displaystyle P_{\rm QCD}^{-}=$
$\displaystyle\int\mathrm{d}x^{-}\mathrm{d}^{2}x^{\perp}\Big{\\{}\frac{1}{2}\bar{\psi}\gamma^{+}\frac{m_{0}^{2}+(i\partial^{\perp})^{2}}{i\partial^{+}}\psi$
$\displaystyle+\frac{1}{2}A_{a}^{i}\left[m_{g}^{2}+(i\partial^{\perp})^{2}\right]A^{i}_{a}+g\bar{\psi}\gamma_{\mu}T^{a}A_{a}^{\mu}\psi$
$\displaystyle+\frac{1}{2}g^{2}\bar{\psi}\gamma^{+}T^{a}\psi\frac{1}{(i\partial^{+})^{2}}\bar{\psi}\gamma^{+}T^{a}\psi\Big{\\}},$
(2)

where $\psi$ and $A^{\mu}$ represent the quark and the gluon fields,
respectively. $T^{a}$ is the half Gell-Mann matrix, $T^{a}=\lambda^{a}/2$. $g$
is the coupling constant. $m_{0}$ and $m_{g}$ are the bare quark mass and the
model gluon mass, respectively. While the gluon mass is zero in QCD, we permit
a phenomenologically motivated gluon mass to fit the mass spectra in our low-
energy model [86]. For the quark mass correction from a higher Fock sector, we
introduce a mass counter term, $\delta m_{q}=m_{0}-m_{q}$, for the quark in
the leading Fock sector, where $m_{q}$ is the renormalized quark mass. In the
vertex interaction, we allow an independent quark mass $m_{f}$ referring to
Ref. [101].

The confinement in the leading Fock sector, which includes transverse and
longitudinal confining potentials, reads as [86, 79],

$\begin{split}&P_{\rm
C}^{-}P^{+}=\kappa^{4}\left\\{x(1-x)\vec{r}_{\perp}^{2}-\frac{\partial_{x}[x(1-x)\partial_{x}]}{(m_{q}+m_{\bar{q}})^{2}}\right\\},\end{split}$
(3)

where $\kappa$ is the strength of the confinement, and
$\vec{r}_{\perp}=\sqrt{x(1-x)}(\vec{r}_{\perp q}-\vec{r}_{\perp\bar{q}})$
represents the holographic variable [102]. Confinement in the
$|q\bar{q}g\rangle$ sector is implemented through the cutoff of the BLFQ basis
functions discussed below. However, since all mass eigenstates of the present
model will have a valence quark-antiquark component, all our solutions will be
confined due to the confining interaction in the quark-antiquark sector.

To compute the Hamiltonian matrix, we follow BLFQ [76] and consider the 2D
harmonic-oscillator (HO) basis functions, $\Phi_{nm}(\vec{p}_{\perp};b)$, with
scale parameter $b$ to explain the transverse degrees of freedom [79]. The
2D-HO wave function carries the radial and the angular quantum numbers denoted
by $n$ and $m$, respectively. We employ a plane-wave basis function to
describe the longitudinal motion of the Fock particle. The longitudinal motion
is confined to a 1D box of length $2L$ with antiperiodic (periodic) boundary
conditions for fermion (boson). The longitudinal momentum is then defined as
$p^{+}=\frac{2\pi}{L}k$, where $k$, an integer (half-integer) for boson
(fermion), is the longitudinal quantum number. We neglect the zero mode for
the boson. With the additional quantum number $\lambda$ for the light-front
helicity, each Fock-particle basis state involves four quantum numbers
$|\alpha_{i}\rangle=|k_{i},n_{i},m_{i},\lambda_{i}\rangle$ and the many-body
basis states are identified as the direct product of the Fock-particle basis
states $|\alpha\rangle=\otimes_{i}|\alpha_{i}\rangle$. In the case of Fock
sectors allowing multiple color-singlet states, we require an additional label
to identify each color-singlet state. Note that this is only necessitated for
the Fock sectors beyond the two that we retain here. In addition, our many-
body basis states have well-defined total angular momentum projection
$M_{J}=\sum_{i}(m_{i}+\lambda_{i})$.

For numerical calculations, we truncate the infinite basis of each Fock sector
by introducing two truncation parameters $K$ and $N_{\rm max}$ in longitudinal
and transverse directions, respectively. The dimensionless variable
$K=\sum_{i}k_{i}$ parameterizes the longitudinal momentum $P^{+}$. Thus, for a
Fock particle $i$, the longitudinal momentum fraction $x$ is defined as
$x_{i}=p^{+}/P^{+}=k_{i}/K$. The variable $K$ manifests as the ‘resolution’ in
the longitudinal direction, and hence a resolution onto the PDFs. The $N_{\rm
max}$ truncation is set by $\sum_{i}(2n_{i}+|m_{i}|+1)\leq N_{\rm max}$ and it
enables factorization of the transverse center of mass motion [78, 77]. The
$N_{\rm max}$ truncation implicitly acts as the infrared (IR) and ultraviolet
(UV) cutoffs. In momentum space, the IR cutoff $\lambda_{\rm IR}\simeq
b/\sqrt{N_{\rm max}}$ and the UV cutoff $\Lambda_{\rm UV}\simeq b\sqrt{N_{\rm
max}}$ [77].

After diagonalizing the full Hamiltonian matrix in the BLFQ framework, we
obtain the mass spectra $M^{2}$ and the corresponding LFWFs in momentum space,

$\displaystyle\Psi^{M_{J}}_{\mathcal{N},\\{\lambda_{i}\\}}({\\{x_{i},\vec{p}_{\perp
i}\\}})$
$\displaystyle=\sum_{\\{n_{i}m_{i}\\}}\psi^{M_{J}}_{\mathcal{N}}({\\{\alpha}_{i}\\})\prod_{i=1}^{\mathcal{N}}\Phi_{n_{i}m_{i}}(\vec{p}_{\perp
i},b)\,,$ (4)

where $\psi^{M_{J}}_{\mathcal{N}=2}(\\{\alpha_{i}\\})$ and
$\psi^{M_{J}}_{\mathcal{N}=3}(\\{\alpha_{i}\\})$ are the components of the
eigenvectors associated with the Fock sectors $|q\bar{q}\rangle$ and
$|q\bar{q}g\rangle$, respectively.

With the truncations $\\{N_{\text{max}},K\\}=\\{14,15\\}$, the model
parameters, summarized in Table 1, have been fixed to generate the mass
spectra of unflavored light mesons [86]. The LFWFs in this model provide a
high-quality description of the electromagnetic form factors and associated
charge radius, decay constant, PDFs for the pion, and pion-nucleus induced
$J/\psi$ production cross sections [86].

Table 1: The model parameters for truncation $\\{N_{\text{max}},K\\}=\\{14,15\\}$ [86]. All are in units of GeV except $g$. $m_{q}$ | $m_{g}$ | $b$ | $\kappa$ | ${m}_{f}$ | $g$
---|---|---|---|---|---
0.39 | 0.60 | 0.29 | 0.65 | 5.69 | 1.92

##  3 Transverse-momentum-dependent parton distributions

For spin-0 mesons, the quark TMDs are parameterized through the quark-quark
correlation function defined as [103],

$\displaystyle\Phi_{q}^{[\Gamma]}(x,k_{\perp})=$
$\displaystyle\frac{1}{2}\int\frac{\mathrm{d}z^{-}\mathrm{d}^{2}\vec{z}_{\perp}}{2(2\pi)^{3}}e^{ik\cdot
z}$ $\displaystyle\times\langle
P|\bar{\psi}(0)\Gamma\mathcal{W}(0,z)\psi(z)|P\rangle|_{z^{+}=0}$ (5)

with $k^{+}=xP^{+}$. $|P\rangle$ is the light-front bound state of the target
meson with mass $M$ and momenta $(P^{+},\vec{P}_{\perp})$, where the
transverse momentum $\vec{P}_{\perp}=\vec{0}$ [104]. The Dirac matrix $\Gamma$
determines the Lorentz structure of the correlator $\Phi^{[\Gamma]}$ and its
‘twist’ $\tau$ [42]. The Wilson line $\mathcal{W}$ preserves the gauge
invariance of the bilocal quark field operators in the correlation function
[105].

For the pion, there are two leading twist (twist-2) TMDs namely the
unpolarized quark TMD, $f^{q}_{1}(x,k_{\perp})$, and the polarized quark TMD,
$h_{1}^{\perp q}(x,k_{\perp})$, also known as the Boer-Mulders function, which
are defined through the parameterizations of the quark-quark correlator with
$\Gamma\equiv\gamma^{+},\,\sigma^{j+}\gamma_{5}$,

$\displaystyle\Phi_{q}^{[\gamma^{+}]}(x,k_{\perp})$
$\displaystyle=f^{q}_{1}(x,k_{\perp}),$ (6)
$\displaystyle\Phi_{q}^{[i\sigma^{j+}\gamma_{5}]}(x,k_{\perp})$
$\displaystyle=-\frac{\varepsilon_{T}^{ij}k_{\perp}^{i}}{M}{h_{1}^{\perp
q}}(x,k_{\perp}),$ (7)

where $k_{\perp}$ is the transverse momentum of a struck quark,
$\epsilon_{T}^{11}=\epsilon_{T}^{22}=0$, and
$\epsilon_{T}^{12}=-\epsilon^{21}=1$. Meanwhile, the subleading twist
(twist-3) quark TMDs are defined as follows

$\displaystyle\Phi^{[1]}_{q}(x,k_{\perp})$
$\displaystyle=\frac{M}{P^{+}}e^{q}(x,k_{\perp}),$ (8)
$\displaystyle\Phi_{q}^{[\gamma^{j}]}(x,k_{\perp})$
$\displaystyle=\frac{M}{P^{+}}\frac{k_{\perp}^{j}}{M}f^{\perp
q}(x,k_{\perp}),$ (9)
$\displaystyle\Phi_{q}^{[\gamma^{j}\gamma_{5}]}(x,k_{\perp})$
$\displaystyle=-\frac{M}{P^{+}}\frac{\varepsilon_{T}^{ij}k_{\perp}^{i}}{M}{g^{\perp
q}}(x,k_{\perp}),$ (10)
$\displaystyle\Phi_{q}^{[i\sigma^{ij}\gamma_{5}]}(x,k_{\perp})$
$\displaystyle=-\frac{M}{P^{+}}\varepsilon_{T}^{ij}{h}^{q}(x,k_{\perp}).$ (11)

The TMDs $h_{1}^{\perp q}(x,k_{\perp})$, $g^{\perp q}(x,k_{\perp})$ and
$h^{q}(x,k_{\perp})$ are T-odd TMDs in nature. The rest are T-even TMDs [103].
Note that the twist-3 TMD correlation functions have a suppression with the
factor $M/P^{+}$. In this present work, we do not consider the effect of the
Wilson line $\mathcal{W}$ and set it to be a unit matrix. With this
approximation: $\mathcal{W}\approx\mathbbm{1}$, only the T-even TMDs survive
[67].

In the light-front field theory [106], the quark field is decomposed into a
‘good’ $\psi_{+}$, and a ‘bad’ $\psi_{-}$ component
($\psi_{\pm}=\frac{1}{4}\gamma^{\mp}\gamma^{\pm}\psi$)111The projection
operators used here are different from the definition in Ref. [106]. We follow
the convention, $x^{\pm}\equiv x^{0}\pm x^{3}$.. The bad component $\psi_{-}$
can be expressed in terms of the good component $\psi_{+}$ and the gauge
fields as,

$\psi_{-}(z)=\frac{\gamma^{+}}{2i\partial^{+}}\left[i(\partial_{j}-igA_{\perp
j}(z))\gamma_{j}+m_{q}\right]\psi_{+}(z),$ (12)

where $m_{q}$ is the quark mass, and $j=1,2$. The above constraint equation
comes from the EOM in the light-cone gauge $A^{+}=0$. Note that in Eqs. (6)
and (7), the Dirac matrices in the leading-twist operators project the
correlator into $\bar{\psi}_{+}\psi_{+}$, while in Eqs. (8)-(11), the Dirac
matrices in the subleading-twist operators project the correlator into
$\bar{\psi}_{+}\psi_{-}+\bar{\psi}_{-}\psi_{+}$. Therefore, the twist-3 TMDs
can not be expressed as the probabilities or the differences of probabilities.
They are related to quark-quark matrix elements and quark-quark-gluon matrix
elements.

The higher-twist operators can be decomposed into the twist-2 and other
operators using the EOM of quark fields [43]. Consequently, the twist-2 (Eqs.
(6) and (7)) and the twist-3 (Eqs. (8)-(11)) TMDs are related via the
following relations [67]

$\displaystyle xe^{q}(x,k_{\perp})$
$\displaystyle=x\tilde{e}^{q}(x,k_{\perp})+\frac{m_{q}}{M}f^{q}_{1}(x,k_{\perp}),$
(13) $\displaystyle xf^{\perp q}(x,k_{\perp})$
$\displaystyle=x\tilde{f}^{\perp q}(x,k_{\perp})+f^{q}_{1}(x,k_{\perp}),$ (14)

where the tilde terms or the so-called genuine twist-3 TMDs are defined as
follows (see appendix for detailed derivation)

$\displaystyle\tilde{e}^{q}(x,k_{\perp})=\frac{i}{Mx}\frac{g}{2}\int\frac{\mathrm{d}z^{-}\mathrm{d}^{2}\vec{z}_{\perp}}{2(2\pi)^{3}}e^{ik\cdot
z}$ $\displaystyle\times\langle
P|\bar{\psi}_{+}(0)\sigma^{j+}[A_{\perp}^{j}(z)-A^{j}_{\perp}(0)]\psi_{+}(z)|P\rangle|_{z^{+}=0},$
(15) $\displaystyle 2xk^{j}_{\perp}\tilde{f}^{\perp
q}(x,k_{\perp})=g\int\frac{\mathrm{d}z^{-}\mathrm{d}^{2}\vec{z}_{\perp}}{2(2\pi)^{3}}e^{ik\cdot
z}$ $\displaystyle\times\langle
P|\bar{\psi}_{+}(0)\gamma^{+}[A_{\perp}^{j}(0)+A_{\perp}^{j}(z)]\psi_{+}(z)|P\rangle|_{z^{+}=0}.$
(16)

Sum rules are particularly important for testing the model consistency. The
twist-2 and the twist-3 integrated TMDs,
$j(x)=\int\mathrm{d}^{2}{k}_{\perp}j(x,k_{\perp})$, where $j(x,k_{\perp})$
denotes the TMDs, must satisfy the following sum rules [67]

$\displaystyle\int\mathrm{d}xf^{q}_{1}(x)$ $\displaystyle=N_{q},$ (17)
$\displaystyle\sum_{q}\int\mathrm{d}xxf_{1}^{q}(x)$
$\displaystyle=1-\int\mathrm{d}xxf_{1}^{g}(x),$ (18)
$\displaystyle\sum_{q}\int\mathrm{d}xe^{q}(x)$
$\displaystyle=\frac{\sigma_{\pi}}{m_{q}},$ (19)
$\displaystyle\int\mathrm{d}xxe^{q}(x)$ $\displaystyle=\frac{m_{q}}{M}N_{q}.$
(20)

Equation (17) represents the number sum rule, where $N_{q}$ is the total
number of valence quarks of a particular species [107]. Equation (18)
corresponds to the momentum sum rule, which states that the constituent quark,
antiquark, and gluon together carry the entire light front momentum of the
meson. Note that due to the explicit $k_{\perp}^{j}$ factor in Eq. (9), there
does not exist any PDF nor sum rules corresponding to the TMD $f^{\perp
q}(x,k_{\perp})$. Since the focus of the present work is on the quark
distributions, we omit the superscript $q$ from the TMDs below for simplicity.

The zeroth moment of $e(x)$, Eq. (19), involves the singular term at $x=0$ and
is related to the scalar form factor of the pion at zero momentum transfer,
$\sigma_{\pi}$ [41, 43]. The second sum rule of $e(x)$, Eq. (20), can be cast
into the first-moment sum rule of $\tilde{e}(x)$ using the relation of TMDs at
different twists in Eq. (13) and the number sum rule given in Eq. (17),

$\int\mathrm{d}x\,x\,\tilde{e}(x)=0.$ (21)

We will discuss this sum rule within our model in the following section.

### Overlap representations of TMDs

Here, we provide an overview of the LFWF overlap representations for the
leading-twist and the subleading-twist TMDs of the pion. In the Fock space,
the meson state in Eq. (3) is expressed by considering all Fock components as

$\displaystyle|P\rangle$
$\displaystyle=\sum_{n}\sum_{\lambda_{1},\lambda_{2},\cdots,\lambda_{n}}\int\prod^{n}_{1}\frac{\left[\mathrm{d}x_{i}\mathrm{d}^{2}\vec{k}_{\perp
i}\right]}{2(2\pi)^{3}\sqrt{x_{i}}}$ $\displaystyle\times
2(2\pi)^{3}\delta\left(1-\sum_{1}^{n}x_{i}\right)\delta^{2}\left(\vec{P}_{\perp}-\sum_{1}^{n}\vec{k}_{\perp
i}\right)$
$\displaystyle\times\Psi_{n,\\{\lambda_{i}\\}}\left(\left\\{x_{i},k_{\perp
i}\right\\}\right)\left|\\{x_{i}P^{+},\vec{k}_{\perp
i}+x_{i}\vec{P}_{\perp},\lambda_{i}\\}\right\rangle,$ (22)

with the following orthonormal relations,

$\langle
P^{\prime}|P\rangle=\sum_{n,n^{\prime}}2P^{+}(2\pi)^{3}\delta^{3}(P-P^{\prime})\delta_{nn^{\prime}}P_{n},$
(23)

where
$\delta^{3}(P-P^{\prime})=\delta(P^{+}-P^{\prime+})\delta^{2}(\vec{P}_{\perp}-\vec{P}_{\perp}^{\prime})$.
The summation is over all the Fock components and $P_{n}$ defines the
probability of finding the $n$-parton state in the meson. In this work, we
have only two LFWFs, $\psi_{2}$ for $|q\bar{q}\rangle$ and $\psi_{3}$ for
$|q\bar{q}g\rangle$. With the model parameters summarized in Table 1, we find
that the probabilities of the $|q\bar{q}\rangle$ and $|q\bar{q}g\rangle$
components are $49.2\%$ and $50.8\%$, respectively.

Using the expression of the meson state, Eq. (22), and the quantized quark and
gluon fields in the Lepage-Brodsky prescription [100, 108], we express the
TMDs $f_{1}(x,k_{\perp})$, $\tilde{e}(x,k_{\perp})$, and
$\tilde{f}^{\perp}(x,k_{\perp})$ in terms of the overlap of the LFWFs as

$\displaystyle f_{1}(x,k_{\perp}$
$\displaystyle)=\sum_{\lambda_{1},\lambda_{2}}\mathrm{d}[12]\Psi^{*}_{2,\lambda_{1}\lambda_{2}}\Psi_{2,\lambda_{1}\lambda_{2}}\delta^{3}(\tilde{p}^{\prime}_{1}-\tilde{k})$
$\displaystyle+\sum_{\lambda_{1},\lambda_{2},\lambda_{g}}$
$\displaystyle\int\mathrm{d}[123]\Psi^{*}_{3,\lambda_{1}\lambda_{2}\lambda_{g}}\Psi_{3,\lambda_{1}\lambda_{2}\lambda_{g}}\delta^{3}(\tilde{p}_{1}-\tilde{k}),$
(24) $\displaystyle\tilde{e}(x,k_{\perp}$
$\displaystyle)=-\frac{gC_{F}\sqrt{2}}{Mx}\sum_{\lambda_{2},\lambda_{g}}\int\mathrm{d}[12]\mathrm{d}[123]\frac{2(2\pi)^{3}}{\sqrt{x_{g}}}$
$\displaystyle\times(\Psi^{*}_{2,-\lambda_{2}}\Psi_{3,+\lambda_{2}\lambda_{g}}-\Psi_{2,+\lambda_{2}}\Psi^{*}_{3,-\lambda_{2}\lambda_{g}})$
$\displaystyle\times[\delta^{3}(\tilde{p}^{\prime}_{1}-\tilde{k})-\delta^{3}(\tilde{p}_{1}-\tilde{k})]\delta^{3}(\tilde{p}^{\prime}_{2}-\tilde{p}_{2}),$
(25) $\displaystyle\tilde{f}^{\perp}(x,k_{\perp}$
$\displaystyle)=\frac{gC_{F}\sqrt{2}}{2xk^{1}_{\perp}}\sum_{\lambda_{1},\lambda_{2},\lambda_{g}}\int\mathrm{d}[12]\mathrm{d}[123]\frac{2(2\pi)^{3}}{\sqrt{x_{g}}}(-)^{\frac{\lambda_{g}+1}{2}}$
$\displaystyle\times[\Psi^{*}_{2,\lambda_{1}\lambda_{2}}\Psi_{3,\lambda_{1}\lambda_{2}\lambda_{g}}\delta^{3}(\tilde{p}_{1}-\tilde{k})$
$\displaystyle+\Psi_{2,\lambda_{1}\lambda_{2}}\Psi^{*}_{3,\lambda_{1}\lambda_{2},\lambda_{g}}\delta^{3}(\tilde{p}^{\prime}_{1}-\tilde{k})]\delta^{3}(\tilde{p}_{2}^{\prime}-\tilde{p}_{2}).$
(26)

Here we define the conventions

$\displaystyle\mathrm{d}[12]\equiv\frac{\prod_{1}^{2}\mathrm{d}x^{\prime}_{i}\mathrm{d}^{2}\vec{p}^{\;\prime}_{\perp
i}}{[2(2\pi)^{3}]^{2}}2(2\pi)^{3}\delta^{3}(\tilde{P}-\sum_{1}^{2}\tilde{p}^{\prime}_{i}),$
(27)
$\displaystyle\mathrm{d}[123]\equiv\frac{\prod_{1}^{3}\mathrm{d}x_{i}\mathrm{d}^{2}\vec{p}_{\perp
i}}{[2(2\pi)^{3}]^{3}}2(2\pi)^{3}\delta^{3}(\tilde{P}-\sum_{1}^{3}\tilde{p}_{i}),$
(28)

where
$\delta^{3}(\tilde{P}-\sum_{1}^{n}\tilde{p}_{i})=\delta(1-\sum_{1}^{n}x_{i})\delta^{2}(\vec{P}_{\perp}-\sum_{1}^{n}\vec{p}_{\perp
i})$, and we omit the arguments
$(\tilde{p}^{\prime}_{1},\tilde{p}^{\prime}_{2})$ from the quark-antiquark
LFWFs and $(\tilde{p}_{1},\tilde{p}_{2},\tilde{p}_{g})$ from the quark-
antiquark-gluon LFWFs, where $\tilde{k}\equiv(x,\vec{k}_{\perp})$. The color
factor, $C_{F}={2}/{\sqrt{3}}$, comes from the quark-quark-gluon matrix
element $\langle
q\bar{q}|\bar{\psi}_{b}A^{\mu,a}T^{a}_{bc}\psi_{c}|q\bar{q}g\rangle$ and its
hermitian conjugate. The unpolarized TMD $f_{1}(x,k_{\perp})$ involves the
individual overlaps of LFWFs of the $|q\bar{q}\rangle$ and the
$|q\bar{q}g\rangle$ Fock sectors. However, the genuine twist-3 TMDs
$\tilde{e}(x,k_{\perp})$ and $\tilde{f}^{\perp}(x,k_{\perp})$ involve the
mixed overlaps of LFWFs between the $|q\bar{q}\rangle$ and the
$|q\bar{q}g\rangle$ Fock sectors. This is interpreted as an interference
between scattering from a coherent quark-gluon pair and from a single quark
[42, 41, 43].

##  4 Numerical results

With the model parameters for the basis truncation
$\\{N_{\text{max}},K\\}=\\{14,15\\}$ summarized in Table 1, we solve the
light-front stationary Schrödinger equation to obtain the two Fock-state LFWFs
of the pion with the mass $M=0.139~{}\mathrm{GeV}$. The LFWFs
$\Psi_{n,\\{\lambda_{i}\\}}(\\{x_{i},\vec{k}_{\perp i}\\})$ are boost-
invariant and depend only on the longitudinal fraction $x_{i}$ and the
relative transverse momentum of the constituent parton of the meson
$\vec{k}_{\perp i}$. The physical transverse momentum of the parton is given
by $\vec{p}_{\perp i}=\vec{k}_{\perp i}+x_{i}\vec{P}_{\perp}$ and the
longitudinal momentum is $p_{i}^{+}=k_{i}^{+}=x_{i}P^{+}$. We convert the
nonperturbative solutions in the single particle coordinates, Eq. (4), to
those in the relative coordinates,
$\Psi_{n,\\{\lambda_{i}\\}}(\\{x_{i},\vec{k}_{\perp i}\\})$, by factorizing
out the center of mass motion [78, 84, 89]. We then employ the resulting LFWFs
with the expressions in Eqs. (24-26) to compute the quark TMDs and PDFs of the
pion.

###  4.1 Twist-2 and twist-3 TMDs

The TMDs computed in BLFQ have oscillations in the transverse momentum
direction due to the oscillatory nature of the HO basis functions applied in
the transverse plane. The transverse truncation $N_{\rm max}$, together with
the HO scale $b$ determines the transverse basis states. The average over the
BLFQ results at different $N_{\rm max}$ decreases these finite basis
artifacts. Following the averaging procedure adopted in Refs. [89, 85], we
take an average of the BLFQ results for three different $N_{\mathrm{max}}$ at
fixed $K$. We first take the average between the BLFQ results at
$N_{\mathrm{max}}=n$ and $N_{\mathrm{max}}=n+2$. We consider another average
of the results obtained at $N_{\mathrm{max}}=n+2$ and $N_{\mathrm{max}}=n+4$.
The final result is then achieved by averaging the previous two averages.
Thus, this two-step averaging method involves the BLFQ results at
$N_{\mathrm{max}}=\\{n,\,n+2,\,n+4\\}$, and at fixed $K$.

Figure 1 shows the pion’s TMDs with different $N_{\mathrm{max}}$ and $K=15$.
The upper panel of Fig. 1 presents the TMDs as functions of $k_{\perp}$ at
fixed longitudinal momentum fraction $x=2.5/15$, whereas the lower panel of
Fig. 1 shows the TMDs as functions of $x$ at fixed $k_{\perp}=0.22$ GeV. We
observe that the BLFQ results at different $N_{\mathrm{max}}$ have
oscillations. The magnitude of the oscillation is larger when $k_{\perp}$ and
$x$ are small. The phases of the oscillations alter with changing
$N_{\mathrm{max}}$. Through the above method, we average the BLFQ results of
different $N_{\mathrm{max}}$, which reduces the oscillation of our results.
This can be seen in Fig. 1 (solid black lines). Meanwhile, the amplitude of
the oscillations is small in the large $k_{\perp}$ and $x>0.1$ region
reflecting better stability of our results in that domain.

The 3D structure of our BLFQ results for the TMDs $f_{1}(x,k_{\perp})$,
$\tilde{e}(x,k_{\perp})$, and $\tilde{f}^{\perp}(x,k_{\perp})$, after
implementing the average, is shown in Fig. 2. We find that the TMDs still have
small oscillations in the $k_{\perp}$ direction after averaging. The TMDs drop
to zero abruptly in the transverse plane after the UV cutoff ($\sim
b\sqrt{N_{\text{max}}}$). Note that the UV cutoff increases with increasing
$N_{\text{max}}$ and thus, the fall-off of the TMDs extends to the relevant
higher-momentum region for larger $N_{\text{max}}$. This behavior of TMDs
within our BLFQ framework has been observed for the physical
spin-$\frac{1}{2}$ particle in QED, i.e., electrons reported in Ref. [89].

Figure 1: The BLFQ results for the twist-2 TMD, $f_{1}(x,k_{\perp})$, and the
genuine twist-3 TMDs, $\tilde{e}(x,k_{\perp})$ and
$\tilde{f}^{\perp}(x,k_{\perp})$ at different truncation parameters. Upper
panel: the TMDs are functions of $k_{\perp}$ for fixed $x=2.5/15$. Lower
panel: the TMDs are functions of $x$ for fixed $k_{\perp}=0.22$ GeV. The solid
black lines represent the TMDs after implementing the averaging scheme with
truncation parameters $N_{\text{max}}=\\{12,14,16\\}$, $K=15$.

Figure 2: 3D plots of the BLFQ results for the twist-2 TMD,
$f_{1}(x,k_{\perp})$, and the genuine twist-3 TMDs, $\tilde{e}(x,k_{\perp})$
and $\tilde{f}^{\perp}(x,k_{\perp})$. Results are all obtained by averaging
over the BLFQ results at $N_{\text{max}}=\\{12,14,16\\}$, $K=15$.

Figure 3 presents the total twist-3 TMDs, $e(x,k_{\perp})$ and
$f^{\perp}(x,k_{\perp})$, using the EOM relations given in Eqs. (13) and (14).
We find that the qualitative behavior of these TMDs is similar to each other
and their magnitudes are larger than that of the twist-2 TMD
$f_{1}(x,k_{\perp})$. The behavior of these TMDs mainly follows from the
function $f_{1}(x,k_{\perp})/x$, which is singular at $x=0$. However, the
contributions of twist-3 TMDs in the cross section are suppressed by a factor
$M/P^{+}$ [109] and thus, the higher twist effects are weakened in high-energy
collisions. It can also be noticed that the twist-3 TMDs fall quickly in the
transverse direction.

Figure 3: 3D plots of the BLFQ results for the twist-3 TMDs, $e(x,k_{\perp})$
and $f^{\perp}(x,k_{\perp})$. Results are all obtained by averaging over the
BLFQ results at $N_{\text{max}}=\\{12,14,16\\}$, $K=15$.

###  4.2 $k_{\perp}$ moments of the TMDs

We consider the interpretation of the quark unpolarized TMD
$f_{1}(x,k_{\perp})$, which describes the probability of having an unpolarized
active quark with the longitudinal momentum fraction $x$ and the relative
transverse momentum $k_{\perp}$ in an unpolarized hadron [21]. We then define
the $x$-dependent mean squared transverse momentum of the quark as

$\langle
k_{\perp}^{2}\rangle_{f_{1}}(x)=\frac{\int\mathrm{d}^{2}k_{\perp}k_{\perp}^{2}f_{1}(x,k_{\perp})}{\int\mathrm{d}^{2}k_{\perp}f_{1}(x,k_{\perp})}.$
(29)

Figure 4 shows the results for the $x$-dependence of the mean squared
transverse momentum $\langle k_{\perp}^{2}\rangle_{f_{1}}(x)$ of the quark
with different Fock components: $|q\bar{q}\rangle$ and
$|q\bar{q}\rangle+|q\bar{q}g\rangle$. Both the distributions have their peaks
at $x=0.5$. We observe that the $\langle k_{\perp}^{2}\rangle$ distribution
from the leading Fock sector is symmetric about $x=0.5$. This occurs due to
the same quark and antiquark masses in the pion. Meanwhile, the $\langle
k_{\perp}^{2}\rangle$ distribution from $|q\bar{q}\rangle+|q\bar{q}g\rangle$
is smaller than that from $|q\bar{q}\rangle$. This is attributed to the fact
that the dynamical gluon in $|q\bar{q}g\rangle$ carries away a part of the
quark (antiquark) momentum.

Although there is no probabilistic interpretation for twist-3 TMDs $e$ and
$f^{\perp}$, we calculate their first and second $k_{\perp}$ moment defined as
[59]

$\langle k_{\perp}^{n}\rangle_{\rm
TMD}=\frac{\int\mathrm{d}x\int\mathrm{d}^{2}{k}_{\perp}k_{\perp}^{n}{\rm
TMD}}{\int\mathrm{d}x\int\mathrm{d}^{2}{k}_{\perp}\mathrm{TMD}}.$ (30)

The corresponding results are summarized in Table 2. Our results are close to
those from the light-front constituent model (LFCM) at a low scale
$\mu_{0}\sim 0.5$ GeV [67]. Note that the latter model’s results are based on
the leading Fock component.

Figure 4: The $x$-dependence of the mean squared transverse momentum $k_{\perp}$ of the unpolarized TMD $f_{1}$ for the leading Fock sector $|q\bar{q}\rangle$ and two Fock sectors $|q\bar{q}\rangle+|q\bar{q}g\rangle$ after averaging with the truncation parameters $N_{\text{max}}=\\{12,14,16\\}$, $K=15$ under the BLFQ framework. Table 2: $\langle k_{\perp}\rangle$ and $\langle k_{\perp}^{2}\rangle$ moments as defined in Eq. (30) for the twist-2 TMD, $f_{1}(x,k_{\perp})$, and the twist-3 TMDs, $e(x,k_{\perp})$ and $f^{\perp}(x,k_{\perp})$. The BLFQ results are all obtained by averaging with the truncation parameters $N_{\text{max}}=\\{12,14,16\\}$, $K=15$. The LFCM results are from Ref. [67]. All are in units of GeV.  | $\langle k_{\perp}\rangle_{\mathrm{BLFQ}}$ | $\langle k^{2}_{\perp}\rangle^{1/2}_{\mathrm{BLFQ}}$ | $\langle k_{\perp}\rangle_{\mathrm{LFCM}}$ | $\langle k^{2}_{\perp}\rangle^{1/2}_{\mathrm{LFCM}}$
---|---|---|---|---
$f_{1}$ | 0.26 | 0.30 | 0.28 | 0.32
$e$ | 0.26 | 0.30 | 0.26 | 0.30
$f^{\perp}$ | 0.25 | 0.29 | 0.26 | 0.30

###  4.3 Twist-2 and twist-3 PDFs

In this subsection, we turn our attention to the integrated TMDs namely the
PDFs. Figure 5 shows the twist-2 PDF $xf_{1}(x)$, and the genuine twist-3 PDFs
$x\tilde{e}(x)$ and $x\tilde{f}^{\perp}(x)$. Interestingly, and in contrast to
our TMD results, the PDFs have no oscillations and are very stable for
different $N_{\text{max}}$ truncations. Clearly, integrating the TMDs over
$k_{\perp}$ eliminates the unphysical oscillations and yields smooth results
for the PDFs. The twist-2 PDF $xf_{1}(x)$ shows a positive distribution, while
the genuine twist-3 functions, $x\tilde{e}(x)$ and $x\tilde{f}^{\perp}(x)$,
exhibit positive behavior when $x<0.5$, but they are negative above $x=0.5$.
The magnitude of the twist-2 distribution is much higher than those of the
twist-3 distributions. We also find that the integrated value of the
distribution $x\tilde{e}(x)$ over $x$ vanishes showing that our BLFQ results
properly satisfy the sum rule defined in Eq. (21). Note that this sum rule is
derived from the EOM relations and thus, it must be independent of the system.
The genuine twist-3 PDF $\tilde{e}(x)$ for the nucleon has been studied in
Ref. [68], where the authors considered $|qqq\rangle$ and $|qqqg\rangle$ Fock
components but their distribution $\tilde{e}(x)$ does not satisfy the sum rule
defined in Eq. (21). In Ref. [71], the author investigated the twist-3 PDF
$e(x)$ and the genuine twist-3 PDF $\tilde{e}(x)$ for a dressed quark at one
loop, where the distribution $\tilde{e}(x)$ satisfies the first-moment sum
rule (21). Using the sum rule defined in Eq. (19), we compute the pion scalar
form factor from our BLFQ results, $\sigma_{\pi}=13.44$ GeV.

Figure 5: The results for the twist-2 unpolarized PDF and the genuine twist-3
PDFs with the different $N_{\text{max}}$ under the BLFQ framework.

Figure 6 shows the total twist-3 PDFs and their decompositions in terms of
twist-2 and genuine twist-3 PDFs. We observe that the twist-2 terms and the
genuine twist-3 PDFs have distinctly different $x$-dependence behaviors. The
former strongly dominates the total twist-3 PDFs.

Figure 6: The results for the total twist-3 PDFs, the genuine twist-3 PDFs,
and the contributions from the twist-2 PDFs after averaging with the
truncation parameters $N_{\text{max}}=\\{12,14,16\\}$, $K=15$ under the BLFQ
framework.

Figure 7: Comparison of the unpolarized PDFs $f_{1}(x)$, $e(x)$ and
$f^{\perp}(x)$ of the LFCM [67] and the results after normalization to the
leading Fork state $|q\bar{q}\rangle$ and after averaging
$N_{\text{max}}=\\{12,14,16\\}$ of the BLFQ framework.

Figure 7 compares our integrated TMDs $f_{1}(x)$, $e(x)$, and $f^{\perp}(x)$
computed from only contributions of the leading Fock sector $|q\bar{q}\rangle$
with the LFCM’s results [67], which demonstrates the pion structure based on a
two-body light-front Fock state $|q\bar{q}\rangle$. For $f_{1}(x)$, both
models produce symmetric distributions about $x=0.5$ due to equal quark and
antiquark mass. However, the $f_{1}(x)$ in the LFCM is wider, while our BLFQ
result is narrower and steeper. The PDF $e(x)$ in the LFCM rises rapidly from
zero in the small $x$ region and decreases quickly in the large $x$ region.
Meanwhile, the BLFQ result shows a different trend for $e(x)$ dropping to a
plateau in the small $x$ region and then continuing to drop to zero. We notice
that $f^{\perp}(x)$ exhibits a somewhat similar trend as observed for $e(x)$.
At large-$x$, both the LFCM and the BLFQ framework provide qualitatively
similar behavior for $f^{\perp}(x)$.

##  5 Summary

In this work, we calculated the twist-2 TMD $f_{1}(x,k_{\perp})$ and the
genuine twist-3 TMDs $\tilde{e}(x,k_{\perp})$ and
$\tilde{f}^{\perp}(x,k_{\perp})$ of the pion at the model scale from its LFWFs
within the basis light-front quantization framework. These wave functions have
been obtained from the eigenvectors of the light-front QCD Hamiltonian in the
light-cone gauge for the light mesons by considering them within
$|q\bar{q}\rangle$ and $|q\bar{q}g\rangle$ Fock spaces, together with a 3D
confinement in the leading Fock sector. We employed the model-independent EOM
relations to obtain the total twist-3 TMDs $e(x,k_{\perp})$ and
$f^{\perp}(x,k_{\perp})$. In this study, the gauge link has been set to unity,
which leaves us only the T-even quark TMDs. We also adopted an averaging
procedure for reducing finite basis artifacts of our BLFQ results for the
TMDs. The integrated TMDs exhibit good consistency between computations with
different basis truncations. Our BLFQ results satisfy the momentum sum rule
related to $e(x)$. For the total twist-3 PDFs ${e}(x)$ and ${f}^{\perp}(x)$,
we found that the contributions from the genuine twist-3 PDFs have a
distinctive $x$-dependence from those of the twist-2 PDF. Due to imprecise
knowledge of the nonperturbative Sudakov-like form factors and complexity of
the evolution of high-twist TMDs [110] at present, it is not feasible yet to
evolve our twist-3 TMDs in order to employ them for predicting experimental
observables, such as spin asymmetries [109].

We compared our results from the leading Fock sector $|q\bar{q}\rangle$ with
the LFCM. Our twist-3 PDFs show a different trend at small $x$ from the
prediction of the LFCM, while at large $x$ both approaches provide similar
qualitative behaviors. Considering that the genuine twist-3 TMDs and the T-odd
TMDs with gauge links have similar operator structures, the results of this
work are encouraging for an extension to the T-odd TMDs. Future studies will
focus on the Boer-Mulders TMD $h_{1}^{\perp}$ of pions and its application to
spin-asymmetries in Drell-Yan processes. In addition, a detailed analysis of
the pion’s gluon TMDs in the BLFQ will be reported in a future study.

We believe that these results including twist-3 TMDs provide predictions for
future experiments, such as precise measurements of the pion-nucleus induced
Drell-Yan processes and future measurements of the 3D structure of the pion at
the EIC and at the EicC [111], which will refine our 3D knowledge of the pion
and its internal correlations based on QCD.

## Acknowledgements

We thank Siqi Xu, Jiatong Wu, Tiancai Peng and Zhe Liu for many helpful
discussions. J. L. is supported by Special Research Assistant Funding Project,
Chinese Academy of Sciences. C. M. is supported by new faculty start up
funding the Institute of Modern Physics, Chinese Academy of Sciences, Grants
No. E129952YR0. C. M. also thanks the Chinese Academy of Sciences Presidents
International Fellowship Initiative for the support via Grants No. 2021PM0023.
X. Z. is supported by new faculty startup funding by the Institute of Modern
Physics, Chinese Academy of Sciences, by Key Research Program of Frontier
Sciences, Chinese Academy of Sciences, Grant No. ZDB-SLY-7020, by the Natural
Science Foundation of Gansu Province, China, Grant No. 20JR10RA067, by the
Foundation for Key Talents of Gansu Province, by the Central Funds Guiding the
Local Science and Technology Development of Gansu Province, Grant No.
22ZY1QA006, and by the Strategic Priority Research Program of the Chinese
Academy of Sciences, Grant No. XDB34000000. J. P. V. is supported by the
Department of Energy under Grants No. DE-FG02-87ER40371, No. DE-SC0018223
(SciDAC4/NUCLEI) and DE-SC0023495 (SciDAC5/NUCLEI). A portion of the
computational resources were also provided by Taiyuan Advanced Computing
Center.

##  Appendix A Derivation of the EOM relations

In this appendix, we show the details of the derivation of the EOM relations
given in Eqs. (13) and (14). In Ref. [59], the authors derived EOM relations
in a free quark model without considering the gauge link. Here, we also ignore
the gauge links, however, our bound state contains a dynamical gluon.

For an arbitrary Dirac matrix $\Gamma$, we insert the EOM in the unintegrated
correlation function, and we obtain

$\displaystyle 0$ $\displaystyle=\int\frac{d^{4}z}{(2\pi)^{4}}e^{ik\cdot
z}\langle P|\bar{\psi}(0)\Gamma(i\cancel{D}(z)-m_{q})\psi(z)|P\rangle$
$\displaystyle=\int\frac{d^{4}z}{(2\pi)^{4}}e^{ik\cdot z}\langle
P|\bar{\psi}(0)\Gamma[\cancel{k}+g\cancel{A}(z)-m_{q}]\psi(z)|P\rangle,$ (31)

where the first term is derived from the surface integral. Adding this formula
and its complex conjugate yields

$\displaystyle 0=\int\frac{d^{4}z}{(2\pi)^{4}}e^{ik\cdot z}\langle P|$
$\displaystyle\bar{\psi}(0)\\{[\cancel{k}+g\cancel{A}(0)-m_{q}]\bar{\Gamma}$
$\displaystyle+$
$\displaystyle\Gamma[\cancel{k}+g\cancel{A}(z)-m_{q}]\\}\psi(z)|P\rangle,$
(32)

where $\bar{\Gamma}=\gamma^{0}\Gamma\gamma^{0}$. When $\Gamma=\gamma^{+}$, we
integrate the unintegrated correlator of Eq. (32) over the light-front energy
component $k^{-}$ and the light-front time component $z^{+}$ in the light-cone
gauge $A^{+}=0$, and we obtain

$\displaystyle 0=\int$ $\displaystyle\frac{d^{3}z}{(2\pi)^{4}}e^{ik\cdot
z}\langle P|\bar{\psi}(0)[2k^{+}-2m_{q}\gamma^{+}$ $\displaystyle-
ig(A^{j}(z)-A^{j}(0))\sigma^{j+}]\psi(z)|P\rangle|_{z^{+}=0}.$ (33)

From the above equation, we can find three kinds of gamma matrix structures,
$\mathbbm{1}$, $\gamma^{+}$ and $\sigma^{j+}$. The three gamma structures, by
definitions in Eqs. (6) and (8), correspond to the twist-3 TMD,
$e(x,k_{\perp})$, the twist-2 TMD, $f_{1}(x,k_{\perp})$, and the genuine
twist-3 TMD, $\tilde{e}(x,k_{\perp})$, respectively. Arranging the formula, we
get the EOM relation,

$\displaystyle xe(x,k_{\perp})$
$\displaystyle=x\tilde{e}(x,k_{\perp})+\frac{m_{q}}{M}f_{1}(x,k_{\perp}),$
(34)

where the genuine twist-3 TMD is defined as follows

$\displaystyle\tilde{e}(x,k_{\perp})=\frac{i}{Mx}\frac{g}{2}\int\frac{\mathrm{d}z^{-}\mathrm{d}^{2}z_{\perp}}{2(2\pi)^{3}}e^{ikz}$
$\displaystyle\times\langle
P|\bar{\psi}_{+}(0)\sigma^{j+}[A_{\perp}^{j}(z)-A^{j}_{\perp}(0)]\psi_{+}(z)|P\rangle|_{z^{+}=0}.$
(35)

We can also get a similar equation by setting $\Gamma=-i\sigma^{j+}$ in Eq.
(32) for the another T-even twist-3 TMD $\tilde{f}^{\perp}$ of spin-0 hadrons.

The above derivation does not take into account the contributions from gauge
links, but future inclusion of the gauge link is feasible yet will leave the
EOM augmented with additional singular terms [43, 59, 68].

## References

  * [1] A. Accardi, et al., Electron Ion Collider: The Next QCD Frontier: Understanding the glue that binds us all, Eur. Phys. J. A 52 (9) (2016) 268\.  arXiv:1212.1701, doi:10.1140/epja/i2016-16268-9.
  * [2] A. Bacchetta, M. Diehl, K. Goeke, A. Metz, P. J. Mulders, M. Schlegel, Semi-inclusive deep inelastic scattering at small transverse momentum, JHEP 02 (2007) 093.  arXiv:hep-ph/0611265, doi:10.1088/1126-6708/2007/02/093.
  * [3] A. V. Belitsky, X. Ji, F. Yuan, Final state interactions and gauge invariant parton distributions, Nucl. Phys. B 656 (2003) 165–198.  arXiv:hep-ph/0208038, doi:10.1016/S0550-3213(03)00121-4.
  * [4] H.-L. Lai, M. Guzzi, J. Huston, Z. Li, P. M. Nadolsky, J. Pumplin, C. P. Yuan, New parton distributions for collider physics, Phys. Rev. D 82 (2010) 074024\.  arXiv:1007.2241, doi:10.1103/PhysRevD.82.074024.
  * [5] J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. M. Nadolsky, W. K. Tung, New generation of parton distributions with uncertainties from global QCD analysis, JHEP 07 (2002) 012.  arXiv:hep-ph/0201195, doi:10.1088/1126-6708/2002/07/012.
  * [6] M. Diehl, Generalized parton distributions, Phys. Rept. 388 (2003) 41–277.  arXiv:hep-ph/0307382, doi:10.1016/j.physrep.2003.08.002.
  * [7] J. Collins, Foundations of perturbative QCD, Vol. 32, Cambridge University Press, 2013.
  * [8] M. Diehl, D. Ostermeier, A. Schafer, Elements of a theory for multiparton interactions in QCD, JHEP 03 (2012) 089, [Erratum: JHEP 03, 001 (2016)].  arXiv:1111.0910, doi:10.1007/JHEP03(2012)089.
  * [9] J. C. Collins, L. Frankfurt, M. Strikman, Factorization for hard exclusive electroproduction of mesons in QCD, Phys. Rev. D 56 (1997) 2982–3006.  arXiv:hep-ph/9611433, doi:10.1103/PhysRevD.56.2982.
  * [10] T. C. Rogers, An overview of transverse-momentum–dependent factorization and evolution, Eur. Phys. J. A 52 (6) (2016) 153.  arXiv:1509.04766, doi:10.1140/epja/i2016-16153-7.
  * [11] G. F. Sterman, Partons, factorization and resummation, TASI 95, in: Theoretical Advanced Study Institute in Elementary Particle Physics (TASI 95): QCD and Beyond, 1995, pp. 327–408.  arXiv:hep-ph/9606312.
  * [12] J. C. Collins, D. E. Soper, G. F. Sterman, Factorization for Short Distance Hadron - Hadron Scattering, Nucl. Phys. B 261 (1985) 104–142.  doi:10.1016/0550-3213(85)90565-6.
  * [13] J. C. Collins, D. E. Soper, The Theorems of Perturbative QCD, Ann. Rev. Nucl. Part. Sci. 37 (1987) 383–409.  doi:10.1146/annurev.ns.37.120187.002123.
  * [14] J. C. Collins, Hard scattering factorization with heavy quarks: A General treatment, Phys. Rev. D 58 (1998) 094002.  arXiv:hep-ph/9806259, doi:10.1103/PhysRevD.58.094002.
  * [15] V. L. Chernyak, A. R. Zhitnitsky, Asymptotic Behavior of Exclusive Processes in QCD, Phys. Rept. 112 (1984) 173.  doi:10.1016/0370-1573(84)90126-1.
  * [16] D. Müller, D. Robaschik, B. Geyer, F. M. Dittes, J. Hořejši, Wave functions, evolution equations and evolution kernels from light ray operators of QCD, Fortsch. Phys. 42 (1994) 101–141.  arXiv:hep-ph/9812448, doi:10.1002/prop.2190420202.
  * [17] B. Lampe, E. Reya, Spin physics and polarized structure functions, Phys. Rept. 332 (2000) 1–163.  arXiv:hep-ph/9810270, doi:10.1016/S0370-1573(99)00100-3.
  * [18] E. V. Shuryak, Quantum Chromodynamics and the Theory of Superdense Matter, Phys. Rept. 61 (1980) 71–158.  doi:10.1016/0370-1573(80)90105-2.
  * [19] O. Hen, G. A. Miller, E. Piasetzky, L. B. Weinstein, Nucleon-Nucleon Correlations, Short-lived Excitations, and the Quarks Within, Rev. Mod. Phys. 89 (4) (2017) 045002.  arXiv:1611.09748, doi:10.1103/RevModPhys.89.045002.
  * [20] R. D. Ball, et al., Parton distributions from high-precision collider data, Eur. Phys. J. C 77 (10) (2017) 663.  arXiv:1706.00428, doi:10.1140/epjc/s10052-017-5199-5.
  * [21] R. L. Jaffe, Parton Distribution Functions for Twist Four, Nucl. Phys. B 229 (1983) 205–230.  doi:10.1016/0550-3213(83)90361-9.
  * [22] X.-D. Ji, Off forward parton distributions, J. Phys. G 24 (1998) 1181–1205.  arXiv:hep-ph/9807358, doi:10.1088/0954-3899/24/7/002.
  * [23] M. V. Polyakov, C. Weiss, Skewed and double distributions in pion and nucleon, Phys. Rev. D 60 (1999) 114017.  arXiv:hep-ph/9902451, doi:10.1103/PhysRevD.60.114017.
  * [24] K. Goeke, M. V. Polyakov, M. Vanderhaeghen, Hard exclusive reactions and the structure of hadrons, Prog. Part. Nucl. Phys. 47 (2001) 401–515.  arXiv:hep-ph/0106012, doi:10.1016/S0146-6410(01)00158-2.
  * [25] A. V. Belitsky, A. V. Radyushkin, Unraveling hadron structure with generalized parton distributions, Phys. Rept. 418 (2005) 1–387.  arXiv:hep-ph/0504030, doi:10.1016/j.physrep.2005.06.002.
  * [26] J. Collins, New definition of TMD parton densities, Int. J. Mod. Phys. Conf. Ser. 4 (2011) 85–96.  arXiv:1107.4123, doi:10.1142/S2010194511001590.
  * [27] R. Angeles-Martinez, et al., Transverse Momentum Dependent (TMD) parton distribution functions: status and prospects, Acta Phys. Polon. B 46 (12) (2015) 2501–2534.  arXiv:1507.05267, doi:10.5506/APhysPolB.46.2501.
  * [28] M. Diehl, Introduction to GPDs and TMDs, Eur. Phys. J. A 52 (6) (2016) 149.  arXiv:1512.01328, doi:10.1140/epja/i2016-16149-3.
  * [29] A. Bacchetta, Where do we stand with a 3-D picture of the proton?, Eur. Phys. J. A 52 (6) (2016) 163.  arXiv:2107.06772, doi:10.1140/epja/i2016-16163-5.
  * [30] M. Burkardt, Impact parameter dependent parton distributions and off forward parton distributions for zeta —$>$ 0, Phys. Rev. D 62 (2000) 071503, [Erratum: Phys.Rev.D 66, 119903 (2002)].  arXiv:hep-ph/0005108, doi:10.1103/PhysRevD.62.071503.
  * [31] M. Burkardt, Impact parameter space interpretation for generalized parton distributions, Int. J. Mod. Phys. A 18 (2003) 173–208.  arXiv:hep-ph/0207047, doi:10.1142/S0217751X03012370.
  * [32] Y. Liu, S. Xu, C. Mondal, X. Zhao, J. P. Vary, Angular momentum and generalized parton distributions for the proton with basis light-front quantization, Phys. Rev. D 105 (9) (2022) 094018.  arXiv:2202.00985, doi:10.1103/PhysRevD.105.094018.
  * [33] X.-D. Ji, Deeply virtual Compton scattering, Phys. Rev. D 55 (1997) 7114–7125.  arXiv:hep-ph/9609381, doi:10.1103/PhysRevD.55.7114.
  * [34] L. Favart, M. Guidal, T. Horn, P. Kroll, Deeply Virtual Meson Production on the nucleon, Eur. Phys. J. A 52 (6) (2016) 158.  arXiv:1511.04535, doi:10.1140/epja/i2016-16158-2.
  * [35] S. Boffi, B. Pasquini, Generalized parton distributions and the structure of the nucleon, Riv. Nuovo Cim. 30 (9) (2007) 387–448.  arXiv:0711.2625, doi:10.1393/ncr/i2007-10025-7.
  * [36] D. Boer, P. J. Mulders, Time reversal odd distribution functions in leptoproduction, Phys. Rev. D 57 (1998) 5780–5786.  arXiv:hep-ph/9711485, doi:10.1103/PhysRevD.57.5780.
  * [37] S. J. Brodsky, D. S. Hwang, I. Schmidt, Initial state interactions and single spin asymmetries in Drell-Yan processes, Nucl. Phys. B 642 (2002) 344–356.  arXiv:hep-ph/0206259, doi:10.1016/S0550-3213(02)00617-X.
  * [38] S. J. Brodsky, D. S. Hwang, I. Schmidt, Final state interactions and single spin asymmetries in semiinclusive deep inelastic scattering, Phys. Lett. B 530 (2002) 99–107.  arXiv:hep-ph/0201296, doi:10.1016/S0370-2693(02)01320-5.
  * [39] P. J. Mulders, R. D. Tangerman, The Complete tree level result up to order 1/Q for polarized deep inelastic leptoproduction, Nucl. Phys. B 461 (1996) 197–237, [Erratum: Nucl.Phys.B 484, 538–540 (1997)].  arXiv:hep-ph/9510301, doi:10.1016/0550-3213(95)00632-X.
  * [40] K. Goeke, A. Metz, M. Schlegel, Parameterization of the quark-quark correlator of a spin-1/2 hadron, Phys. Lett. B 618 (2005) 90–96.  arXiv:hep-ph/0504130, doi:10.1016/j.physletb.2005.05.037.
  * [41] R. L. Jaffe, X.-D. Ji, Chiral odd parton distributions and Drell-Yan processes, Nucl. Phys. B 375 (1992) 527–560.  doi:10.1016/0550-3213(92)90110-W.
  * [42] R. L. Jaffe, X.-D. Ji, Chiral odd parton distributions and polarized Drell-Yan, Phys. Rev. Lett. 67 (1991) 552–555.  doi:10.1103/PhysRevLett.67.552.
  * [43] A. V. Efremov, P. Schweitzer, The Chirally odd twist 3 distribution e(a)(x), JHEP 08 (2003) 006.  arXiv:hep-ph/0212044, doi:10.1088/1126-6708/2003/08/006.
  * [44] M. Burkardt, Transverse force on quarks in deep-inelastic scattering, Phys. Rev. D 88 (2013) 114502.  arXiv:0810.3589, doi:10.1103/PhysRevD.88.114502.
  * [45] H. Avakian, et al., Jlab experiment (2008).
  * [46] H. Avakian, et al., Jlab experiment (2006).
  * [47] S. Pisano, et al., Jlab experiment (2014).
  * [48] R. Abdul Khalek, et al., Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report, Nucl. Phys. A 1026 (2022) 122447\.  arXiv:2103.05419, doi:10.1016/j.nuclphysa.2022.122447.
  * [49] R. Abdul Khalek, et al., Snowmass 2021 White Paper: Electron Ion Collider for High Energy Physics (3 2022).  arXiv:2203.13199.
  * [50] S. Amoroso, et al., Snowmass 2021 whitepaper: Proton structure at the precision frontier (3 2022).  arXiv:2203.13923.
  * [51] D. P. Anderle, et al., Electron-ion collider in China, Front. Phys. (Beijing) 16 (6) (2021) 64701.  arXiv:2102.09222, doi:10.1007/s11467-021-1062-0.
  * [52] D. Boer, et al., Gluons and the quark sea at high energies: Distributions, polarization, tomography (8 2011).  arXiv:1108.1713.
  * [53] R. Jakob, P. J. Mulders, J. Rodrigues, Modeling quark distribution and fragmentation functions, Nucl. Phys. A 626 (1997) 937–965.  arXiv:hep-ph/9704335, doi:10.1016/S0375-9474(97)00588-5.
  * [54] Z. Lu, I. Schmidt, T-odd quark-gluon-quark correlation function in the diquark model, Phys. Lett. B 712 (2012) 451–455.  arXiv:1202.0700, doi:10.1016/j.physletb.2012.05.023.
  * [55] W. Mao, Z. Lu, Beam spin asymmetries of charged and neutral pion production in semi-inclusive DIS, Eur. Phys. J. C 73 (2013) 2557.  arXiv:1306.1004, doi:10.1140/epjc/s10052-013-2557-9.
  * [56] W. Mao, Z. Lu, B.-Q. Ma, Transverse single-spin asymmetries of pion production in semi-inclusive DIS at subleading twist, Phys. Rev. D 90 (1) (2014) 014048\.  arXiv:1405.3876, doi:10.1103/PhysRevD.90.014048.
  * [57] A. I. Signal, Calculations of higher twist distribution functions in the MIT bag model, Nucl. Phys. B 497 (1997) 415–434.  arXiv:hep-ph/9610480, doi:10.1016/S0550-3213(97)00231-9.
  * [58] H. Avakian, A. V. Efremov, P. Schweitzer, F. Yuan, The transverse momentum dependent distribution functions in the bag model, Phys. Rev. D 81 (2010) 074035\.  arXiv:1001.5467, doi:10.1103/PhysRevD.81.074035.
  * [59] C. Lorcé, B. Pasquini, P. Schweitzer, Unpolarized transverse momentum dependent parton distribution functions beyond leading twist in quark models, JHEP 01 (2015) 103.  arXiv:1411.2550, doi:10.1007/JHEP01(2015)103.
  * [60] P. Schweitzer, The Chirally odd twist three distribution function e**alpha(x) in the chiral quark soliton model, Phys. Rev. D 67 (2003) 114010.  arXiv:hep-ph/0303011, doi:10.1103/PhysRevD.67.114010.
  * [61] M. Wakamatsu, Comparative analysis of the transversities and the longitudinally polarized distribution functions of the nucleon, Phys. Lett. B 653 (2007) 398–403.  arXiv:0705.2917, doi:10.1016/j.physletb.2007.08.013.
  * [62] M. Wakamatsu, Y. Ohnishi, The Nonperturbative origin of delta function singularity in the chirally odd twist three distribution function e(x), Phys. Rev. D 67 (2003) 114011.  arXiv:hep-ph/0303007, doi:10.1103/PhysRevD.67.114011.
  * [63] Y. Ohnishi, M. Wakamatsu, pi N sigma term and chiral odd twist three distribution function e(x) of the nucleon in the chiral quark soliton model, Phys. Rev. D 69 (2004) 114002.  arXiv:hep-ph/0312044, doi:10.1103/PhysRevD.69.114002.
  * [64] C. Cebulla, J. Ossmann, P. Schweitzer, D. Urbano, The Twist-3 parton distribution function e**a(x) in large-N(c) chiral theory, Acta Phys. Polon. B 39 (2008) 609–640.  arXiv:0710.3103.
  * [65] J. Balla, M. V. Polyakov, C. Weiss, Nucleon matrix elements of higher twist operators from the instanton vacuum, Nucl. Phys. B 510 (1998) 327–364.  arXiv:hep-ph/9707515, doi:10.1016/S0550-3213(98)00638-5.
  * [66] B. Dressler, M. V. Polyakov, On the twist - three contribution to h(L) in the instanton vacuum, Phys. Rev. D 61 (2000) 097501.  arXiv:hep-ph/9912376, doi:10.1103/PhysRevD.61.097501.
  * [67] C. Lorcé, B. Pasquini, P. Schweitzer, Transverse pion structure beyond leading twist in constituent models, Eur. Phys. J. C 76 (7) (2016) 415.  arXiv:1605.00815, doi:10.1140/epjc/s10052-016-4257-8.
  * [68] B. Pasquini, S. Rodini, The twist-three distribution $e^{q}(x,k_{\perp})$ in a light-front model, Phys. Lett. B 788 (2019) 414–424.  arXiv:1806.10932, doi:10.1016/j.physletb.2018.11.033.
  * [69] M. Burkardt, Y. Koike, Violation of sum rules for twist three parton distributions in QCD, Nucl. Phys. B 632 (2002) 311–329.  arXiv:hep-ph/0111343, doi:10.1016/S0550-3213(02)00263-8.
  * [70] R. Kundu, A. Metz, Higher twist and transverse momentum dependent parton distributions: A Light front Hamiltonian approach, Phys. Rev. D 65 (2002) 014009\.  arXiv:hep-ph/0107073, doi:10.1103/PhysRevD.65.014009.
  * [71] A. Mukherjee, Twist Three Distribution e(x): Sum Rules and Equation of Motion Relations, Phys. Lett. B 687 (2010) 180–183.  arXiv:0912.1446, doi:10.1016/j.physletb.2010.03.023.
  * [72] A. Accardi, A. Bacchetta, W. Melnitchouk, M. Schlegel, What can break the Wandzura-Wilczek relation?, JHEP 11 (2009) 093.  arXiv:0907.2942, doi:10.1088/1126-6708/2009/11/093.
  * [73] E. Ydrefors, W. de Paula, T. Frederico, G. Salmè, Unpolarized transverse-momentum dependent distribution functions of a quark in a pion with Minkowskian dynamics (1 2023).  arXiv:2301.11599.
  * [74] M. Cerutti, L. Rossi, S. Venturini, A. Bacchetta, V. Bertone, C. Bissolotti, M. Radici, Extraction of pion transverse momentum distributions from Drell-Yan data, Phys. Rev. D 107 (1) (2023) 014014.  arXiv:2210.01733, doi:10.1103/PhysRevD.107.014014.
  * [75] P. C. Barry, L. Gamberg, W. Melnitchouk, E. Moffat, D. Pitonyak, A. Prokudin, N. Sato, Tomography of pions and protons via transverse momentum dependent distributions (2 2023).  arXiv:2302.01192.
  * [76] J. P. Vary, H. Honkanen, J. Li, P. Maris, S. J. Brodsky, A. Harindranath, G. F. de Teramond, P. Sternberg, E. G. Ng, C. Yang, Hamiltonian light-front field theory in a basis function approach, Phys. Rev. C 81 (2010) 035205.  arXiv:0905.1411, doi:10.1103/PhysRevC.81.035205.
  * [77] X. Zhao, H. Honkanen, P. Maris, J. P. Vary, S. J. Brodsky, Electron g-2 in Light-Front Quantization, Phys. Lett. B 737 (2014) 65–69.  arXiv:1402.4195, doi:10.1016/j.physletb.2014.08.020.
  * [78] P. Wiecki, Y. Li, X. Zhao, P. Maris, J. P. Vary, Basis Light-Front Quantization Approach to Positronium, Phys. Rev. D 91 (10) (2015) 105009.  arXiv:1404.6234, doi:10.1103/PhysRevD.91.105009.
  * [79] Y. Li, P. Maris, X. Zhao, J. P. Vary, Heavy Quarkonium in a Holographic Basis, Phys. Lett. B 758 (2016) 118–124.  arXiv:1509.07212, doi:10.1016/j.physletb.2016.04.065.
  * [80] S. Jia, J. P. Vary, Basis light front quantization for the charged light mesons with color singlet Nambu–Jona-Lasinio interactions, Phys. Rev. C 99 (3) (2019) 035206.  arXiv:1811.08512, doi:10.1103/PhysRevC.99.035206.
  * [81] W. Qian, S. Jia, Y. Li, J. P. Vary, Light mesons within the basis light-front quantization framework, Phys. Rev. C 102 (5) (2020) 055207.  arXiv:2005.13806, doi:10.1103/PhysRevC.102.055207.
  * [82] S. Tang, Y. Li, P. Maris, J. P. Vary, Heavy-light mesons on the light front, Eur. Phys. J. C 80 (6) (2020) 522.  arXiv:1912.02088, doi:10.1140/epjc/s10052-020-8081-9.
  * [83] C. Mondal, S. Xu, J. Lan, X. Zhao, Y. Li, D. Chakrabarti, J. P. Vary, Proton structure from a light-front Hamiltonian, Phys. Rev. D 102 (1) (2020) 016008\.  arXiv:1911.10913, doi:10.1103/PhysRevD.102.016008.
  * [84] S. Xu, C. Mondal, J. Lan, X. Zhao, Y. Li, J. P. Vary, Nucleon structure from basis light-front quantization, Phys. Rev. D 104 (9) (2021) 094036.  arXiv:2108.03909, doi:10.1103/PhysRevD.104.094036.
  * [85] S. Nair, C. Mondal, X. Zhao, A. Mukherjee, J. P. Vary, Basis light-front quantization approach to photon, Phys. Lett. B 827 (2022) 137005.  arXiv:2201.12770, doi:10.1016/j.physletb.2022.137005.
  * [86] J. Lan, K. Fu, C. Mondal, X. Zhao, j. P. Vary, Light mesons with one dynamical gluon on the light front, Phys. Lett. B 825 (2022) 136890.  arXiv:2106.04954, doi:10.1016/j.physletb.2022.136890.
  * [87] L. Adhikari, C. Mondal, S. Nair, S. Xu, S. Jia, X. Zhao, J. P. Vary, Generalized parton distributions and spin structures of light mesons from a light-front Hamiltonian approach, Phys. Rev. D 104 (11) (2021) 114019.  arXiv:2110.05048, doi:10.1103/PhysRevD.104.114019.
  * [88] C. Mondal, S. Nair, S. Jia, X. Zhao, J. P. Vary, Pion to photon transition form factors with basis light-front quantization, Phys. Rev. D 104 (9) (2021) 094034.  arXiv:2109.02279, doi:10.1103/PhysRevD.104.094034.
  * [89] Z. Hu, S. Xu, C. Mondal, X. Zhao, J. P. Vary, Transverse structure of electron in momentum space in basis light-front quantization, Phys. Rev. D 103 (3) (2021) 036005.  arXiv:2010.12498, doi:10.1103/PhysRevD.103.036005.
  * [90] J. Lan, C. Mondal, M. Li, Y. Li, S. Tang, X. Zhao, J. P. Vary, Parton Distribution Functions of Heavy Mesons on the Light Front, Phys. Rev. D 102 (1) (2020) 014020.  arXiv:1911.11676, doi:10.1103/PhysRevD.102.014020.
  * [91] J. Lan, C. Mondal, S. Jia, X. Zhao, J. P. Vary, Pion and kaon parton distribution functions from basis light front quantization and QCD evolution, Phys. Rev. D 101 (3) (2020) 034024.  arXiv:1907.01509, doi:10.1103/PhysRevD.101.034024.
  * [92] J. Lan, C. Mondal, S. Jia, X. Zhao, J. P. Vary, Parton Distribution Functions from a Light Front Hamiltonian and QCD Evolution for Light Mesons, Phys. Rev. Lett. 122 (17) (2019) 172001.  arXiv:1901.11430, doi:10.1103/PhysRevLett.122.172001.
  * [93] S. Xu, C. Mondal, X. Zhao, Y. Li, J. P. Vary, Nucleon spin decomposition with one dynamical gluon (9 2022).  arXiv:2209.08584.
  * [94] T. Peng, Z. Zhu, S. Xu, X. Liu, C. Mondal, X. Zhao, J. P. Vary, Basis light-front quantization approach to $\Lambda$ and $\Lambda$c and their isospin triplet baryons, Phys. Rev. D 106 (11) (2022) 114040.  arXiv:2208.00355, doi:10.1103/PhysRevD.106.114040.
  * [95] Z. Hu, S. Xu, C. Mondal, X. Zhao, J. P. Vary, Transverse momentum structure of proton within the basis light-front quantization framework, Phys. Lett. B 833 (2022) 137360.  arXiv:2205.04714, doi:10.1016/j.physletb.2022.137360.
  * [96] J. C. Collins, D. E. Soper, Back-To-Back Jets: Fourier Transform from B to K-Transverse, Nucl. Phys. B 197 (1982) 446–476.  doi:10.1016/0550-3213(82)90453-9.
  * [97] S. Catani, D. de Florian, M. Grazzini, Universality of nonleading logarithmic contributions in transverse momentum distributions, Nucl. Phys. B 596 (2001) 299–312.  arXiv:hep-ph/0008184, doi:10.1016/S0550-3213(00)00617-9.
  * [98] G. Bozzi, S. Catani, D. de Florian, M. Grazzini, Transverse-momentum resummation and the spectrum of the Higgs boson at the LHC, Nucl. Phys. B 737 (2006) 73–120.  arXiv:hep-ph/0508068, doi:10.1016/j.nuclphysb.2005.12.022.
  * [99] G. Bozzi, S. Catani, G. Ferrera, D. de Florian, M. Grazzini, Transverse-momentum resummation: A Perturbative study of Z production at the Tevatron, Nucl. Phys. B 815 (2009) 174–197.  arXiv:0812.2862, doi:10.1016/j.nuclphysb.2009.02.014.
  * [100] S. J. Brodsky, H.-C. Pauli, S. S. Pinsky, Quantum chromodynamics and other field theories on the light cone, Phys. Rept. 301 (1998) 299–486.  arXiv:hep-ph/9705477, doi:10.1016/S0370-1573(97)00089-6.
  * [101] S. D. Glazek, R. J. Perry, Special example of relativistic Hamiltonian field theory, Phys. Rev. D 45 (1992) 3740–3754.  doi:10.1103/PhysRevD.45.3740.
  * [102] S. J. Brodsky, G. F. de Teramond, H. G. Dosch, J. Erlich, Light-Front Holographic QCD and Emerging Confinement, Phys. Rept. 584 (2015) 1–105.  arXiv:1407.8131, doi:10.1016/j.physrep.2015.05.001.
  * [103] S. Meissner, A. Metz, M. Schlegel, K. Goeke, Generalized parton correlation functions for a spin-0 hadron, JHEP 08 (2008) 038.  arXiv:0805.3165, doi:10.1088/1126-6708/2008/08/038.
  * [104] J. C. Collins, Fragmentation of transversely polarized quarks probed in transverse momentum distributions, Nucl. Phys. B 396 (1993) 161–182.  arXiv:hep-ph/9208213, doi:10.1016/0550-3213(93)90262-N.
  * [105] A. Bacchetta, F. G. Celiberto, M. Radici, P. Taels, Transverse-momentum-dependent gluon distribution functions in a spectator model, Eur. Phys. J. C 80 (8) (2020) 733.  arXiv:2005.02288, doi:10.1140/epjc/s10052-020-8327-6.
  * [106] J. B. Kogut, D. E. Soper, Quantum Electrodynamics in the Infinite Momentum Frame, Phys. Rev. D 1 (1970) 2901–2913.  doi:10.1103/PhysRevD.1.2901.
  * [107] M. D. Schwartz, Quantum Field Theory and the Standard Model, Cambridge University Press, 2014.
  * [108] G. P. Lepage, S. J. Brodsky, Exclusive Processes in Perturbative Quantum Chromodynamics, Phys. Rev. D 22 (1980) 2157.  doi:10.1103/PhysRevD.22.2157.
  * [109] A. Bacchetta, G. Bozzi, M. G. Echevarria, C. Pisano, A. Prokudin, M. Radici, Azimuthal asymmetries in unpolarized SIDIS and Drell-Yan processes: a case study towards TMD factorization at subleading twist, Phys. Lett. B 797 (2019) 134850.  arXiv:1906.07037, doi:10.1016/j.physletb.2019.134850.
  * [110] S. Rodini, A. Vladimirov, Definition and evolution of transverse momentum dependent distribution of twist-three, JHEP 08 (2022) 031, [Erratum: JHEP 12, 048 (2022)].  arXiv:2204.03856, doi:10.1007/JHEP08(2022)031.
  * [111] J. M. M. Chávez, V. Bertone, F. De Soto Borrero, M. Defurne, C. Mezrag, H. Moutarde, J. Rodríguez-Quintero, J. Segovia, Accessing the Pion 3D Structure at US and China Electron-Ion Colliders, Phys. Rev. Lett. 128 (20) (2022) 202501.  arXiv:2110.09462, doi:10.1103/PhysRevLett.128.202501.