# Proton image and momentum distributions from light-front dynamics

E. Ydrefors<EMAIL_ADDRESS>Instituto Tecnológico de Aeronáutica, DCTA,
12228-900 São José dos Campos, Brazil  Institute of Modern Physics, Chinese
Academy of Sciences, Lanzhou 730000, China  T. Frederico<EMAIL_ADDRESS>Instituto Tecnológico de Aeronáutica, DCTA, 12228-900 São José dos Campos,
Brazil

###### Abstract

We apply a dynamical three-constituent quark light-front model to study the
proton. The dynamics is based on the notion of a diquark (bound or virtual) as
the dominant interaction channel, which parameterize a contact interaction
between the quarks in order to build the three-body Faddeev Bethe-Salpeter
equations for the valence state, and we focus on the totally symmetric part of
the wave function. The Dirac electromagnetic form factor is used to fix the
model parameters, and the valence wave function is obtained. From that we
investigate its Ioffe-time image, non-polarized longitudinal and transverse
momentum distributions, and the double momentum distribution.

Bethe-Salpeter equation, pion valence state, proton structure, Ioffe-time
distribution, momentum distribution

##  I Introduction

The complex nucleon wave function on the null-plane $(x^{+}=t+z=0)$ expressed
in the Fock space in terms of its constituent degrees of freedom, namely
quarks and gluons at a given scale $\mu$ and strongly interacting, ultimately
provides the image through the associated probability densities [1, 2, 3]. The
relevant degrees of freedom at the hadronic scale are the dressed
constituents, which carry the complex infrared (IR) physics, namely, the
confinement and spontaneous chiral symmetry breaking ($\chi$SB). In
particular, the dynamical $\chi$SB is reflected in the large dressing of the
light-flavored quarks, and also in the nucleon mass [3]. It is also well known
the large IR dressing of the gluon as computed with lattice QCD (for a recent
discussion of the gluon propagator in the Landau gauge see [4]).

The light-front (LF) wave function is an eigenstate of the mass squared
operator and the compatible operators $\vec{P}$ (momentum), $J^{2}$ (squared
angular momentum), $J_{z}$, and other compatible operators like the parity.
However, $P_{z}$, $\vec{J}_{\perp}$, and parity are not diagonal in the Fock
space, i.e. they contain the interaction [1, 2]. Although simple to state,
still today the connection between QCD in Euclidean space, their infrared (IR)
properties, with the LF wave function, and its Fock decomposition is yet a
challenge for our understanding, beyond the large momentum behavior that
thanks to the asymptotic freedom is well known, like the counting rules (see
e.g. [1]) as well as the ultraviolet (UV) behavior of the Fock amplitudes [5].

Ideally the LF nucleon wave function onto the null-plane should have an
infinite number of Fock-state components that evolves with the renormalization
scale $\mu$. The wave function can be decomposed into Fock components each one
associated with a probability amplitude
$\Psi_{n}(x_{1},\vec{k}_{1\perp},x_{2},\vec{k}_{2\perp},...;\mu)$ for $n\geq
3$ partons, which is invariant under LF kinematical boosts. The probability
corresponding to each Fock component is given by:

$\displaystyle
P_{n}(\mu)=\Bigg{\\{}\prod_{i=1}^{n}\int\frac{d^{2}k_{i\perp}}{(2\pi)^{2}}\int^{1}_{0}dx_{i}\Bigg{\\}}\,\delta\left(1-\sum_{i=1}^{n}x_{i}\right)$
(1)
$\displaystyle\times\delta\left(\sum_{i=1}^{n}\vec{k}_{i\perp}\right)|\Psi_{n}(x_{1},\vec{k}_{1\perp},x_{2},\vec{k}_{2\perp},...;\mu)|^{2}\,,$

where the transverse momentum of the constituent is $\vec{k}_{i\perp}$ and its
longitudinal momentum fraction $x_{i}$. We observe that the probability
$P_{n}(\mu)$ is invariant under LF kinematical boosts including translations
on the null-plane hypersurface, so that we can choose, in particular, the
frame where the transverse momentum vanishes.

The valence component corresponds to $n=3$. For simplicity, we have not
depicted the dependence on the polarization state of the nucleon, as well as
its constituents. The probability for each Fock component is $P_{n}(\mu)$,
where we kept the scale dependence. At the hadron scale $(\sim\Lambda_{QCD})$
the dominant component is the valence one, and for example in the pion case it
amounts to about 70% as it has been computed recently in a Bethe-Salpeter (BS)
framework [6]. The total normalization is $\sum_{i\geq
3}^{\infty}P_{n}(\mu)=1$.

Each Fock amplitude can be written in the configuration space associated with
the null-plane, where the three-dimensional position coordinates for each
constituent are $\\{b^{-}_{i},\vec{b}_{i\perp}\\}$, namely the light-like
coordinate $(b^{-}_{i}=t-z)$ and transverse position $(\vec{b}_{i\perp})$,
conjugate to $k^{+}_{i}=x_{i}p^{+}$ and $\vec{k}_{i\perp}$, respectively. The
Fock component of the wave function on the null-plane is obtained by a Fourier
transform:

$\displaystyle\tilde{\Psi}_{n}(\tilde{x}_{1},\vec{b}_{1\perp},...;\mu)=$ (2)
$\displaystyle\Bigg{\\{}\prod_{i=1}^{n}\int\frac{d^{2}k_{i\perp}}{(2\pi)^{2}}\int^{1}_{0}\frac{dx_{i}}{2\pi}\text{e}^{i\tilde{x}_{i}x_{i}-i\vec{b}_{i\perp}\cdot\vec{k}_{i\perp}}\Bigg{\\}}$
$\displaystyle\times\,\delta\left(1-\sum_{i=1}^{n}x_{i}\right)\,\delta\left(\sum_{i=1}^{n}\vec{k}_{i\perp}\right)\Psi_{n}(x_{1},\vec{k}_{1\perp}...;\mu)\,,$

where the dependence on the Ioffe time $\tilde{x}_{i}=b^{-}_{i}p^{+}$ [7, 8,
9] was given in the probability amplitude instead of the light-like position
onto the null-plane. The probability densities
$|\tilde{\Psi}_{n}(\tilde{x}_{1},\vec{b}_{1\perp},...;\mu)|^{2}$ build an
image of the nucleon on the null-plane, where the light-like coordinate, shows
the relevance of the Ioffe time to complete the image of the nucleon (see a
recent general discussion in [10] and in the case of the pion in [6]).
Importantly, the Ioffe-time is a Lorentz invariant quantity which is related
to the spatial distance between the struck quark and the spectators. Using the
Ioffe-time one can through the inverse Fourier transform construct frame-
independent PDFs. It should also be noted that the Ioffe-time representation
of the PDF can be related at small space-like separations to the so-called
Ioffe-time pseudo-distribution [11], which has been used to obtain the parton
distribution of the pion from Lattice QCD [12].

Different parton probability densities, namely one-, two- and N-body ones can
be defined given the LF wave function and reveal the multifaceted structure of
the nucleon, which are associated with different observables being of interest
not only for the present hadron facilities but also for the physics cases of
the future Electron Ion Collider [13]. In particular, we mention the
electromagnetic ones, as for example the elastic form factor, and the parton
distributions which are associated with one-body probability densities.

Other quantities which can be computed from the LF wave function are the
generalized parton distributions (GPDs), being the non-perturbative objects
entering the cross sections for deeply virtual Compton scattering (DVCS), and
the transverse momentum distributions associated with semi-inclusive deeply
inelastic scattering (SIDIS) [14]. The PDFs extracted from inclusive deep
inelastic scattering give only information about the longitudinal momentum
fraction of the parton, i.e. simply a one-dimensional view of the hadron. The
GPDs and transverse momentum distributions (TMDs) provide a more complete
image of the hadronic structure, in particular regarding the distribution of
spin and orbital momentum in hadrons. That also allows a three-dimensional
nucleon tomography in mixed position-momentum space (see e.g. [15]). However,
the most complete image is obtained through the six-dimensional Wigner
distributions, and their Fourier transforms that are related to the
generalized TMDs (GTMDs), which can appear in the representation of hard QCD
processes [16, 17, 18, 15]. GTMDs are associated to matrix elements of bi-
local partonic field operators with separation in all three light-front
coordinates defined onto the null-plane hypersurface. In general, they are
off-forward matrix elements between hadron states, which depend on the partons
longitudinal transverse momentum components. In particular, the GTMDs
correlate hadronic states with the same parton longitudinal momentum, namely
for vanishing skewness, and different relative transverse distance between the
struck partons’ initial and final states. It is worth to mention that, the
GTMDs contain the bi-local correlators that define both the GPDs, TMDs, PDFs
and as well as the electromagnetic form factor, which are obtained by taking
certain limits or performing integrations, see e.g. [15].

In particular, we remind that the space-like electromagnetic form factor can
be obtained from the celebrated Drell-Yan-West formula [1] using the ”plus”
component of the current, with momentum $q^{+}=0$ and
$q^{2}=-\vec{q}_{\perp}^{2}=-Q^{2}$, which is diagonal in the Fock space:

$\displaystyle F($ $\displaystyle
Q^{2})=\sum_{n=3}^{\infty}F_{(n)}(Q^{2})=\sum_{n=3}^{\infty}\sum_{j=1}^{n}e_{j}$
(3)
$\displaystyle\times\Bigg{\\{}\prod_{i=1}^{n}\int\frac{d^{2}k_{i\perp}}{(2\pi)^{2}}\int^{1}_{0}dx_{i}\Bigg{\\}}\,\delta\left(1-\sum_{i=1}^{n}x_{i}\right)\,$
$\displaystyle\times\delta\left(\sum_{i=1}^{n}\vec{k}^{\text{f}}_{i\perp}\right)\Psi_{n}^{\dagger}(x_{1},\vec{k}^{\text{f}}_{1\perp},x_{2},\vec{k}^{\text{f}}_{2\perp},...,x_{j},\vec{k}^{\text{f}}_{j\perp},...)$
$\displaystyle\times\Psi_{n}(x_{1},\vec{k}^{\text{i}}_{1\perp},x_{2},\vec{k}^{\text{i}}_{2\perp},...,x_{j},\vec{k}^{\text{i}}_{j\perp},...)\,,$

where the number of constituents in the Fock components are $n\geq 3$, and
$e_{j}$ is the constituent charge in units of the fundamental charge. The
partial contribution to the form factor from each Fock component of the wave
function is $F_{(n)}(Q^{2})$. In the adopted frame the pair-creation
contribution to the plus component of the current are suppressed, which is
important since the present model of the proton is limited to the valence
component, i.e. $n=3$ in (3). The quark momenta obtained via the LF boost from
the Breit-frame to the rest-frame of the initial (i) and final (f)
nucleon states are given by:

$\displaystyle\vec{k}^{\text{i}}_{i\perp}=\vec{k}_{i\perp}+\frac{\vec{q}_{\perp}}{2}x_{i},\quad\vec{k}^{\text{f}}_{i\perp}=\vec{k}_{i\perp}-\frac{\vec{q}_{\perp}}{2}x_{i};\quad
i\neq j\,$ (4)
$\displaystyle\quad\text{and}\quad\vec{k}^{\text{f(i)}}_{j\perp}=(-)\frac{\vec{q}_{\perp}}{2}(1-x_{j})-\sum_{i\neq
j}\vec{k}_{i\perp}\,,$

with the transverse momentum of the quark that absorbed the virtual photon
being: $\vec{k}_{j\perp}=\pm\frac{\vec{q}_{\perp}}{2}-\sum^{n}_{i\neq
j}\vec{k}_{i\perp}\,,$ with $+$ and $-$ meaning the momentum in the initial
and final hadron states, respectively. For each Fock component of the LF wave
function the transverse momenta add up to
$\sum_{j=1}^{n}\vec{k}^{\text{i}}_{j\perp}=\sum_{j=1}^{n}\vec{k}^{\text{f}}_{j\perp}=\vec{0}_{\perp}$,
for the rest-frames of the initial and final hadron states. The normalized
proton wave function gives $F(0)=1$.

Beyond the electromagnetic processes, proton-
proton collisions performed at the Large Hadron Collider in its high-
luminosity phase requires a detailed consideration of the nucleon structure
for the understanding of the observed data, associated with multiple parton
interactions (MPIs), which are required for the description of hadronic final
states (see e.g. [19]).

MPIs become more important for high-energy collisions as the parton flux
increases, while the parton momentum fractions decrease, as the nucleon
momentum is shared among more participants. Therefore, the search for new
physics demands the consideration of MPIs in the dedicated experimental
analysis (see the review book on Multiple Parton Interactions at the LHC
[20]). One example, is the double parton scattering (DPS) in hadron-hadron
collisions, where two independent hard-scattering processes happen between
partons from a parton pair in each hadron. It receives contributions from all
LF Fock-components of each hadron wave function, and such information is
encoded in the double parton distribution function [21].

The DPS cross section which depends on the double parton distribution
functions (dPDFs) contain contributions from all Fock-component of the wave
function, and it is written as [21]:

$\displaystyle D($ $\displaystyle
x_{1},x_{2},\vec{\eta}_{\perp})=\sum_{n=3}^{\infty}D_{n}(x_{1},x_{2},\vec{\eta}_{\perp})=$
(5)
$\displaystyle\sum_{n=3}^{\infty}\int\frac{d^{2}k_{1\perp}}{(2\pi)^{2}}\frac{d^{2}k_{2\perp}}{(2\pi)^{2}}\,\Bigg{\\{}\prod_{i\neq
1,2}\int\frac{d^{2}k_{i\perp}}{(2\pi)^{2}}\int^{1}_{0}dx_{i}\Bigg{\\}}\,$
$\displaystyle\times\delta\left(1-\sum_{i=1}^{n}x_{i}\right)\,\delta\left(\sum_{i=1}^{n}\vec{k}_{i\perp}\right)$
$\displaystyle\times\Psi_{n}^{\dagger}(x_{1},\vec{k}_{1\perp}+\vec{\eta}_{\perp},x_{2},\vec{k}_{2\perp}-\vec{\eta}_{\perp},...)$
$\displaystyle\times\Psi_{n}(x_{1},\vec{k}_{1\perp},x_{2},\vec{k}_{2\perp},...)\,,$

where $\vec{\eta}_{\perp}$ is the transverse momentum shift and for simplicity
we have not depicted the polarization states for neither the constituents nor
the hadron itself. The Fourier transform of
$D(x_{1},x_{2},\vec{\eta}_{\perp})$ in $\vec{\eta}_{\perp}$ gives the
probability of finding the constituents 1 and 2 with momentum fraction $x_{1}$
and $x_{2}$ at a relative distance in the transverse direction
$\vec{y}_{\perp}$ within the hadron state. Here $\vec{y}_{\perp}$ is the
Fourier conjugate of $\vec{\eta}_{\perp}$. The quantity
$D_{n}(x_{1},x_{2},\vec{\eta}_{\perp})$ is contribution from a given Fock-
component of the LF wave function to the double parton distribution.

In terms of operator product the double distribution reads [22]:

$\displaystyle\tilde{D}(x_{1},x_{2},\vec{y}_{\perp})=2p^{+}\int
dy^{-}\int\frac{dz^{-}_{1}}{2\pi}\frac{dz^{-}_{2}}{2\pi}$ (6)
$\displaystyle\times\text{e}^{i(x_{1}z^{-}_{1}+x_{2}z^{-}_{2})p^{+}}\sum_{\lambda}\langle
p,\lambda|\mathcal{O}(y,z_{1})\mathcal{O}(0,z_{2})|p,\lambda\rangle\,,$

which has been obtained for the nucleon by recent LQCD calculations for
different operator structures [22]. Despite such efforts, it is useful to
obtain the dPDFs at the nucleon scale and identifying properties of the LF
wave function, as for example using AdS/QCD approach [23] and LF constituent
quark models (see e.g. [24]). Among such efforts to detail the LF wave
function within a dynamical framework, we should mention the Basis Light-Front
Quantization applied to QCD [25] and recently used to study the nucleon [26].

Motivated by the above discussion, our goal in this work is to explore the
consequences of the relativistic LF three-body dynamics in the structure of
the valence state of the nucleon, by studying one- and two-quark momentum
distributions, where effectively the interaction is dominated by a strong
scalar quark-quark correlation. This model relies on the use of the contact
interaction between the constituents within the Faddeev Bethe-Salpeter
approach on the light-front [27, 28, 29, 30], and more recently the model was
extended beyond the valence state in Euclidean space [31] and in Minkowski
space [32, 33]. For practical applications, the profile of the IR momentum
dependence of the valence wave function, as for instance computed from the
transverse amplitude, obtained by both the LF truncation and the full four-
dimensional approach are essentially equivalent, once the bound state binding
energies are close [31], which will be enough for the present study. As a
note, we observe that short-range correlations between two quarks are present
in the model, in analogy with the nucleon-nucleon short-range correlations
(see e.g. [34]), which have also its counterpart in the relativistic three-
body wave function [33], with the proviso that the UV behavior has to be
viewed with caution as the scaling laws [5] from QCD are not built-in.

The model emphasizes the IR dynamics of constituent quarks with a dominant
scalar diquark correlation. Indeed, one main feature that the continuum
approaches to QCD have been teaching us is that two-quark correlations,
namely, diquarks, which are not asymptotic states, are known to play a
relevant role in the structure and dynamics of the nucleon (see e.g. [35] and
the recent review [36]).

We should remind that the successful Nambu-Jona-Lasinio model applied to
investigate phenomenological aspects of QCD in the IR region [37], embodies
the dynamical chiral symmetry breaking by producing massive constituent
($m\sim 300$ MeV) for the $u$ and $d$ quarks and pions/kaons as Goldstone
bosons, bringing in addition diquarks, with the favored one being the scalar
color antitriplet $([ud]^{\overline{3}_{c}}_{0^{+}})$ state. We see the
renewed interest from LQCD groups in determining the properties of diquarks in
a gauge invariant way [38] gives at the physical pion mass a difference of
319(1) MeV between the mass of the lightest diquark,
$[ud]^{\overline{3}_{c}}_{0^{+}}$, and the light antiquark, and a size of
about one fm. The low-energy diquark effective degree of freedom has also been
invoked to smooth the transition between the hadron to quark phases of dense
matter (see Ref. [39]).

The model adopted in this work considers a bound or a virtual state pole in
the quark-quark transition amplitude as the main dynamical characteristic,
which is in line with modern evidences of the relevance of the
$[ud]^{\overline{3}_{c}}_{0^{+}}$ state in the IR properties of the quark-
quark effective interaction within the nucleon. We aim to explore the proton
bound-state structure in terms of constituent quarks degrees of freedom by
calculating the valence LF wave function, where our focus is to study its
Ioffe time representation, as well as the different one- and two-quark
momentum distributions.

The rest of this work is organized as follows. A brief presentation of the LF
three-quark model is given in section II, containing the description of the
homogeneous LF Faddeev integral equation and numerical results for the vertex
Faddeev component. The results for the distribution amplitude and Ioffe-time
image of the proton are given in section III. The calculations of the valence
Dirac form factor of the proton are discussed in IV. The results for the
momentum distributions, namely valence parton distribution, valence double
parton distribution and transverse momentum distributions for single and two
quarks in the forward limit and integrated in longitudinal momenta are shown
in section V. The main points of our work are summarized in section VI. The
work is completed by two appendices: in Appendix A it is given the derivation
of the main dynamical integral equation of the model and in Appendix B the
adopted numerical method to solve it.

##  II Brief presentation of the LF three-quark model

The effective LF three-body model [27] that will be applied in the present
work to study the proton, was originated as an attempt to generalize
Weinberg’s infinite momentum frame realization of the two-boson BS equation
[40] to the three-body problem. Weinberg’s original proposal kept a close
relation of the three-dimensional dynamics in the light-front to the Minkowski
space BS equation for the bound state. Later on, within the framework of LF
quantization such an equation was realized to be the lowest-order equation
with the kernel expressing the coupling of the two- and three-particle LF
Fock-states (see e.g. [41]). The full equivalence of the two-body BS equation
and its LF three-dimensional representation have to take into account, besides
the valence component, an infinite number of Fock-states. In principle,
covariance under kinematical boosts is guaranteed even working with a finite
truncation of the Fock-space, however covariance under dynamical LF boosts,
which are non-diagonal in the LF Fock-space, requires the dynamics to involve
an infinite set of these states. Therefore, the BS equation has built-in
dynamically an infinite set of LF Fock-states. One possibility of projecting
the BS equation to LF was done using the quasi-potential approach [42], where
all the dynamics is buried in an effective interaction which contains the
virtual propagation of the system in an infinite number of Fock-states, in
close relation to the ”iterated resolvent method” [1] to reduce the QCD
dynamics in a hadron to its valence component.

At the three-body level, relevant to study the nucleon structure, the
counterpart of the Weinberg’s equation for the contact interaction was
proposed in [27]. It was performed the projection to the LF of the Minkowski
space Faddeev-Bethe-Salpeter (FBS) equation for the three-boson vertex keeping
only the valence contribution. In principle, the interaction kernel of the LF-
FBS equation contains contributions beyond the valence component, appearing as
effective LF three-body forces (see e.g. [31]). Indeed, the projection of the
Minkowski space FBS equation onto the LF was done via the quasi-potential
approach in [43] and further developed in [44, 45]. We outline in Appendix A
the derivations including the three-body valence wave function and LF
dynamical equation within the three-body BS framework.

As we have mentioned above, the model adopted to investigate the Ioffe time
representation of the wave function and also the double parton distribution,
is based on the contact interaction between the constituent quarks, where the
spin degree of freedom is not taken into account, as it is our aim to to study
the spatial non-polarized distribution of the quarks in the valence state. In
the model we consider only the totally symmetric momentum part of the the
colorless three-quark wave function corresponding to the valence nucleon
state, as we are interested for the time being on the investigation of the
properties associated with the momentum distributions and the image of the
nucleon onto the null-plane. The valence LF wave function is given by [33]:

$\displaystyle\Psi_{3}(x_{1},\vec{k}_{1\perp},x_{2},\vec{k}_{2\perp},x3,\vec{k}_{3\perp})=$
(7)
$\displaystyle\frac{\Gamma(x_{1},k_{1\perp})+\Gamma(x_{2},k_{2\perp})+\Gamma(x_{3},k_{3\perp},)}{\sqrt{x_{1}x_{2}x_{3}}(M_{N}^{2}-M^{2}_{0}(x_{1},\vec{k}_{1\perp},x_{2},\vec{k}_{2\perp},x_{3},\vec{k}_{3\perp}))},$

with $\Gamma(x_{i},k_{i\perp})$, where $k_{i\perp}=|\vec{k}_{i\perp}|$, being
the Faddeev component of the vertex function for the bound state,
$x_{1}+x_{2}+x_{3}=1$,
$\vec{k}_{1\perp}+\vec{k}_{2\perp}+\vec{k}_{3\perp}=\vec{0}_{\perp}$ and

$\displaystyle M^{2}_{0}(x_{1}$
$\displaystyle,\vec{k}_{1\perp},x_{2},\vec{k}_{2\perp},x_{3},\vec{k}_{3\perp})=$
(8)
$\displaystyle\frac{\vec{k}_{1\perp}^{2}+m^{2}}{x_{1}}+\frac{\vec{k}_{2\perp}^{2}+m^{2}}{x_{2}}+\frac{\vec{k}_{3\perp}^{2}+m^{2}}{x_{3}},$

is the free three-body squared mass for on-mass-shell constituents. The
factorized form of the valence wave function, namely with a vertex function
depending solely on the bachelor quark LF momenta, is a consequence of the
effective contact interaction between the constituent quarks, which is an
idealized model resembling the successful Nambu-Jona-Lasinio model applied to
model QCD [37]. It should be understood as an effective low-energy model which
is meant to have significance in the IR region where constituent quarks are
massive and bound forming the nucleon. We show in Appendix A the derivation of
the valence LF wave function starting from the three-legs Bethe-Salpeter
amplitude.

###  II.1 Homogeneous LF Faddeev integral equation

The Faddeev equation for the vertex component of the valence LF wave function
is given by [27, 30]:

$\displaystyle\Gamma(x,k_{\perp})=$
$\displaystyle\frac{\mathcal{F}(M^{2}_{12})}{(2\pi)^{3}}\int_{0}^{1-x}\frac{dx^{\prime}}{x^{\prime}(1-x-x^{\prime})}$
(9)
$\displaystyle\times\int_{0}^{\infty}d^{2}k^{\prime}_{\perp}\frac{\Gamma(x^{\prime},k^{\prime}_{\perp})}{\widehat{M}_{0}^{2}-M_{N}^{2}},$

where

$\widehat{M}_{0}^{2}=M^{2}_{0}(x,\vec{k}_{\perp},x^{\prime},\vec{k}^{\prime}_{\perp},1-x-x^{\prime},-(\vec{k}_{\perp}+\vec{k}^{\prime}_{\perp}))\,.$
(10)

For the sake of completeness, we provide in Appendix A a derivation of Eq. (9)
starting from the three-boson BS equation by projecting it onto the LF via the
quasi-potential technique.

The two-quark amplitude has the expression [33]:

$\mathcal{F}(M^{2}_{12})=\frac{\Theta(-M_{12}^{2})}{\frac{1}{16\pi^{2}y}\log\frac{1+y}{1-y}-\frac{1}{16\pi
ma}}\\\
+\frac{\Theta(M_{12}^{2})\,\,\Theta(4m^{2}-M_{12}^{2})}{\frac{1}{8\pi^{2}y^{\prime}}\arctan
y^{\prime}-\frac{1}{16\pi ma}}\,,$ (11)

with its argument, the effective off-shell mass of the two-quark subsystem
squared, given by

$\displaystyle
M^{2}_{12}=(1-x)M^{2}_{N}-\frac{k^{2}_{\perp}+(1-x)m^{2}}{x}\,,$ (12)
$\displaystyle y=\frac{M_{12}}{\sqrt{4m^{2}-M^{2}_{12}}}\quad\text{and}\quad
y^{\prime}=\frac{\sqrt{-M^{2}_{12}}}{\sqrt{4m^{2}-M^{2}_{12}}}\,.$

Additionally, in Eq. (11), $\Theta(x)$ denotes the Heaviside theta function.

The kernel of the LF Faddeev equation (9) contains the quark exchange
mechanism expressed by the presence of the three-quark LF resolvent, namely
the operator $[\widehat{M}^{2}_{0}-M_{N}^{2}]^{-1}$. Consistently with the
adopted model, it is well known [35] that the four-dimensional formulation of
the three-quark BS equation presents the quark exchange kernel, when the
diquarks dominates the quark-quark interaction. We should emphasize the
physical significance of the present model in the context of nucleon models
formulated commonly within the BS approach in Euclidean space [35]. Our model
provides directly the LF wave function allowing to access momentum
distributions and keeps the straight relation with the Bethe-Salpeter
framework, in contrast with commonly used Euclidean BS approaches.
Furthermore, it incorporates the main physics of more sophisticated Euclidean
formulations, as the quark exchange kernel, and a pole in the quark-quark
amplitude representing a bound or virtual diquark state.

The quark-quark scattering amplitude, $\mathcal{F}(M^{2}_{12})$, weights the
Faddeev LF integral equation and carries the pole of the bound or virtual
diquark, that depends on the scattering length $a$, which can be either
positive or negative. If $a>0$, the quark-quark system is bound and the
nucleon will be described as a quark-diquark system. On the contrary if $a$ is
negative no physical two-body bound-state exists and the nucleon is thus a
Borromean state. In both cases, $\mathcal{F}(M^{2}_{12})$ has a pole, for
$a>0$ in the physical complex-energy sheet and for $a<0$ in the
$2^{\text{nd}}$ sheet, meaning the virtual state. Therefore, in either case
the strong diquark correlation is present in the model and should be
interpreted as dominating the IR properties of the nucleon. Both of these two
cases will be investigated in this paper. An earlier study of the nucleon
performed with a truncated form of Eq. (9) was performed in [28].

In Appendix B it is explained the adopted numerical method to solve the
integral equation by using a bicubic spline expansion. The condition $F(0)=1$,
has been adopted to normalize the solution, where $F(Q^{2})$ is the valence
Dirac form factor, which will be discussed in Sec. IV.

Model | $m$ | $a$ | $M_{2}$ | $M_{N}$ | $r_{F_{1}}$
---|---|---|---|---|---
| [MeV] | $[m^{-1}]$ | [MeV] | [$m$] | [fm]
I | 317 | -1.84 | - | 2.97 | 0.97
II | 362 | 3.60 | 681 | 2.60 | 0.72
Table 1: Constituent mass in MeV, scattering length, diquark mass, and three-
body mass for the two considered models. Our diquark mass for model II is
slightly smaller than the scalar diquark mass of 691 MeV obtained in [46].
Also shown is the radius defined as $r_{F_{1}}=\hbar
c\sqrt{-6\frac{dF_{1}}{dQ^{2}}|_{Q^{2}=0}}$ and the corresponding experimental
value is $0.757$ fm [47].

###  II.2 Vertex Faddeev component

The structure of the three-quark valence state is encoded in the vertex
function $\Gamma(x,k_{\perp})$, which was computed with the two parameter sets
from Table 1. The constituent quark masses are 317 MeV (model I) and 362 MeV
(model II) to be compared with about 350 MeV from a recent LQCD calculation in
the Landau gauge [48]. We choose two possibilities for diquarks, namely an
unbound one for $a<0$ and a bound one for $a>0$, with a diquark mass of 681
MeV. These parameters are found by reproducing qualitatively the space-like
Dirac form factor up to about 1 GeV2, as it will be shown later on. We observe
that, the diquark mass of 681 MeV, which has a difference of 319 MeV with
respect to the quark mass, coincidentally matches the gauge invariant result
from the LQCD calculation [38] of 319(1) MeV at the physical pion mass.

The results for the vertex function $\Gamma(x,k_{\perp})$ are shown in Fig. 1
for models I (lower panel) and II (upper panel). The vertex function for both
models has characteristic transverse momentum around the IR scale of
$\sim\Lambda_{QCD}$, which drives the decreasing behaviour with $k_{\perp}$.
In addition, the vertex function peaks between $x\sim 0.35-0.4$, and the peak
evolves to somewhat larger values of $x$ with $k_{\perp}$, as a consequence of
the dominance of $k^{2}_{\perp}/x$ in the free squared mass operator, which
comes with the quark exchange kernel, as should be a general feature of the
diquark (bound or virtual) dominance in the quark-quark interaction.

It is seen in the upper panel of Fig. 1 that for model II with $a>0$, there is
one node around $x=0.8$ and also one for small $x$. As studied in detail in
Ref. [31], for $a>0$ the physical ground state of (9) is not the lowest energy
solution of the equation. That is, it exists another unphysical solution with
$M^{2}_{N}<0$. This state is the relativistic analog of the well-known Thomas
collapse in non-relativistic three-body systems with zero-range interaction
[49]. For example, at $1/(a\,m)=0.26$ one has $M_{N}^{2}=-69\>m^{2}$, so it is
a very deep state. In principle, it should be possible to remove this state by
a momentum cut-off $\Lambda$ of the order of 1 GeV, which would weaken the
interaction in the short-range region.

Figure 1: Vertex function, $\Gamma(x,k_{\perp})$, for the models I (lower
panel) and II (upper panel).

##  III Distribution amplitude and Ioffe-time image

###  III.1 Distribution amplitude

The distribution amplitude (DA), $\phi(x_{1},x_{2})$ for the nucleon is
defined as

$\displaystyle\phi($ $\displaystyle x_{1},x_{2},x_{3})=\int
d^{2}k_{1\perp}d^{2}k_{2\perp}d^{2}k_{3\perp}$ (13)
$\displaystyle\times\delta(\vec{k}_{1\perp}+\vec{k}_{2\perp}+\vec{k}_{3\perp})\Psi_{3}(x_{1},\vec{k}_{1\perp},x_{2},\vec{k}_{2\perp},x_{3},\vec{k}_{3\perp}),$

with $x_{3}=1-x_{1}-x_{2}$, and obeys the symmetry relation

$\phi(x_{1},x_{2},x_{3})=\phi(x_{2},x_{1},x_{3}).$ (14)

It gives the dependence of the wave function on the longitudinal momentum
fraction when the quarks share the same transverse position. In Fig. 2, the
calculated DA is shown for the two models considered in the present work. In
the figure the DA was normalized so that

$\displaystyle\int_{0}^{1}dx_{1}\int_{0}^{1-x_{1}}dx_{2}\int_{0}^{1}dx_{3}\delta(1-x_{1}-x_{2}-x_{3})$
(15) $\displaystyle\times\phi(x_{1},x_{2},x_{3})=1.$

The two different models give similar results with a slightly wider
distribution for model II, which reflects the wave function behaviour close to
$x_{i}\sim 0$. Further insight comes with the Fourier transform as discussed
in what follows.

Figure 2: Distribution amplitude as a function of $x_{1}$ and $x_{2}$, for the
model I (lower panel) and for model II (upper panel).

###  III.2  Ioffe-time image of the valence state

For the study of the space-time structure of the proton it is of interest to
obtain the wave function in terms of the Ioffe-times ($\tilde{x}_{1}$ and
$\tilde{x}_{2}$) and the transverse coordinates ($\vec{b}_{1\perp}$ and
$\vec{b}_{2\perp}$), which provides the image of the proton on the null-plane
$x^{+}=0$. Such study has been performed recently for the pion [6]. This is
accomplished through the Fourier transform of
$\Psi_{3}(x_{1},\vec{k}_{1\perp};x_{2},\vec{k}_{2\perp};x_{3},\vec{k}_{3\perp})$.
For simplicity, we consider here the particular case:

$\displaystyle\Phi(\tilde{x}_{1}$
$\displaystyle,\tilde{x}_{2})\equiv\tilde{\Psi}_{3}(\tilde{x}_{1},\vec{0}_{\perp},\tilde{x}_{2},\vec{0}_{\perp})=$
(16)
$\displaystyle\int_{0}^{1}dx_{1}\,\text{e}^{i\tilde{x}_{1}\,x_{1}}\int_{0}^{1-x_{1}}dx_{2}\int_{0}^{1}dx_{3}$
$\displaystyle\times\delta(1-x_{1}-x_{2}-x_{3})\,\text{e}^{i\tilde{x}_{2}\,x_{2}}\,\phi(x_{1},x_{2},x_{3})\,,$

where the configuration space wave function is computed at the origin
$\vec{b}_{1\perp}=\vec{b}_{2\perp}=\vec{0}_{\perp}$.

Figure 3: Upper panel: Squared modulus of the Ioffe-time distribution as a
function of $\tilde{x}_{1}$ and $\tilde{x}_{2}$, for the model I. Lower panel:
Squared modulus of the Ioffe-time distribution as a function of
$\tilde{x}_{1}$ for two fixed values of $\tilde{x}_{2}$, namely
$\tilde{x}_{2}$ = 0 (solid line) and $\tilde{x}_{2}$ = 10 (dashed line).
Results shown for the model I (blue line) and model II (red line).

In Fig. 3 we present our results for the squared modulus of the Ioffe-time
distribution given by Eq. (16). In the upper panel is shown the 3D plot of the
distribution in terms of the variables $\tilde{x}_{1}$ and $\tilde{x}_{2}$ for
model I. It is clear the preference of quarks to minimize the relative
distance in Ioffe time, as also observed along $\tilde{x}_{1}=\tilde{x}_{2}$.
The decrease along the just reflects the presence of the third quark that
recoils as the center of mass is at rest. Notably, there is no perceptible
difference in this plot between model I and II.

Then, in the lower panel of Fig. 3 we show for both models the Ioffe-time
distribution as a function of $\tilde{x}_{1}$ for two fixed values of
$\tilde{x}_{2}$, namely $\tilde{x}_{2}=0$ and $\tilde{x}_{2}=10$. It is seen
that the results obtained with the two parameter sets are almost identical for
$\tilde{x}_{1}<17$. In addition, we observe the equality between the
$|\tilde{\phi}(0,10)|=|\tilde{\phi}(10,0)|=|\tilde{\phi}(10,10)|$, which comes
from the permutation symmetry of the wave function:

$|\tilde{\phi}(\tilde{x},0)|=|\tilde{\phi}(0,\tilde{x})|=|\tilde{\phi}(\tilde{x},\tilde{x})|\,,$
(17)

which is a general characteristic of the model. This explains also the rough
flat behavior of $|\tilde{\phi}(\tilde{x}_{1},10)|$ for $0<\tilde{x}_{1}<10$
(dashed line).

For $\tilde{x}_{2}=10$ (dashed lines in Fig. 3) a sizable decrease of the
magnitude is observed at $\tilde{x}_{1}>10$ and for larger values an
oscillatory behavior is seen. This reflects the size of the proton of about 1
fm and a mass of 1 GeV, with their dimensionless product taking into account
the factor 1/2 from the adopted metric $\lambda^{-1}\sim 0.1$ with the
characteristic oscillatory pattern having a wave length in Ioffe time of about
10, which is observed in the figure. This is also roughly the dimensionless
scale, which governs the decrease of the wave function in the two quarks
relative separation in Ioffe time.

Both trends, namely oscillation and damping of the wave function are
essentially the same for models I and II, as their proton charge radius are
somewhat close, besides the same proton mass, as seen in Table 1. It suggests
that this general behavior should be quite model independent. Notably, a
similar qualitative behavior of the Ioffe-time distribution for the pion was
obtained in [6]. We observe an exponential damping of the probability density
with the relative separation between the Ioffe time of the two quarks, and the
damping is expected to be more sizable if confinement is incorporated as it is
effective at large distances.

##  IV  Valence Dirac form factor of the proton

In the three-body null-plane model, i.e. only taking into account the valence
contribution $(n=3)$ in the form factor formula (3), the Dirac form factor is
given by:

$\displaystyle F_{\text{1}}$
$\displaystyle(Q^{2})=\Bigg{\\{}\prod_{i=1}^{3}\int\frac{d^{2}k_{i\perp}}{(2\pi)^{2}}\int^{1}_{0}dx_{i}\Bigg{\\}}\,\delta\left(1-\sum_{i=1}^{3}x_{i}\right)\,$
(18)
$\displaystyle\times\delta\left(\sum_{i=1}^{3}\vec{k}^{\text{f}}_{i\perp}\right)\Psi_{3}^{\dagger}(x_{1},\vec{k}^{\text{f}}_{1\perp},...)\Psi_{3}(x_{1},\vec{k}^{\text{i}}_{1\perp},...),$

with $Q^{2}=\vec{q}_{\perp}\cdot\vec{q}_{\perp}$. In Eq. (4) the transverse
momentum of the quarks in the Breit frame are given.

One has in Eq. (18) that
$d^{2}k_{i\perp}=|\vec{k}_{i\perp}|d|\vec{k}_{i\perp}|d\theta_{i}$ ($i=1,2$)
with
$\vec{k}_{i\perp}\cdot\vec{q}_{\perp}=|\vec{k}_{i\perp}||\vec{q}_{\perp}|\cos\theta_{i}$.
Additionally, the needed magnitudes of the transverse momenta are given by

$\displaystyle\Bigl{|}\vec{k}^{\text{f(i)}}_{i\perp}\Bigr{|}^{2}$
$\displaystyle=\Bigl{|}\vec{k}_{i\perp}\pm\frac{\vec{q}_{\perp}}{2}x_{i}\Bigr{|}^{2}$
(19) $\displaystyle=\vec{k}_{i\perp}^{2}+\frac{Q^{2}}{4}x_{i}^{2}\pm
x_{i}|\vec{k}_{i\perp}||\vec{q}_{\perp}|\cos\theta_{i}\,,$

with $-$ for f and $+$ for i, in addition we have

$\displaystyle\Bigl{|}$
$\displaystyle\vec{k}^{\text{f(i)}}_{3\perp}\Bigr{|}^{2}=\Bigl{|}\pm\frac{\vec{q}_{\perp}}{2}(x_{3}-1)-\vec{k}_{1\perp}-\vec{k}_{2\perp}\Bigr{|}^{2}=$
(20)
$\displaystyle(1-x_{3})^{2}\frac{Q^{2}}{4}+\vec{k}_{1\perp}^{2}+2|\vec{k}_{1\perp}||\vec{k}_{2\perp}|\cos(\theta_{1}-\theta_{2})$
$\displaystyle+\vec{k}_{2\perp}^{2}\pm(1-x_{3})|\vec{q}_{\perp}|\bigl{(}|\vec{k}_{1\perp}|\cos\theta_{1}+|\vec{k}_{2\perp}|\cos\theta_{2}\bigr{)}\,.$

Figure 4:  Computed $F_{1}(Q^{2})$ (solid and dot-dashed lines) compared with
the empirical fit (dashed line) obtained in Ref. [50].

In Fig. 4, the computed Dirac form factor, $F_{1}(Q^{2})$, for the two
parameter sets listed in Table 1, is compared with the global fit to
experimental data by Ye et al [50]. It is seen that model II with $a=3.60/m$
gives a quite good agreement with the experimental data. This thus favors the
description of the nucleon as a quark-diquark system. The computed values of
the radius $r_{F_{1}}=\hbar c\sqrt{-6\frac{dF_{1}}{dQ^{2}}|_{Q^{2}=0}}$ for
the two considered models are also listed in Table 1. The model II gives a
radius of 0.72 fm which is about 5% lower than the experimental value of
$0.757$ fm [47] from the charge form factor. On the contrary, for the first
model with $a<0$ a rather large radius of 0.97 fm was obtained. We could have
attempted to fit the charge radius from $F_{1}$, by changing the constituent
quark mass and scattering length. However, we choose to keep the qualitative
reproduction of the form factor up to $Q^{2}\sim 1$ GeV2, which should be the
scale of our model.

##  V  Momentum distributions

###  V.1 Valence parton distribution

We study next the decomposition of the single parton distribution function
(PDF), obtained from the integrand of Eq. (18) of the Dirac form factor:

$\displaystyle{f}_{1}(x_{1})$
$\displaystyle=\frac{1}{(2\pi)^{6}}\int_{0}^{1-x_{1}}dx_{2}\int_{0}^{1}dx_{3}\delta(1-x_{1}-x_{2}-x_{3})$
(21) $\displaystyle\times\int
d^{2}k_{1\perp}d^{2}k_{2\perp}d^{2}k_{3\perp}\delta(\vec{k}_{1\perp}+\vec{k}_{2\perp}+\vec{k}_{3\perp})$
$\displaystyle\times|\Psi_{3}(x_{1},\vec{k}_{1\perp},x_{2},\vec{k}_{2\perp},x_{3},\vec{k}_{3\perp})|^{2}$
$\displaystyle=I_{11}+I_{22}+I_{33}+I_{12}+I_{13}+I_{23}.$

where the contributions to the PDF are defined for $i=1,2,3$ as:

$\displaystyle
I_{ii}=\frac{1}{(2\pi)^{6}}\int_{0}^{1-x_{1}}dx_{2}\int_{0}^{1}\frac{dx_{3}}{x_{1}x_{2}x_{3}}\delta(1-x_{1}-x_{2}-x_{3})$
(22) $\displaystyle\times\int
d^{2}k_{1\perp}d^{2}k_{2\perp}d^{2}k_{3\perp}\delta(\vec{k}_{1\perp}+\vec{k}_{2\perp}+\vec{k}_{3\perp})$
$\displaystyle\times\frac{[\Gamma(x_{i},\vec{k}_{i\perp})]^{2}}{(M_{N}^{2}-M^{2}_{0}(x_{1},\vec{k}_{1\perp},x_{2},\vec{k}_{2\perp},x_{3},\vec{k}_{3\perp}))^{2}}\,,$

and for $i\neq j$:

$\displaystyle
I_{ij}=\frac{2}{(2\pi)^{6}}\int_{0}^{1-x_{1}}dx_{2}\int_{0}^{1}\frac{dx_{3}}{x_{1}x_{2}x_{3}}\delta(1-x_{1}-x_{2}-x_{3})$
(23) $\displaystyle\times\int
d^{2}k_{1\perp}d^{2}k_{2\perp}d^{2}k_{3\perp}\delta(\vec{k}_{1\perp}+\vec{k}_{2\perp}+\vec{k}_{3\perp})$
$\displaystyle\times\frac{\Gamma(x_{i},\vec{k}_{i\perp})\Gamma(x_{j},\vec{k}_{j\perp})}{(M_{N}^{2}-M^{2}_{0}(x_{1},\vec{k}_{1\perp},x_{2},\vec{k}_{2\perp},x_{3},\vec{k}_{3\perp}))^{2}}\,.$

Due to the symmetries of the three-body wave function under exchange of
particles 2 and 3, it follows that $I_{22}=I_{33}$ and $I_{12}=I_{13}$.

Figure 5:  Contributions to the PDF for the models I (lower panel) and II
(upper panel)

The contributions to the PDF at vanishing $Q^{2}$ are presented for the two
considered models in the middle and upper panels of Fig. 4. The total PDF is
also shown in each panel with a thick solid line. For both models a maximum of
the PDF is seen at $x\approx 1/3$. As is seen the right panel, the model II
with a positive scattering length gives an almost flat behavior around $x=0.8$
for the PDF. Larger differences in the behavior of the contributions can also
be observed for this set of parameters.

Interesting to observe that all the contributions have about the same size,
and peaks around 0.35, despite we are measuring the PDF for the quark labeled
by 1 with momentum fraction $x_{1}$. More variation of the peaks position are
seen for $a>0$ where the vertex function has a node for model II (see Fig. 1),
while the interplay with the denominator of the wave function where the
smallest virtuality in the mass squared leads to fixed positions in all
contributions around $1/3$. The contribution from $I_{11}$ corresponding to a
configuration, where quark 1 is picked up while the pair of quarks interacts,
does not dominate, meaning that the symmetrization of the momentum component
of the wave function is crucial for the proton PDF.

###  V.2  Valence double parton distribution

Following Eq. (5), we write the valence contribution to the double quark
distribution function as:

$\displaystyle
D_{3}(x_{1},x_{2};\vec{\eta}_{\perp})=\frac{1}{(2\pi)^{6}}\int_{0}^{1}dx_{3}\delta(1-x_{1}-x_{2}-x_{3})$
(24) $\displaystyle\times\int
d^{2}k_{1\perp}d^{2}k_{2\perp}d^{2}k_{3\perp}\delta(\vec{k}_{1\perp}+\vec{k}_{2\perp}+\vec{k}_{3\perp})$
$\displaystyle\times\Psi_{3}^{\dagger}(x_{1},\vec{k}_{1\perp}+\vec{\eta}_{\perp};x_{2},\vec{k}_{2\perp}-\vec{\eta}_{\perp};x_{3},\vec{k}_{3\perp})$
$\displaystyle\times\Psi_{3}(x_{1},\vec{k}_{1\perp};x_{2},\vec{k}_{2\perp};x_{3},\vec{k}_{3\perp}).$

Our results for the DPDF calculated for $\vec{\eta}_{\perp}=\vec{0}_{\perp}$
are shown for the two considered models in Fig. 6. For this particular value
of transverse momentum $D_{3}(x_{1},x_{2},\vec{0}_{\perp})$ the double
distribution is the probability density for finding quarks with momentum
fraction $x_{1}$ and $x_{2}$. In the upper panel it is seen that for model II,
a strong suppression of the DPDF is seen for $x_{1}>0.6$ as for the PDF. The
model with $a<0$ gives a slightly more narrow DPDF. Observe the different
shapes of the boundaries of the double quark distribution, giving
complementary information with respect to the two-quark transverse momentum
distribution, which is sensitive to the size of the proton, as we are going to
discuss. The boundary for the higher probability density region for model I
has an isosceles triangle shape, while for model II it has an isosceles
trapezoid shape. The totally symmetric character of the wave function leads to
the symmetry properties of the boundaries.

Figure 6: Double quark distribution function,
$D_{3}(x_{1},x_{2},\vec{0}_{\perp})$, for the models I (lower panel) and II
(upper panel).

The triangular shaped boundary for model I, could be anticipated from the peak
of largest probability for $x_{i}\sim 0.3$, being visible for 0.2 and above in
Fig. 6, up to the boundary $x_{1}+x_{2}=1-x_{3}$. The trapezoid shaped
boundary observed for model II can be associated with the strong damping of
the PDF above $x\sim 0.6$ seen in Fig. 5 and to its peak around $x\sim 0.35$,
these two properties compete to provide the form seen in the upper panel of
Fig. 6. Model II corresponds to an excited state and the nodes appearing in
the vertex function for $x$ around 0.2 and 0.6 provides such peculiar boundary
form. What is noticeable is the sensitivity of the double PDF to the detail of
the vertex function, while the valence transverse distributions are
essentially sensitive to the size of the three quark configuration. Therefore,
it is quite interesting to see that radially excited states have its
particular imprints on the double quark distribution, as well as on the PDF.

###  V.3  Transverse momentum densities

The single quark transverse momentum distribution in the forward limit [18]
and integrated in the longitudinal momentum is associated with the probability
density to find a quark with momentum $k_{\perp}$:

$\displaystyle L$
$\displaystyle{}_{1}(k_{1\perp})=\frac{k_{1\perp}}{(2\pi)^{6}}\int_{0}^{1}dx_{1}\int_{0}^{1-x_{1}}dx_{2}\int_{0}^{1}dx_{3}$
(25)
$\displaystyle\times\delta(1-x_{1}-x_{2}-x_{3})\int_{0}^{2\pi}d\theta_{1}\int
d^{2}k_{2\perp}d^{2}k_{3\perp}$
$\displaystyle\times\delta(\vec{k}_{1\perp}+\vec{k}_{2\perp}+\vec{k}_{3\perp})|\Psi_{3}(x_{1},\vec{k}_{1\perp},x_{2},\vec{k}_{2\perp},x_{3},\vec{k}_{3\perp})|^{2}$
$\displaystyle\times\,|\Psi_{3}(x_{1},\vec{k}_{1\perp},x_{2},\vec{k}_{2\perp},x_{3},\vec{k}_{3\perp})|^{2}\,,$

and the two-quark one reads

$\displaystyle
L_{2}(k_{1\perp},k_{2\perp})=\frac{k_{1\perp}k_{2\perp}}{(2\pi)^{6}}\int_{0}^{1}dx_{1}\int_{0}^{1-x_{1}}dx_{2}\int_{0}^{1}dx_{3}$
(26) $\displaystyle\times\delta(1-x_{1}-x_{2}-x_{3})\int
d^{2}k_{3\perp}\delta(\vec{k}_{1\perp}+\vec{k}_{2\perp}+\vec{k}_{3\perp})$
$\displaystyle\times\int_{0}^{2\pi}d\theta_{1}\int_{0}^{2\pi}d\theta_{2}|\Psi_{3}(x_{1},\vec{k}_{1\perp},x_{2},\vec{k}_{2\perp},x_{3},\vec{k}_{3\perp})|^{2}.$

Figure 7: One-quark transverse momentum density vs $k_{1\perp}$ for the models
I $(a=-1.84/m)$ and II $(a=3.60/m)$.

Figure 8: Two-quark transverse momentum density versus $k_{1\perp}$ and
$k_{2\perp}$, for the model I (lower panel) and II (upper panel).

In Fig. 7 the single quark transverse momentum density from Eq. (25) is shown
for models I and II. As expected, for model I with $a=-1.84/m$, the momentum
distribution is narrower than for $a=3.6/m$, as the radius for the former case
is larger. The peak of the momentum distribution is about 0.08 GeV for model I
and 0.12 GeV for model II, reflecting the larger size of the proton in model I
compared to model II.

In Fig. 8 the two-quark transverse momentum density is shown for the model I
(lower panel) and model II (upper panel). The more compact configuration of
model II is reflected in the wider distribution, and the probability density
peak is consistent with Fig. 7 to what was observed for the quark
distribution.

##  VI  Summary

We study the Ioffe-time image, non-polarized longitudinal and transverse
momentum distributions, and the double momentum distribution of the proton
relying on a dynamical constituent quark light-front three-body model. The
dynamics is based on the prevalence of the scalar diquark channel in the
quark-quark interaction. We assume a minimal structure of quark-quark contact
interaction, resembling the Nambu-Jona-Lasinio model, with the simplified
assumption of factorization of the spin degree of freedom, and focusing on the
totally symmetric momentum component of the light-front wave function. We
study two possibilities, namely a diquark in a bound or a virtual state.

The three-body light-front Faddeev Bethe-Salpeter equations for the valence
state was solved in the presence of virtual (model I) or bound (model II)
diquark states, for positive and negative scattering lengths, respectively.
The contact interaction allows to simplify the integral equations for the
Faddeev components of the vertex function, which depend only on the spectator
quark longitudinal momentum fraction and transverse momentum. We have not used
a momentum cut-off in the model as originally introduced in [27] and adopted
the no-cutoff version [30]. This simplified dynamical model allowed to
investigate only non-polarized quantities.

The adopted dynamical light-front model has two parameters: the constituent
quark mass and scattering length. To determine these parameters the proton
mass was fixed to its experimental value, and the binding energy as well as
the scattering length were varied to have a qualitative fit of the Dirac form
up to about 1 GeV2. It reproduces the Dirac form factor radius somewhat close
to 1 fm, having the case of bound diquark a more compact configuration than
the case with the virtual diquark state. These two possibilities produces
quite different proton properties for the model with no cut-off. The bound
diquark produces a deep three-quark non-physical state $(M^{2}<0)$, and the
one associated with the proton is in this case an excited state, where the
vertex function has nodes in the $x$ dependence, while for the virtual diquark
the nucleon is the ground state of the model. These two distinct natures of
the valence state for model I and II allowed to study their observable
consequences in the different momentum distributions and image.

Specifically, we computed several proton non-
polarized quantities for the models (I) and (II):

  1. (i)

the distribution amplitude (DA); it corresponds to the probability amplitude
to find the quarks with given momentum fraction at the same transverse
position. The model I has a narrow distribution on the $(x_{1},x_{2})$ plane
when compared to model II.

  2. (ii)

the Ioffe-time image; it corresponds to the
Fourier-transform of the DA in the Ioffe time, and gives the probability
density of finding quarks along the light-like direction for quarks at the
same transverse position. The quarks tend to have close Ioffe-time positions,
and exhibit characteristic oscillations reflecting the size and mass of the
system, besides the symmetry of the configuration space wave function.

  3. (iii)

the quark distribution function; it corresponds to the probability density to
find a quark with a given momentum fraction at the nucleon scale, peaked
around $x\sim 0.3$, but distinguishing model I and II, the last one having
nodes in the spectator function, and presenting a more localized distribution.

  4. (iv)

the double quark distribution function; for $\vec{\eta}_{\perp}=0$; it
corresponds to the probability density to find the quarks with given momentum
fractions at the nucleon scale. The shape of the boundary with large
probability to observe the momentum fractions distinguish model I and II, with
an isosceles triangular and trapezoid shapes, for respectively, the ground-
state and excited-state configurations.

  5. (v)

the single quark transverse momentum density; is associated with the
probability density to find a quark with a given transverse momentum and
mainly sensitive to the size of the three-quark configuration and found peaked
around 0.1 GeV.

  6. (vi)

the double quark transverse momentum density; is associated with the
probability density to find the quarks with a given transverse momentum and,
again, mainly sensitive to the size of the three-quark configuration, a
smaller region in momentum is found for the larger size of the proton for
model I compared to the more compact configuration of the quarks in model II.

Future challenges for improving the nucleon effective model to be taken: the
computation of the Bethe-Salpeter amplitude in the four-dimensional Minkowski
space, which includes an infinite number of Fock-components, the introduction
of a cut-off and the spin degree of freedom, which we expect will provide more
insights into the nucleon structure.

###### Acknowledgements.

This work is a part of the project INCT-FNA proc. No. 464898/2014-5. This
study was financed in part by Conselho Nacional de Desenvolvimento Científico
e Tecnológico (CNPq) under the grant 308486/2015-3 (TF). E.Y. thanks for the
financial support of the grants #2016/25143-7 and #2018/21758-2 from FAPESP.
We thank the FAPESP Thematic Projects grants #13/26258-4 and #17/05660-0.

##  Appendix A Derivation of the LF wave function and vertex equation

The derivations made in this Appendix are based on the LF projection technique
of the BS equation and corresponding amplitude based on the Quasi-Potential
expansion developed in [42] for the two-boson problem and in Refs. [43, 44,
45] for the three-boson case. For the sake of completeness, we sketch here the
main steps in deriving the valence wave function (7) and the associated
Faddeev equation for the vertex function (9).

The starting point is the three-boson BS amplitude, which is defined as:

$\Psi_{M}(y_{1},y_{2},y_{3};p)=\langle
0\left|T(\varphi(y_{1})\varphi(y_{2})\varphi(y_{3}))\right|p\rangle\,,$ (27)

where $y_{i}$ is the space-time position of particle $i$, $\varphi(y)$ is the
bosonic field operator and $p$ the total momentum.

The valence LF wave function comes from the projection of the BS amplitude
onto the null-plane:

$\psi_{3}(\vec{k}_{1\perp},x_{1};\vec{k}_{2\perp},x_{2};\vec{k}_{3\perp},x_{3})=(p^{+})^{2}\,(x_{1}x_{2}x_{3})^{\frac{1}{2}}\\\
\times\chi_{3}(\vec{k}_{1\perp},x_{1};\vec{k}_{2\perp},x_{2};\vec{k}_{3\perp},x_{3})\,,$
(28)

with

$\chi_{3}(\vec{k}_{1\perp},x_{1};...)=\int
dk^{-}_{1}\,dk^{-}_{2}\,\Phi_{M}(k_{1},k_{2},k_{3};p)\,,$ (29)

where $\Phi_{M}$ is the momentum representation of the Minkowski space BS
amplitude $\Psi_{M}$. We have introduced the auxiliary LF amplitude $\chi_{3}$
for the convenience of the derivations done in what follows.

The elimination of the relative LF time for the three-body BS equation and
associated amplitude requires an integration over two independent momenta
$k^{-}$, due to four-momentum conservation, and we introduce the following
operation for a quantity defined in Minkowski space:

$\displaystyle|A$ $\displaystyle:=$ $\displaystyle\int
dk_{1}^{-}dk_{2}^{-}\langle k_{1}^{-}k_{2}^{-}|A,$ $\displaystyle A|$
$\displaystyle:=$ $\displaystyle\int
dk_{1}^{-}dk_{2}^{-}A|k_{1}^{-}k_{2}^{-}\rangle,$ (30)

with $A$ being a matrix element of an operator that has matrix elements which
are functions of two independent momenta after the center of mass motion is
factorized. With the above notation the LF wave function is:

$|\chi_{3}\rangle=\big{|}|\Phi_{M}\rangle=\big{|}G_{0}|\Gamma_{M}\rangle$ (31)

where the BS equation for the vertex function is

$\displaystyle|\Gamma_{M}\rangle$ $\displaystyle=$ $\displaystyle
VG_{0}|\Gamma_{M}\rangle.$ (32)

Explicitly the three-particle free Green’s function is given by

$\langle
k_{1}^{-},k_{2}^{-}|G_{0}|k_{1}^{\prime-},k_{2}^{\prime-}\rangle=\frac{-i}{(2\pi)^{2}}\frac{\delta(k_{1}^{-}-k_{1}^{\prime-})}{\hat{k}_{1}^{+}\hat{k}_{2}^{+}(p^{+}-{\hat{k}_{1}}^{+}-\hat{k}_{2}^{+})}\\\
\times\frac{\delta(k_{2}^{-}-k_{2}^{\prime-})}{(k_{1}^{-}-\hat{k}_{1on}^{-})(k_{2}^{-}-\hat{k}_{2on}^{-})(p^{-}-k_{1}^{-}-k_{2}^{-}-\hat{k}_{3on}^{-})}~{},$
(33)

where the hat means operator character and the on-minus-shell momentum
$\hat{k}^{-}_{ion}=(\hat{\vec{k}}^{2}_{i\perp}+m^{2}_{i})/\hat{k}_{i}^{+}$.
The momentum conservation applies to the kinematical components of the
momentum of particle 3, such that
$\hat{k}^{+}_{3}=p^{+}-\hat{k}^{+}_{1}-\hat{k}^{+}_{2}$ and the analogous
expression for the transverse components. By performing the LF projection
using Eq. (30), the free LF Green’s function becomes $g_{0}=|G_{0}|,$ being
the free light-front resolvent, explicitly written as:

$g_{0}(\underline{k}_{1},\underline{k}_{2})=\frac{i\theta(p^{+}-k_{1}^{+}-k_{2}^{+})\theta(k_{1}^{+})\theta(k_{2}^{+})}{k_{1}^{+}k_{2}^{+}(K^{+}-{k_{1}}^{+}-k_{2}^{+})}\\\
\times\frac{1}{p^{-}-k_{1on}^{-}-k_{2on}^{-}-(p-k_{1}-k_{2})_{on}^{-}}~{},$
(34)

where $\underline{k}_{i}\equiv\\{k^{+}_{i},\vec{k}_{i\perp}\\}$.

For the auxiliary Green’s function $\tilde{G}_{0}$ one makes the choice

$\widetilde{G}_{0}=G_{0}|g^{-1}_{0}|G_{0},$ (35)

and with that the four-dimensional BS equation for the vertex function is also
a solution of:

$|\Gamma_{M}\rangle=W\widetilde{G}_{0}|\Gamma_{M}\rangle=WG_{0}|g^{-1}_{0}|G_{0}|\Gamma_{M}\rangle\,,$
(36)

where the quasi-potential $W$ is given by

$W=V+V\Delta_{0}W$ (37)

with $\Delta_{0}=G_{0}-\widetilde{G}_{0}\,.$

The LF auxiliary amplitude turns out to be:

$|\chi_{3}\rangle=\big{|}G_{0}|\Gamma_{M}\rangle=\big{|}G_{0}WG_{0}\big{|}g^{-1}_{0}|G_{0}\big{|}\Gamma_{M}\rangle\,,$
(38)

and introducing the LF interaction

$w=g_{0}^{-1}|G_{0}WG_{0}|g_{0}^{-1}\,,$ (39)

and substituting in the expressions for $|\chi_{3}\rangle$:

$|\chi_{3}\rangle=g_{0}\,w|G_{0}\big{|}\Gamma_{M}\rangle=g_{0}\,w|\chi_{3}\rangle\,,$
(40)

where the LF vertex function is a solution of

$|\Gamma_{LF}\rangle=w\,g_{0}|\Gamma_{LF}\rangle\,.$ (41)

We introduce in what follows the Faddeev decomposition to solve the above
equation for the vertex function. The potential is built from the two-body
ones as

$V=\sum_{i=1}^{3}V_{i},\quad\text{and}\quad V_{i}=V_{(2)i}S^{-1}_{i}$ (42)

where $S_{i}$ is the propagator of the particle $i$ and $V_{(2)i}$ is the
interaction between the particles $j$ and $k$.

The Faddeev decomposition of $W$ reads

$W=\sum_{i}W_{i}\quad\text{and}\quad w=\sum_{i}w_{i},$ (43)

where $w_{i}=g_{0}^{-1}|G_{0}W_{i}G_{0}|g_{0}^{-1}\,,$ and formally the LF
three-body wave function is:

$|\chi_{3}\rangle=g_{0}\,\sum_{i}|v_{i}\rangle\quad\text{with}\quad|v_{i}\rangle=w_{i}|\chi_{3}\rangle\,,$
(44)

where $|v_{i}\rangle$ are the Faddeev components of the LF vertex function,
namely $|\Gamma_{LF}\rangle=\sum_{i}|v_{i}\rangle$.

In terms of the pairwise interaction, $V_{i}$, one has

$\displaystyle W_{i}$ $\displaystyle=V_{i}+V_{i}\Delta_{0}(V_{i}+V_{j}+V_{k})$
(45)
$\displaystyle+V_{i}\Delta_{0}(V_{i}+V_{j}+V_{k})\Delta_{0}(V_{i}+V_{j}+V_{k})+\cdots$

At the lowest (LO) order the effective potential is

$w^{\text{LO}}_{i}=g_{0}^{-1}|G_{0}V_{i}G_{0}|g_{0}^{-1}\,,$ (46)

and the Faddeev equations for the components of the vertex read

$\displaystyle|v^{\text{LO}}_{i}\rangle$ $\displaystyle=$ $\displaystyle
t^{\text{LO}}_{i}g_{0}(|v^{\text{LO}}_{j}\rangle+|v^{\text{LO}}_{k}\rangle)\,,$
(47)

with the LF T-matrix being a solution of

$t^{\text{LO}}_{i}=w^{\text{LO}}_{i}+w^{\text{LO}}_{i}g_{0}\,t^{\text{LO}}_{i}\,.$
(48)

Furthermore, it should be noted that the LF resolvent and potential are
immersed in the three-body system.

For the contact interaction the matrix element of the potential $V_{i}$ is

$\langle
k_{j},k_{k}|V_{i}|k^{\prime}_{j},k^{\prime}_{k}\rangle=\lambda(2\pi)^{2}\delta(k_{i}-k^{\prime}_{i})(k_{i}^{2}-m^{2})\,.$
(49)

By introducing it in Eq. (46), solving the LF T-matrix equation (48), and also
taken into that the two-body system is immersed in the three-body one, it is
found that:

$\Bigl{\langle}\underline{k}_{j},\underline{k}_{k}\Bigl{|}t^{\text{LO}}_{i}\Bigr{|}\underline{k}^{\prime}_{i},\underline{k}^{\prime}_{j}\Bigr{\rangle}=-i\mathcal{F}(M^{2}_{jk})k^{+}_{i}\delta(\underline{k}_{i}-\underline{k}^{\prime}_{i}),$
(50)

where $M^{2}_{jk}=(p-k_{i,on})^{2}$ and momentum conservation implies that
$\underline{p}=\underline{k}_{i}+\underline{k}_{j}+\underline{k}_{k}=\underline{k}^{\prime}_{i}+\underline{k}^{\prime}_{j}+\underline{k}^{\prime}_{k}\,.$

Furthermore, for the zero-range interaction we have that
$V_{i}\Delta_{0}V_{i}=0$ and it then follows that $W_{i}=V_{i}$ at the valence
order. We thus obtain

$v^{\text{LO}}_{i}(\underline{k}_{i})=-2i\mathcal{F}(M^{2}_{jk})\int
d\underline{k}^{\prime}_{j}k^{+}_{i}g_{0}(\underline{k}^{\prime}_{i},\underline{k}^{\prime}_{j})v^{\text{LO}}_{j}(\underline{k}^{\prime}_{j})\,,$
(51)

where the measure is
$d\underline{k}\equiv\frac{dk^{+}d^{2}k_{\perp}}{2(2\pi)^{3}}$ and the factor
of two comes from the symmetrization of the total vertex function with respect
to the exchange of the bosons, which also implies that

$v^{\text{LO}}_{i}(\underline{k}_{i})=v^{\text{LO}}_{j}(\underline{k}_{j})=v^{\text{LO}}_{k}(\underline{k}_{k})\,.$

Finally, identifying $\Gamma(x_{i},k_{i\perp})\equiv
v^{\text{LO}}_{i}(\underline{k}_{i})$, with $x_{i}=k^{+}_{i}/K^{+}$ one finds
Eq. (9) from Eq. (51).

##  Appendix B Numerical methods

In the present work the homogeneous integral equation (9) was solved for given
values of $a$ and $M_{3}=M_{N}/m$ by expanding the vertex function
$\Gamma(x,k_{\perp})$ in a bicubic basis, on the domain
$\Omega=I_{k_{\perp}}\times I_{x}=[0,\infty]\times[0,1]$, of the form

$\Gamma(x,k_{\perp})=\sum_{i=0}^{2N_{k_{\perp}}}\sum_{j=0}^{2N_{x}-1}A_{ij}S_{i}(k_{\perp})S_{j}(x).$
(52)

. In the calculations the intervals for $k_{\perp}$ and $x$ were partitioned
into $N_{k_{\perp}}$ and $N_{x}$ subintervals, respectively.

The equation (9) can then be turned into a generalized eigenvalue problem of
the form

$\sum_{i^{\prime}j^{\prime}}P_{iji^{\prime}j^{\prime}}A_{i^{\prime}j^{\prime}}=\lambda(M_{3})\sum_{i^{\prime}j^{\prime}}V_{iji^{\prime}j^{\prime}}A_{i^{\prime}j^{\prime}},$
(53)

where
$P_{iji^{\prime}j^{\prime}}=S_{i^{\prime}}(k^{(i)}_{\perp})S_{j^{\prime}}(x^{(j)})$
and $V_{iji^{\prime}j^{\prime}}$ is the
right-hand side of (9) with $\Gamma$ replaced by
$S_{i^{\prime}}(k^{(i)}_{\perp})S_{j^{\prime}}(x^{(j)})$. The value of the
three-body mass $M_{3}$ were found by iteratively solving the non-linear
equation $\lambda(M_{3})=1$ and the corresponding coefficients from the
solution of (53). The obtained solution for the vertex function was
subsequently normalized by imposing the condition $F(0)=1$, where $F(Q^{2})$
is the valence LF Dirac form factor which is discussed in Sec. IV.

## References

  * Brodsky  _et al._ [1998] S. J. Brodsky, H.-C. Pauli, and S. S. Pinsky, Quantum chromodynamics and other field theories on the light cone, Phys. Rep. 301, 299 (1998), arXiv:hep-ph/9705477 [hep-ph] .
  * Bakker  _et al._ [2014] B. Bakker, A. Bassetto, S. Brodsky, W. Broniowski, S. Dalley, T. Frederico, S. Głazek, J. Hiller, C.-R. Ji, V. Karmanov, _et al._ , Light-Front Quantum Chromodynamics: A framework for the analysis of hadron physics, Nucl. Phys. B Proc. Suppl. 251, 165 (2014).
  * Arrington  _et al._ [2021] J. Arrington, C. A. Gayoso, P. C. Barry, V. Berdnikov, D. Binosi, L. Chang, M. Diefenthaler, M. Ding, R. Ent, T. Frederico, and et al., Revealing the structure of light pseudoscalar mesons at the electron–ion collider, J. Phys. G: Nucl. Part. Phys. 48, 075106 (2021).
  * Li  _et al._ [2020] S. W. Li, P. Lowdon, O. Oliveira, and P. J. Silva, The generalised infrared structure of the gluon propagator, Phys. Lett. B 803, 135329 (2020).
  * Ji  _et al._ [2003] X.-D. Ji, J.-P. Ma, and F. Yuan, Generalized counting rule for hard exclusive processes, Phys. Rev. Lett. 90, 241601 (2003), arXiv:hep-ph/0301141 .
  * de Paula  _et al._ [2021] W. de Paula, E. Ydrefors, J. H. Alvarenga Nogueira, T. Frederico, and G. Salmè, Observing the
Minkowskian dynamics of the pion on the null-plane, Phys. Rev. D 103, 014002
(2021), arXiv:2012.04973 [hep-ph] .

  * Gribov  _et al._ [1966] V. Gribov, B. Ioffe, and I. Pomeranchuk, What is the range of interactions at high-energies,
Sov. J. Nucl. Phys. 2, 549 (1966).

  * Ioffe [1969] B. Ioffe, Space-time picture of photon and neutrino scattering and electroproduction cross-section asymptotics, Phys. Lett. B 30, 123 (1969).
  * Braun  _et al._ [1995] V. Braun, P. Gornicki, and L. Mankiewicz, Ioffe - time distributions instead of parton momentum distributions in description of deep inelastic scattering, Phys. Rev. D 51, 6036 (1995), arXiv:hep-ph/9410318 .
  * Miller and Brodsky [2020] G. A. Miller and S. J. Brodsky, The Frame-
Independent Spatial Coordinate $\tilde{z}$: Implications for Light-Front Wave
Functions, Deep Inelastic Scattering, Light-Front Holography, and Lattice QCD
Calculations, Phys. Rev. C 102, 022201 (2020),
arXiv:1912.08911 [hep-ph] .

  * Radyushkin [2017] A. Radyushkin, Quasi-parton distribution functions, momentum distributions, and pseudo-parton distribution functions, Phys. Rev. D 96, 034025 (2017), arXiv:1705.01488 [hep-ph] .
  * Joó  _et al._ [2019] B. Joó, J. Karpie, K. Orginos, A. V. Radyushkin, D. G. Richards, R. S. Sufian, and S. Zafeiropoulos, Pion valence structure from Ioffe-time parton pseudodistribution functions, Phys. Rev. D 100, 114512 (2019), arXiv:1909.08517 [hep-lat] .
  * Abdul Khalek and et al. [2021] R. Abdul Khalek and et al., Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report (2021), arXiv:2103.05419 [physics.ins-det] .
  * Diehl [2003] M. Diehl, Generalized parton distributions,
Phys. Rep. 388, 41 (2003), arXiv:hep-ph/0307382 [hep-ph] .

  * Lorcé and Pasquini [2013] C. Lorcé and B. Pasquini, Structure analysis of the generalized correlator of quark and gluon for a spin-1/2 target, JHEP 09, 138 (2013).
  * Meißner  _et al._ [2008] S. Meißner, K. Goeke, A. Metz, and M. Schlegel, Generalized parton correlation functions for a spin-0 hadron, JHEP 08, 038 (2008).
  * Meißner  _et al._ [2009] S. Meißner, A. Metz, and M. Schlegel, Generalized parton correlation functions for a spin-1/2 hadron, JHEP 08, 056 (2009), arXiv:0906.5323 [hep-ph] .
  * Lorcé  _et al._ [2011] C. Lorcé, B. Pasquini, and M. Vanderhaeghen, Unified framework for generalized and transverse-momentum dependent parton distributions within a 3Q light-cone picture of the nucleon, JHEP 05, 041 (2011).
  * Treleani and Calucci [2018] D. Treleani and G. Calucci, Inclusive and Exclusive Cross-Sections, Sum Rules, Adv. Ser. Direct. High Energy Phys. 29, 29 (2018), arXiv:1707.00271 [hep-ph] .
  * Bartalini and Gaunt [2019] P. Bartalini and J. R. Gaunt, eds.,  _Multiple Parton Interactions at the LHC_, Vol. 29 (WSP, 2019).
  * Blok  _et al._ [2011] B. Blok, Y. Dokshitzer, L. Frankfurt, and M. Strikman, Four-jet production at LHC and Tevatron in QCD, Phys. Rev. D 83, 071501(R) (2011).
  * G. S. Bali, M. Diehl, B. Gläßle, A. Schäfer and C. Zimmermann [2021] G. S. Bali, M. Diehl, B. Gläßle, A. Schäfer and C. Zimmermann, Double parton distributions in the nucleon from lattice QCD, JHEP 09 (2021), 106, arXiv:2106.03451 [hep-lat] .
  * Traini  _et al._ [2017] M. Traini, M. Rinaldi, S. Scopetta, and V. Vento, The effective cross section for double parton scattering within a holographic AdS/QCD approach, Phys. Lett. B 768, 270 (2017), arXiv:1609.07242 [hep-ph] .
  * Rinaldi and Ceccopieri [2019] M. Rinaldi and F. A. Ceccopieri, Double parton scattering and the proton transverse structure at the LHC, JHEP 09, 097 (2019), arXiv:1812.04286 [hep-ph] .
  * Vary  _et al._ [2010] J. Vary, H. Honkanen, J. Li, P. Maris, S. Brodsky, A. Harindranath, G. de Teramond, P. Sternberg, E. Ng, and C. Yang, Hamiltonian light-front field theory in a basis function approach, Phys. Rev. C 81, 035205 (2010), arXiv:0905.1411 [nucl-th] .
  * Mondal  _et al._ [2020] C. Mondal, S. Xu, J. Lan, X. Zhao, Y. Li,
D. Chakrabarti, and J. P. Vary (BLFQ Collaboration), Proton structure from a
light-front Hamiltonian, Phys. Rev. D 102, 016008 (2020).

  * Frederico [1992] T. Frederico, Null-plane model of three bosons with zero-range interaction, Phys. Lett. B 282, 409 (1992).
  * de Araújo  _et al._ [1995] W. R. B. de Araújo, J. P. B. C. de Melo, and T. Frederico, Faddeev null-plane model of the nucleon, Phys. Rev. C 52, 2733 (1995).
  * Suisso  _et al._ [2002] E. F. Suisso, J. P. B. C. de Melo, and T. Frederico, Relativistic dynamics of Qqq systems, Phys. Rev. D 65, 094009 (2002).
  * Carbonell and Karmanov [2003] J. Carbonell and V. A. Karmanov, Three-boson relativistic bound states with zero-range two-body interaction, Phys. Rev. C 67, 037001 (2003).
  * Ydrefors  _et al._ [2017] E. Ydrefors, J. H. Alvarenga Nogueira,
V. Gigante, T. Frederico, and V. A. Karmanov, Three-body bound states with
zero-range interaction in the Bethe–Salpeter approach, Phys. Lett. B 770, 131
(2017), arXiv:1703.07981 [nucl-th] .

  * Ydrefors  _et al._ [2019] E. Ydrefors, J. H. Alvarenga Nogueira,
V. A. Karmanov, and T. Frederico, Solving the three-body bound-state Bethe-
Salpeter equation in Minkowski space, Phys. Lett. B 791, 276 (2019),
arXiv:1903.01741 [hep-ph] .

  * Ydrefors  _et al._ [2020] E. Ydrefors, J. H. Alvarenga Nogueira,
V. A. Karmanov, and T. Frederico, Three-boson bound states in Minkowski space
with contact
interactions, Phys. Rev. D
101, 096018 (2020), arXiv:2005.07943 [hep-ph] .

  * Hen  _et al._ [2017] O. Hen, G. A. Miller, E. Piasetzky, and L. B. Weinstein, Nucleon-nucleon correlations, short-lived excitations, and the quarks within, Rev. Mod. Phys. 89, 045002 (2017).
  * Eichmann  _et al._ [2016] G. Eichmann, H. Sanchis-Alepuz, R. Williams,
R. Alkofer, and C. S. Fischer, Baryons as relativistic three-quark bound
states, Prog. Part. Nucl. Phys. 91, 1 (2016), arXiv:1606.09602 [hep-ph] .

  * Barabanov  _et al._ [2021] M. Y. Barabanov, M. A. Bedolla, W. K. Brooks, G. D. Cates, C. Chen, Y. Chen, E. Cisbani, M. Ding, G. Eichmann, R. Ent, J. Ferretti, R. W. Gothe, T. Horn, S. Liuti, C. Mezrag, A. Pilloni, A. J. R. Puckett, C. D. Roberts, P. Rossi, G. Salmé, E. Santopinto, J. Segovia, S. N. Syritsyn, M. Takizawa, E. Tomasi-Gustafsson, P. Wein, and B. B. Wojtsekhowski, Diquark correlations in hadron physics: Origin, impact and evidence, Prog. Part. Nucl. Phys. 116, 103835 (2021), arXiv:2008.07630 [hep-ph] .
  * Klevansky [1992] S. P. Klevansky, The Nambu-Jona-Lasinio model of quantum chromodynamics, Rev. Mod. Phys. 64, 649 (1992).
  * Francis  _et al._ [2021] A. Francis, P. de Forcrand, R. Lewis, and K. Maltman, Diquark properties from full QCD lattice simulations, arXiv:2106.09080 [hep-lat] (2021).
  * Fukushima  _et al._ [2021] K. Fukushima, J. M. Pawlowski, and N. Strodthoff, Emergent Hadrons and Diquarks, arXiv:2103.01129 [hep-ph] (2021).
  * Weinberg [1966] S. Weinberg, Dynamics at infinite momentum, Phys. Rev. 150, 1313 (1966).
  * Lepage and Brodsky [1980] G. Lepage and S. J. Brodsky,
Exclusive Processes in Perturbative Quantum Chromodynamics, Phys. Rev. D 22,
2157 (1980).

  * Sales  _et al._ [2000] J. H. O. Sales, T. Frederico, B. V. Carlson,
and P. U. Sauer, Light front Bethe-Salpeter equation, Phys. Rev. C 61, 044003
(2000), arXiv:nucl-th/9909029 [nucl-th] .

  * Marinho and Frederico [2008] J. A. O. Marinho and T. Frederico, Next-to-leading order light-front three-body dynamics, PoS LC2008, 036 (2008).
  * Frederico and Salmè [2011] T. Frederico and G. Salmè, Projecting the Bethe-Salpeter Equation onto the Light-Front and back: A Short Review, Few Body Syst. 49, 163 (2011), arXiv:1011.1850 [nucl-th] .
  * Guimarães  _et al._ [2014] K. S. F. F. Guimarães, O. Lourenço, W. de Paula, T. Frederico, and A. C. dos Reis, Final state interaction in $D^{+}\to K^{-}\pi^{+}\pi^{+}$ with $K\pi$ I = 1/2 and 3/2 channels, JHEP 08, 135 (2014), arXiv:1404.3797 [hep-ph] .
  * Ferretti [2019] J. Ferretti, Effective Degrees of Freedom in Baryon and Meson Spectroscopy, Few-Body Syst. 60, 17 (2019).
  * Xiong  _et al._ [2019] W. Xiong _et al._ , A small proton charge radius from an electron–proton scattering experiment, Nature 575, 147 (2019).
  * Oliveira  _et al._ [2019] O. Oliveira, P. J. Silva, J.-I. Skullerud, and A. Sternbeck, Quark propagator with two flavors of O(a)-improved Wilson fermions, Phys. Rev. D 99, 094506 (2019), arXiv:1809.02541 [hep-lat] .
  * Thomas [1935] L. H. Thomas, The Interaction Between a Neutron and a Proton and the Structure of ${\mathrm{H}}^{3}$, Phys. Rev. 47, 903 (1935).
  * Ye  _et al._ [2018] Z. Ye, J. Arrington, R. J. Hill, and G. Lee,
Proton and Neutron Electromagnetic Form Factors and Uncertainties, Phys. Lett.
B 777, 8 (2018), arXiv:1707.09063 [nucl-ex] .