Proton image and momentum distributions from light-front dynamics

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 paramaterizes 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.


I. INTRODUCTION
The complex nucleon wave function on the nullplane (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 µ 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 (χSB). In particular, the dynamical χ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 P (momentum), J 2 (squared angular momentum), J z , and other compatible operators like the parity. However, P z , J ⊥ , 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 nullplane should have an infinite number of Fock-state components that evolves with the renormalization scale µ. The wave function can be decomposed into Fock components each one associated with a probability amplitude Ψ n (x 1 , k 1⊥ , x 2 , k 2⊥ , ...; µ) for n ≥ 3 partons, which is invariant under LF kinematical boosts. The probability corresponding to each Fock component is given by: where the transverse momentum of the constituent is k i⊥ and its longitudinal momentum fraction x i . We observe that the probability P n (µ) 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 (µ), where we kept the scale dependence. At the hadron scale (∼ Λ 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 ∞ i≥3 P n (µ) = 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 , b i⊥ }, namely the light-like coordinate (b − i = t − z) and transverse position ( b i⊥ ), conjugate to k + i = x i p + and k i⊥ , respectively. The Fock component of the wave function on the null-plane is obtained by a Fourier transform: where the dependence on the Ioffe timex 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 |Ψ n (x 1 , b 1⊥ , ...; µ)| 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 [15][16][17][18]. 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 = − q 2 ⊥ = −Q 2 , which is diagonal in the Fock space: where the number of constituents in the Fock components are n ≥ 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 Breitframe to the rest-frame of the initial (i) and final (f) nucleon states are given by: with the transverse momentum of the quark that absorbed the virtual photon being: k j⊥ = ± q ⊥ 2 − n i =j k i⊥ , 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 n j=1 k i j⊥ = n j=1 k f j⊥ = 0 ⊥ , for the rest-frames of the initial and final hadron states. The normalized proton wave function gives F (0) = 1.
Beyond the electromagnetic processes, protonproton 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]: where η ⊥ 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 , η ⊥ ) in η ⊥ 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 y ⊥ within the hadron state. Here y ⊥ is the Fourier conjugate of η ⊥ . The quantity D n (x 1 , x 2 , η ⊥ ) is contribution from a given Fockcomponent of the LF wave function to the double parton distribution.
In terms of operator product the double distribution reads [22]: 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 oneand 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 nucleonnucleon 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 twoquark 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 ∼ 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] 3c 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] 3c 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] 3c 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 nondiagonal 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 threebody 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]: with Γ(x i , k i⊥ ), where k i⊥ = | k i⊥ |, being the Faddeev component of the vertex function for the bound state, 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.

A. Homogeneous LF Faddeev integral equation
The Faddeev equation for the vertex component of the valence LF wave function is given by [27,30]: where 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]: with its argument, the effective off-shell mass of the two-quark subsystem squared, given by Additionally, in Eq. (11), Θ(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 op- 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 quarkquark amplitude representing a bound or virtual diquark state.
The quark-quark scattering amplitude, 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 quarkdiquark 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, F(M 2 12 ) has a pole, for a > 0 in the physical complex-energy sheet and for a < 0 in the 2 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.

B. Vertex Faddeev component
The structure of the three-quark valence state is encoded in the vertex function Γ(x, k ⊥ ), which was computed with the two parameter sets from  [46].
Also shown is the radius defined as rF 1 = c −6 dF 1 dQ 2 | Q 2 =0 and the corresponding experimental value is 0.757 fm [47].

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 spacelike Dirac form factor up to about 1 GeV 2 , 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 Γ(x, k ⊥ ) 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 ∼ Λ QCD , which drives the decreasing behaviour with k ⊥ . In addition, the vertex function peaks between x ∼ 0.35 − 0.4, and the peak evolves to somewhat larger values of x with k ⊥ , as a consequence of the dominance of k 2 ⊥ /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 2 N = −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 Λ of the order of 1 GeV, which would weaken the interaction in the short-range region.

A. Distribution amplitude
The distribution amplitude (DA), φ(x 1 , x 2 ) for the nucleon is defined as with x 3 = 1−x 1 −x 2 , and obeys the symmetry relation 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 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 ∼ 0. Further insight comes with the Fourier transform as discussed in what follows.

B. 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 (x 1 andx 2 ) and the transverse coordinates ( b 1⊥ and b 2⊥ ), 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 Ψ 3 (x 1 , k 1⊥ ; x 2 , k 2⊥ ; x 3 , k 3⊥ ). For simplicity, we consider here the particular case: where the configuration space wave function is computed at the origin b 1⊥ = b 2⊥ = 0 ⊥ .
In Fig. 3 we present our results for the squared modulus of the Ioffe-time distribution given by Eq. (15). In the upper panel is shown the 3D plot of the distribution in terms of the variablesx 1 andx 2 for model I. It is clear the preference of quarks to minimize the relative distance in Ioffe time, as also observed along x 1 =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 ofx 1 for two fixed values ofx 2 , namelyx 2 = 0 andx 2 = 10. It is seen that the results obtained with the two parameter sets are almost identical for x 1 < 17. In addition, we observe the equality between the |φ(0, 10)| = |φ(10, 0)| = |φ(10, 10)|, which comes from the permutation symmetry of the wave function: which is a general characteristic of the model. This explains also the rough flat behavior of |φ(x 1 , 10)| for 0 <x 1 < 10 (dashed line). Forx 2 = 10 (dashed lines in Fig. 3) a sizable decrease of the magnitude is observed atx 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 λ −1 ∼ 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 I. 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: with Q 2 = q ⊥ · q ⊥ . In Eq. (4) the transverse momentum of the quarks in the Breit frame are given. One has in Eq. (17) that d 2 k i⊥ = | k i⊥ |d| k i⊥ |dθ i (i = 1, 2) with k i⊥ · q ⊥ = | k i⊥ || q ⊥ | cos θ i . Additionally, the needed magnitudes of the transverse momenta are given by with − for f and + for i, in addition we have In Fig. 4, the computed Dirac form factor, F 1 (Q 2 ), for the two parameter sets listed in Table I, 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 F1 = c −6 dF1 dQ 2 | Q 2 =0 for the two considered models are also listed in Table I. 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 ∼ 1 GeV 2 , which should be the scale of our model.

A. Valence parton distribution
We study next the decomposition of the single parton distribution function (PDF), obtained from the integrand of Eq. (17) of the Dirac form factor: where the contributions to the PDF are defined for i = 1, 2, 3 as: and for i = j:  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 ≈ 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.

B.
Valence double parton distribution Following Eq. (5), we write the valence contribution to the double quark distribution function as: (23) Our results for the DPDF calculated for η ⊥ = 0 ⊥ are shown for the two considered models in Fig. 6. For this particular value of transverse momentum D 3 (x 1 , x 2 , 0 ⊥ ) 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. The triangular shaped boundary for model I, could be anticipated from the peak of largest probability for x i ∼ 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 ∼ 0.6 seen in Fig. 5 and to its peak around x ∼ 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.

C.
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 ⊥ : and the two-quark one reads In Fig. 7 the single quark transverse momentum density from Eq. (24) 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 threebody model. The dynamics is based on the prevalence of the scalar diquark channel in the quark-quark interaction. We assume a minimal structure of quarkquark 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 nonpolarized 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 GeV 2 . 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 nonpolarized quantities for the models (I) and (II): (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.
(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.
(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 ∼ 0.3, but distinguishing model I and II, the last one having nodes in the spectator function, and presenting a more localized distribution.
(iv) the double quark distribution function; for η ⊥ = 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.
(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.
(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 Fockcomponents, the introduction of a cut-off and the spin degree of freedom, which we expect will provide more insights into the nucleon structure. 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 twoboson 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: where y i is the space-time position of particle i, ϕ(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: where Φ M is the momentum representation of the Minkowski space BS amplitude Ψ M . We have introduced the auxiliary LF amplitude χ 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: 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: where the BS equation for the vertex function is Explicitly the three-particle free Green's function is given by where the hat means operator character and the on- The momentum conservation applies to the kinematical components of the momentum of particle 3, such that k + 3 = p + −k + 1 −k + 2 and the analogous expression for the transverse components. By performing the LF projection using Eq. (A4), the free LF Green's function becomes g 0 = |G 0 |, being the free light-front resolvent, explicitly written as: g 0 (k 1 , k 2 ) = iθ(p + − k + 1 − k + 2 )θ(k + 1 )θ(k + 2 ) k + 1 k + 2 (K + − k 1 where k i ≡ {k + i , k i⊥ }. For the auxiliary Green's functionG 0 one makes the choice and with that the four-dimensional BS equation for the vertex function is also a solution of: where the quasi-potential W is given by The LF auxiliary amplitude turns out to be: and introducing the LF interaction and substituting in the expressions for |χ 3 : where the LF vertex function is a solution of 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 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 where w i = g −1 0 |G 0 W i G 0 |g −1 0 , and formally the LF three-body wave function is: where |v i are the Faddeev components of the LF vertex function, namely |Γ LF = i |v i . In terms of the pairwise interaction, V i , one has A19) At the lowest (LO) order the effective potential is and the Faddeev equations for the components of the vertex read with the LF T-matrix being a solution of 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 k j , k k |V i |k j , k k = λ(2π) 2 δ(k i −k i )(k 2 i −m 2 ) . (A23) By introducing it in Eq. (A20), solving the LF Tmatrix equation (A22), and also taken into that the two-body system is immersed in the three-body one, it is found that: where M 2 jk = (p−k i,on ) 2 and momentum conservation implies that p = k i + k j + k k = k i + k j + k k .
Furthermore, for the zero-range interaction we have that V i ∆ 0 V i = 0 and it then follows that W i = V i at the valence order. We thus obtain where the measure is dk ≡ dk + d 2 k ⊥ 2(2π) 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 , with x i = k + i /K + one finds Eq. (9) from Eq. (A25).

Appendix B: Numerical methods
In the present work the homogeneous integral equation (9) . In the calculations the intervals for k ⊥ and x were partitioned into N k ⊥ and N x subintervals, respectively. The equation (9) can then be turned into a generalized eigenvalue problem of the form where P iji j = S i (k (i) ⊥ )S j (x (j) ) and V iji j is the right-hand side of (9) with Γ replaced by S i (k (i) ⊥ )S j (x (j) ). The value of the three-body mass M 3 were found by iteratively solving the non-linear equation λ(M 3 ) = 1 and the corresponding coefficients from the solution of (B2). 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.