Residual mean field model of valence quarks in the nucleon

Abstract

We develop a non-perturbative model for valence parton distribution functions (PDFs) based on the mean field interactions of valence quarks in the nucleonic interior. The main motivation for the model is to obtain a mean field description of the valence quarks as a baseline to study the short range quark–quark interactions that generate the high x tail of PDFs. The model is based on the separation of the valence three-quark cluster and residual system in the nucleon. Then the nucleon structure function is calculated within the effective light-front diagrammatic approach introducing nonperturbative light-front valence quark and residual wave functions. Within the model a new relation is obtained between the position, $x_p$, of the peak of $xq_V(x)$ distribution of the valence quark and the effective mass of the residual system, $m_R$, in the form: $x_{p} \approx {1\over 4} (1-{m_R\over m_N})$ at starting $Q^2$. This relation explains the difference in the peak positions for d- and u-quarks through the expected difference of residual masses for valence d- and u-quark distributions. The parameters of the model are fixed by fitting the calculated valence quark distributions to the phenomenological PDFs. This allowed us to estimate the overall mean field contribution in baryonic and momentum sum rules for valence d- and u-quarks. Finally, the evaluated parameters of the non-perturbative wave functions of valence 3q-cluster and residual system can be used in calculation of other quantities such as nucleon form factors, generalized partonic and transverse momentum distributions.

1 Introduction

Understanding the dynamics of the valence quarks in the nucleon is one of the most interesting aspects of quantum chromodynamics (QCD). The problem is not only how the valence PDFs can be calculated from the first principles in QCD, but also how in general the three relativistic Fermions become bound. The special feature of valence quarks is that they define the baryonic quantum number of the nucleon and one expects them to be a bridge between baryonic spectroscopy and the QCD structure of the nucleon. Since baryonic spectroscopy adheres reasonably well to the SU(6) spin-flavor symmetry, one expects such a symmetry to be found also in the valence quark distributions in the nucleon and the problem is to understand how and in which parts of the valence quark momentum distribution this symmetry is broken.

Historically, the first attempts to model valence quark structure of nucleon were based on the non-relativistic picture of constituent quarks in which SU(6) symmetry was preserved for the constituents carrying approximately one third of the nucleon mass [1, 2]. These models were successful in describing the phenomenologies of baryonic spectroscopy. However, experimental extraction of the valence PDFs indicates that SU(6) is apparently violated, especially at large Bjorken x. Thus, one of the unresolved issues is how to reconcile the apparent success of SU(6) in baryonic spectroscopy and its breaking down in partonic distributions.

Another unique property of valence quarks is that their distribution weighted by Bjorken x has a well defined peak, even if the shape of the PDF is not an observable. The peaking of the momentum distribution for the bound system implies the importance of mean-field dynamics for the individual valence fermions in the nucleon, while the position of the peak is sensitive to the dynamical aspects of interaction in the mean field

The history of the modeling of partonic distributions is very rich, ranging from non-relativistic [1] and relativistic [3] constituent quark models, bag models [4,5,6], models combining partonic and pionic cloud picture of the nucleon [7,8,9] as well as models based on the di-quark picture of the nucleon [10,11,12,13,14] following from the spin-dependent part of the one-gluon exchange in the non-relativistic approximation [15]. All these models are non-perturbative in their nature and their predictions varied widely for the general characteristics of the valence quark distribution such as position of the peak and relative strength of the u- and d-quark distributions. While calculations based on lattice QCD reproduce the general characteristics of valence PDFs, their complexity does not always allow an understanding of the dominant mechanisms of interactions.

In this and the following papers we develop a new model for valence quark distributions of the nucleon based on the multi-quark correlation picture. The validity of such a model is based on the fact that even though the number of the quarks is not conserved in the nucleon, the number of valence quarks is effectively conserved and therefore it is possible to describe them in the framework used for the description of a bound system of finite number of fermions. This approach is similar to highly successful multi-nucleon correlation model of calculation of momentum distribution of nucleons in the nuclei [16,17,18,19,20,21] as well as the calculation of the momentum distribution in ultra-cold atomic Fermi gas with contact interactions (see e.g. [22, 23]). In our approach we consider three distinct interaction dynamics of valence quarks which are: mean field three valence quark (3q) cluster, two-quark and three-quark short range interactions, all of which are expected to dominate at different momentum fraction ranges of the valence quarks. We demonstrate that such a framework in the description of the valence quark dynamics brings a new insight and methodology in analyzing the partonic distributions in the nucleon.

Why modeling? The advantage of the proposed framework is that it creates a new ground for making predictions for different QCD processes involving nucleons, since in this case one can make unique predictions based on whether the process under study is dominated by the interaction of quarks in the mean field or in two-/three-quark correlations. As experience from nuclear physics shows, eventually the ab-initio calculations [24, 25] (lattice calculation in the case of QCD) will reproduce the phenomena observed based on the model description of the processes [26,27,28]. However, ab-initio calculations, because of their general approach, are not always well positioned for making experimentally verifiable specific predictions.

The article is organized as follows: in Sect. 2 we first present the justification for the mean field, two- and three-quark short range interaction picture of partonic distributions. We then discuss the phenomenology of partonic distributions, which are almost universal for the all recent PDFs extracted from the analysis of different high energy scattering data of electro-production and weak interaction processes (e.g. deep inelastic scattering and Drell–Yan). Our focus is on the dynamical characteristics of the x weighted valence PDFs, such as the position and height of the peaks as well as the ratio of d- to u-quark distributions. We examine how the position and the height of the peak changes due to QCD evolution while the difference between the peak positions for u- and d-quark distribution is largely $Q^2$ independent. Another observation is the approximate validity of SU(6) symmetry in the region close to the peak of the partonic distributions. The latter observation is used as a justification for the assumption of approximate SU(6) symmetry for mean-field valence quarks. In Sect. 3 the mean field model of valence quark distributions is presented. Based on this model we calculate the $x q_V(x,Q^2)$ in leading order (LO), which is then fitted to the phenomenological distributions (in Sect. 4) to evaluate the parameters defining the non-perturbative wave functions. The latter allowed us to ascertain the expected strength of the two- and three-quark correlation effects in the normalization and momentum sum rule of valence d- and u-quark distributions. The evaluated parameters of the wave functions can also be used for estimation of mean field contributions of the different QCD processes sensitive to the valence quark dynamics at moderate Bjorken x $\approx 0.2$. Sect. 5 summarizes our results, discusses the predictions that follow from the model and the possibility of their experimental verification. We also discuss the limitations of the model and outline the second part of the work in which the mechanism of quark–quark short range interaction is included to calculate the high momentum component of valence quark distributions. In Appendix A we summarize the Light-Front effective diagrammatic rules on which the calculation in the paper is based. In Appendix B we present the mathematical details of the derivation related to the integration in the transverse momentum space.

2 The model of valence quark distributions and phenomenology of PDFs

The uniqueness of the valence quarks as carriers of the baryonic quantum number of nucleons is that even if the number of quarks in the bound nucleon is not conserved their “effective” number is conserved. This situation allows us to introduce a concept of potential energy and apply a rather well known theoretical framework for the description of a bound system consisting of a finite number of fermions. In this framework the bulk of the momentum distribution is defined by the mean field dynamics while, if interaction is short range, the fermion–fermion short-range correlations are responsible for the high momentum part of the same distribution. As it was mentioned in the introduction, this approach has been successfully applied in nuclear physics and physics of ultra-cold atomic Fermi gases.

Fig. 1

The mean field (a), two-valence (b) and three-valence (c) quark short-range correlation contributions to the deep-inelastic scattering off the nucleon in the partonic picture

2.1 Mean field and quark correlation model of valence quark distributions in the nucleon

Any finite size, bound system with a fixed number of fermions with nonzero interaction length will exhibit a mean field interaction of its constituents [29]. Hereafter, we call a mean field an approximation in which the interaction of a given constituent with other constituents in the system is characterized by a single effective potential representing sum of all possible non-perturbative interactions between constituents. (Such a condition is similar to one in which Hartree–Fock approximation is justified.) In quantum mechanics if the binding of a constituent is due to the minimum of the mean field potential, the latter can be approximated by a harmonic oscillator potential near the vicinity of the minimum. Thus the momentum distribution of the constituent in such a mean field can be described in exponential form with the coefficient of the exponent characterizing the size of the system (or orbit). The same mean field dynamics define also the characteristic average momentum carried by the constituents which can be related to the position of the peak of the momentum weighted distributions. For the system in which interaction strength between constituents is not negligible at distances much smaller than the bulk size of the system (like the contact interaction, see e.g. [23]), the high momentum part of the momentum distributions is a result of the short range pair-wise interactions between constituents, often referred to as short-range correlations.

In the present work we apply such a framework for the calculation of the light-front momentum distribution of valence quarks in the nucleon [30]. The approach is based on the argument, that the “effective” number of valence quarks is conserved and the strength of the quark–quark interactions is sufficiently large in the volume they occupy. It is assumed that the volume the valence quarks occupy is smaller than the actual size of the nucleon. Thus we arrive at the picture of a non-perturbative valence quark cluster embedded in the nucleon. We assume the following interaction dynamics for valence quarks in the nucleon: the non-perturbative interaction among the three valence quarks provides the confinement (referred hereafter as mean field interaction), the short range perturbative quark–quark interactions through hard gluon exchanges generate the high momentum (high Bjorken x) tail of the valence quark distribution (hereafter, referred to as quark–quark short range correlations) and the interaction of three-valence quark cluster with the residual nucleon system, R, influences the final momentum distribution of valence quarks in the nucleon.

Such a framework in leading order, is represented by the three light-cone time-ordered diagrams in Fig. 1 and it is assumed that they should describe the major properties of partonic distribution at Bjorken $x>0.1$. For the smaller x the valence quark distributions are defined predominantly by Regge dynamics dominated by $\alpha = {1\over 2}$ trajectory, which is beyond the scope of this paper (see e.g. Ref.  [31]).

Before proceeding with the calculations within this framework, one can evaluate the valence quark distributions parametrically to validate the model at least on a qualitative level. First, the mean-field dynamics is largely responsible for the bulk properties of the bound three-quark system such as its size and average momentum of its constituents. The functional form of the latter can be described by an exponential function with the exponent being proportional to the spatial extension of the valence quarks.

To identify the range of Bjorken x that one expects to be generated by mean field dynamics we note that the phenomenology of PDFs is well consistent with the average transverse momentum of valence quarks being about $\langle k_t\rangle \sim 200$ MeV/c [32]. Using the approximate spherical symmetry of momentum distribution, we estimate that the Bjorken x range relevant to the mean-field to be $x\lesssim {2\langle k_t\rangle \over m_N} \sim 0.4$. Thus for massless valence quarks one expects an exponential type distribution for the range of $0.1 \lesssim x\lesssim 0.4$.

For the highest end of the Bjorken x range, the phenomenological observation is that for $x\ge 0.8$ the valence PDF’s behave as $(1-x)^3$ which is consistent with the dynamics of two hard gluon exchanges between three valence quarks [33,34,35] . Hereafter we refer them as three-quark (3q) correlations. The short range nature of such a correlation will require that all three quarks have momenta exceeding the above mentioned $\langle k_t\rangle$ and the interacting quark balances the other two spectator quarks. The latter results in $x > rsim {4\langle k_t\rangle \over m_n}\sim 0.8$.

Establishing the regions where one expects the dominance of mean-field and 3q-correlations our next conjecture is that the transition between these two regions happens through the two-quark short range interactions which we will refer to as 2q-correlation. If such two-quark correlations are due to short range gluon interactions then they are dominated among two quarks with opposite helicites resulting in a functional form of the partonic distribution $(1-x/B)^2$, where B is a parameter characterizing the total momentum fraction carried by the center of mass of the 2q-correlated pair.

Fig. 2

The graph shows the comparison of the x weighted valence u-quark distribution fitted by three analytic functions corresponding to the mean field, 2q- and 3q-correlations

In Fig. 2 we fitted functional forms following from the above discussed scenarios of mean-field, 2q- and 3q-correlations to the x weighted valence u-quark distribution. As the figure shows the exponential mean field as well as $(1-{x\over B})^2$ and $(1-x)^3$ forms for 2q- and 3q-correlations reproduce the shape of the valence quark distribution surprisingly well starting at $x > rsim 0.1$. Note that in principle the same functions should reproduce the $x<0.1$ part of the distributions, however, as it was mentioned earlier, starting at $x<0.1$ valence quark distributions are dominated by Regge behavior which is not included in the framework considered here.

2.2 Phenomenological properties of valence quark distributions

One of the most distinguished properties of valence quarks is that while their distributions do not exhibit any peaking structure, their x weighted distribution, $h(x,t)\equiv x q_V(x,t)$ has a clear peak at around $x\sim 0.2$ which makes them qualitatively different from the sea quark distributions (Fig. 3). Note that the peaking property is universal for leading and higher order approximation despite the shape of PDFs not being a physical observable with certain degree of arbitrariness in their definition. Here we define $t = \log {Q^2}$.

Fig. 3

The up and down valence PDFs (left) and x weighted PDFs (right) at $Q^2 = 4~\hbox {GeV}^2$. The 68% confidence level error bands for PDF parametrization are shown. The PDF sets CJ15nlo and NNPDF3.1-nnlo [36, 37] are used in the calculation. Similar features are also observed for other modern PDF parameterizations

One expects that the $Q^2$ dependence of the position and the height of the peak originate from QCD evolution, in which case with an increase of $Q^2$ the gluon radiations of valence quarks shift the peak towards smaller x and diminishes the height of the maximum. This can be seen in Fig. 4 where the valence d- and u-quarks structure functions are presented for $Q^2$ ranging from 5 to 100 GeV$^2$.

Fig. 4

The x dependence of $x q_V(x,Q^2)$ function for d- and u-valence quarks, calculated for CT14nnlo [38] PDFs at various $Q^2$

In the recent work [39] we found that the correlation between the height of the maximum of the $h(x_p,t)$ function and its position, $x_p$ has a universal exponential form: $h(x_p,t) = C e^{Dx_p}$, where constants, C and D saturate with the increase of the order of approximation in coupling constant. As it can bee seen from Fig. 5 the exponential behavior extends to practically the entire range of the coupling constant $\alpha$ being covered in various deep-inelastic processes.

The observed “height-position” correlation results in the following analytic form of the valence PDFs at the vicinity of $x_p$ [39] :

$$\begin{aligned} h(x,t) \equiv x q_V(x,t) \approx C + CDe x(1-x)^{1-x_p(t) \over x_p(t)}. \end{aligned}$$

(1)

It can be checked that for the above relation the h(x, t) function peaks at $x =x_p$. Equation (1) indicates that the exponent of the $(1-x)$ term, ${1-x_p(t) \over x_p(t)}$, is defined by the position of the peak $x_p$, which changes continuously with t due to QCD evolution. However as Fig. 4 (see also Fig. 7b (below) shows, already at starting $Q_0^2=1.69$ GeV$^2$ (before the onset of QCD evolution) the peak-positions $x_p$ are different for valence d and u quarks, resulting in different x-dependencies of PDFs at the vicinity of peaks. The observation that the exponent of $(1-x)$ term is flavor-dependent indicates a more complex dynamics in the generation of valence PDFs at $x\sim x_p$ than one expects, for example, from the mechanism of a fixed number gluon exchanges that will result in a same exponent for valence u- and d-quark PDFs. A flavor-dependent exponent can be generated by considering an effective nonperturbative potential in Weinberg type equations for relativistic bound states [40].

Fig. 5

The peak height and position correlation described by an exponential function: $h(x_p,t) =Ce^{Dx_p}$, from Ref.  [39]. CT14lo and CT14nnlo corresponds to the CT14 valence PDF parameterization in leading and next to next to leading order approximation [38]. Other PDF parameterizations also result in an exponential form of the peak height-position correlation. The CT14nnlo is the Hessian error at 68% confidence level, while the CT14lo shows just the central value

Since the mean-field part of the valence quark wave function will dominate in the calculations related to the bulk (integrated) characteristics of baryons (e.g. mass spectrum [1]), then to explain the approximate SU(6) symmetry of baryonic spectroscopy one should expect the same symmetry to hold for the mean-field part of the valence PDFs. For SU(6) symmetry one expects that the ratio of valence d- to u-quark distributions to scale as 0.5. As Fig. 6 shows all phenomenological valence quark distributions produce a $d_V/u_V$ ratio close to 0.5 in the region, $x\sim x_p \sim 0.2$, where one expects the dominance of mean-field dynamics. It is quite intriguing that some of these distributions such as CJ12 [41] and MMHT2014 [42] produce an almost constant scaling for $d_V/u_V$ ratios close to 0.5 at $x\approx x_p$.¹ The figure also shows the decrease of the $d_V/u_V$ ratio with an increase of x (>0.3), which is consistent with the increasingly large violation of SU(6) symmetry. One such mechanism that can result in such a violation is the onset of quark-correlation dynamics in which case one expects the violation of SU(6) symmetry due to the helicity selection rule in quark–quark interactions through the hard vector (gluon) exchanges (see e.g. Ref.  [34]).

Fig. 6

The ratio of the d- to u-valence quark distributions for various PDF sets [36, 38, 41,42,43] calculated at the lowest point of $Q_0^2 = M_c^2 = 1.69$ GeV$^2$ for which all the PDFs are defined

In Fig. 7a we show the $Q^2$ dependence of the $d_V/u_V$ ratio estimated at the peak position of the $h(x,t) = x\cdot d_V(x,Q^2)$ distribution. As the figure shows the ratio $\approx 0.5$ is largely $Q^2$ independent, which indicates that such a ratio reflects an underlying property of the quark wave function of the nucleon rather than being due to QCD evolution. The QCD evolution only slightly modifies this ratio due to the migration of possible quark correlation effects from large to small x region.

Finally, another phenomenological feature is that for valence d-quarks the position of the peak of the h(x, t) function is systematically lower than that of the valence u-quarks (see Fig. 3b). If such a difference in the peak positions is due to the dynamical property of valence d- and u-quarks in the nucleon, then one expects that it should be insensitive to QCD evolution and be largely independent on $Q^2$. Such an expectation is justified by the weak $Q^2$ dependence of the ratio of peak-positions of valence u- and d-quark distributions for CT class of PDF parameterizations shown in Fig. 7b. The CJ15 parameterization shows $Q^2$ independence starting at $Q^2\approx 10$ GeV$^2$ while the MMHT14 PDFs reach $Q^2$ independence at $Q^2 > 10^3$ GeV$^2$. Even though the magnitude of the gap between peaks of $xd_V(x,t)$ and $xu_V(x,t)$ varies noticeably between different groups of PDFs its $Q^2$ dependence is largely flat, not showing any systematic $Q^2$ dependence.

Fig. 7

a The $Q^2$ dependence of the ratio of d- to u-valence quark distributions at the peak position of h(x, t) distribution for d quark. b The $Q^2$ dependence of the ratio of the peak positions of the d- and u-valence quark h(x, t) distributions. Different PDF sets are used from Refs.  [36, 38, 41,42,43]

Summarizing this section we would like to emphasize the dedicated study of the above discussed characteristics of valence quark PDFs, such as exponential correlation of Fig. 5, the approximate 0.5 scaling of d- to u-valence quark PDF ratios at $x\sim x_p\sim 0.2$ (Figs. 6 and 7a) as well as the gap between d- and u-valence quark peaking positions will give new venues in understanding the valence quark structure of the nucleon. Finally, even though in above discussion we considered PDF in higher order approximation (nlo and nnlo), the observed PDF characteristics are largely the same also for leading order PDFs.

2.3 Main assumptions of the model

The above discussed phenomenology of valence PDFs serves us as a ground for several assumptions for the model outlined below:

• Dynamics: The main assumption of the model is that the mean field, two- and three-quark short-range correlations define the dynamics of the valence quarks in the range of $0.1 \lesssim x\lesssim 1$. In the considered model, lepton scattering from the valence quark, in leading order, proceeds according to the diagrams presented in Fig. 1. Contributions from scattering off the constituents in the residual system is neglected since their momentum fractions populate mainly at $x\lesssim 0.1$ [44]. (see also discussion in Sect. 4.2).

• Valence Quarks in the Nucleon: The model assumes an existence of almost massless three-quark valence cluster, V, in the nucleon. The cluster is compact with the transverse separation between any pairs of quarks $b_{qq} \lesssim 0.3$ Fm which indicates that the 3q system will occupy a size of up to 0.6Fm. Such a separation is consistent with both phenomenology and theoretical evaluations made in all cases in which the 3q core of the nucleon is considered (see e.g. [45,46,47,48]). In the model the valence quark system defines the baryonic number, but not necessarily the total isospin of the nucleon. It can have total isospin $I_V = {1\over 2}$ or ${3\over 2}$, each of them corresponding to different excitations or masses of the residual nucleon system. (For the lowest mass of the recoil system one expects the 3q system to have the same isospin and its projection that the considered nucleon has).

• Residual Structure: Introducing the residual structure of the nucleon with the spectrum of mass, $m_R$, represents an introduction of spectral function formalism in the description of the nucleon structure. The model assumes a certain universality of the residual structure, R, entering in all three mechanisms of generation of valence quark distribution according to Fig. 1. This universality is reflected in the fact that one can fix its main properties within one of the approaches (say within mean field model) and apply it in the calculation of 2q- and 3q-short range correlation contributions. The mass spectrum of the residual system is continuous and effectively depends on whether u- or d-valence quarks are probed. For example, a given valence quark in the proton, can originate from the 3-valence quark cluster having the total isospin and its projection the same as the proton but it can also originate from the 3-quark cluster with isospin, $I_V={1\over 2}$ and projection $I_V^3=-{1\over 2}$ or isospin, $I_V={3\over 2}$, corresponding to the higher mass of the residual system. Thus in the case of the proton one can expand the residual nucleon mass in the form:

  $$\begin{aligned} m_R(u/d)= & {} \alpha _{u/d}\cdot m_R\left( I_V={1\over 2},I^3_V= {1\over 2} \right) \nonumber \\&+ \beta _{u/d} \cdot m_R\left( I_V={1\over 2},I^3_V= -{1\over 2} \right) \nonumber \\&+ \gamma _{(u/d)} \sum \limits _{I^3_V=-{3\over 2}}^{3\over 2} m_R\left( I_V={3\over 2},I^3_V \right) + \cdots \end{aligned}$$

  (2)

  where $I^3_V$ is the isospin projection of the valence quark cluster and “$\cdots$” accounts for the possible higher mass spectrum of the residual system. Here the $\alpha$, $\beta$ and $\gamma$ factors are estimating the probabilities that u- or the d-valence quarks emerge from different isospin configurations.

  From the above representations one expects that the minimal mass in the residual system comes from the contribution of $I_V={1\over 2},I^3_V= {1\over 2}$ term that has the same isospin numbers that the proton has. In the current model we will assume the higher mass contributions to be negligible. A more elaborate model will require introducing mass distribution with different radial light front wave functions for the different isospin contributions.

  As it was mentioned above, for the cases of total isospin ($I_V={1\over 2}$, $I_V^3=-{1\over 2}$) or $I_V={3\over 2}$ the residual system will correspond to higher mass excitations. With this assumption, we observe that, for the proton having uud configuration, (per one valence quark) the u-valence quark will be accompanied with lesser residual mass than that of the d-quark. Thus one expects that for proton in general:

  $$\begin{aligned} m_R(u) < m_R(d), \end{aligned}$$

  (3)

  since for the d-valence quarks one expects relatively larger contribution from higher residual masses. For example, the second $\beta$ term contributes more in the d-valence quark case since it has udd configurations. It is interesting that such a scenario is in qualitative agreement with violation of Gottfried sum rule [49] or the “SeaQuest” result [50], since large $\beta$ factor for d-valence quark will generate an excess of ${\bar{d}}$ quarks due to $(ddu)({\bar{d}}u)$ configuration in the proton.

  In the current paper we will estimate the magnitude of $m_R$ for u- and d-quarks by fitting the calculated valence quark distributions to the empirical PDFs. In principle it could be possible to extract $m_R$ experimentally, with dedicated measurements of tagged structure functions in semi-inclusive DIS processes.

• Mean Field Model: In the considered model we assume that the mean-field dynamics (Fig. 1a) is largely responsible for the position and the peak of the valence quark structure function, h(x, t). We assume that mean-field dynamics preserves the SU(6) flavor-spin symmetry of valence quark-cluster for both the $I_V = {1\over 2}$ and ${3\over 2}$ cases. Expecting that mean the field contribution will be relevant for the region of x up to as much as 0.55–0.65 and requiring the produced final mass in DIS is $\ge 2$ GeV, we set the starting $Q^2 = 4$ GeV$^2$.

• 2q- and 3q-Correlations: Starting at $x > rsim 0.4$ the model assumes the onset of 2q-short range correlations transitioning to 3q correlations at $x > rsim 0.7$. All such correlations are taken place within the valence quark cluster. We assume that such short range correlations are generated by vector exchanges, which in the pQCD limit will correspond to the hard gluon exchanges. In the calculation of such interactions the finite range of correlations will be introduced by the momentum transfer cut-off in the gluon exchanged diagrams. Because of the selection rule, in which gluon exchanges are dominated for quarks with opposite helicities in the massless limit, we expect that with the onset of quark correlations the SU(6) symmetry of the 3q system will break down.

Using the above assumptions in developing the quark correlation model of valence quark distributions we also bring forward several experimentally verifiable predictions. These include the prediction for universality of the recoil system, R, for almost entire range of $x>0.1$ as well as the prediction for the onset of 2q- and 3q-short range quark correlations which will be distinguished by specific correlation properties between struck and recoil valence quarks. These predictions can be quantified once the model parameters are fixed by comparing the calculated valence quark distributions with PDFs extracted from DIS measurements.

It is worth mentioning that the presented approach can be applied also for the studies of QCD processes in the nuclear medium [51,52,53], in which case an additional vertex will be added to the diagrams of Fig. 1 describing the transition of nuclear target to the “nucleon-residual nucleus” intermediate state.

3 Mean-field model of valence quark distributions

Our focus now is on modeling the mean-field dynamics of valence quark interaction and calculating the valence quark distribution functions that can be compared with empirical PDFs. The mean field picture that we adopt assumes that the nucleon core consists of interacting thee-valence quark cluster embedded in the residual nucleon system consisting of sea quarks, gluons, pions and other possible hadrons. The main assumptions of the mean field model are as follows:

• The valence 3q system occupies a space in which the separation of any given quark pair in the impact parameter space is about 0.3 Fm. The interactions among these quarks is described by coupled relativistic three-dimensional harmonic oscillators, thus satisfying confinement condition.

• Valence quarks are almost massless with the invariant energy of the 3q system contributing to the nucleon mass.

• The residual system generates an external field for 3q valence system and occupies a volume less or equal to the nucleon volume.

Fig. 8

The nucleon transition into a valence and residual systems. The diagram is arranged by light-front $\tau = x^+ = z + t$ ordering. The two dashed vertical lines represent the two intermediate states where the light-front wave functions of residual and 3q systems are defined

To calculate the amplitude of virtual photon interaction with the valence quark in the above described system we notice that, in the leading order, the light-front time, $\tau$ evolution of scattering process proceeds according to the diagram in Fig. 8. Here the first vertex (from the left) describes the interaction of the 3q valence quark-cluster in the field of the residual system and the second vertex, the interaction among three valence quarks.

The scattering amplitude corresponding to the diagram of Fig. 8 can be calculated by applying effective light-front perturbation theory in which one introduces effective vertices that can be absorbed into the definition of light-cone wave functions. The present approach represents the light-front projection of the covariant amplitude which formally can be calculated within the effective Feynman diagrammatic approach [54].

3.1 General considerations: reference frame, kinematics and the cross section of the reaction

The calculations here will be done using light-cone (LC) coordinates with the 4-momenta and four-products defined as follows:

$$\begin{aligned}&k^\mu = (k^+, k^-, {\mathbf {k}}_\perp ), \ \ k^\pm = E \pm k_z, \ \ {\mathbf {k}}_\perp = (k_x, k_y), \ \ \nonumber \\&k_1 \cdot k_2 = {k_{1}^-k_{2}^+\over 2} + {k_{1}^+k_{2}^-\over 2} - \mathbf{k_{1, \perp }\cdot k_{2, \perp }}. \end{aligned}$$

(4)

We consider a reference frame where the four-momenta of the nucleon, $p^\mu _N$ and the virtual photon, $q^\mu$ are:

$$\begin{aligned}&p_N^\mu = \left( p^+_N, \frac{m_N^2}{p^+_N}, {\mathbf {0}}_\perp \right) , \ \ q^\mu = \left( 0, \frac{2p \cdot q}{p^+_N}, {\mathbf {q}}_\perp \right) , \ \ \nonumber \\&Q^2 = - q^2 = |{\mathbf {q}}_\perp |^2, \end{aligned}$$

(5)

where $m_N$ is the mass of the nucleon and we use the standard definition of Bjorken variable: $x_B = \frac{Q^2}{2 p \cdot q}$.

The important kinematical condition of the chosen reference frame is that we assume:

$$\begin{aligned} p_{N}^+ \gg m_N, k_{i}^-,k_{i,\perp }, \end{aligned}$$

(6)

where $k_i^-$ and $k_{i,\perp }$ are the “–” and transverse components of the momenta of constituents in the nucleon.

In calculating the cross section of deep inelastic inclusive scattering from the nucleon it is convenient to introduce the nucleonic tensor, $W^{\mu \nu }$, through which in the one-photon exchange approximation the differential cross section can be presented as follows:

$$\begin{aligned} {d\sigma \over dQ^2 dx} = {2\pi \alpha ^2\over Q^4 }{y^2\over Q^2}m_NL_{\mu \nu }W_N^{\mu \nu }, \end{aligned}$$

(7)

where $y = {p_N \cdot q\over p_N \cdot k}$ and the spin averaged leptonic tensor is defined as:

$$\begin{aligned} L_{\mu \nu } = 4\left[ \left( k_\mu - {q_\mu \over 2}\right) \left( k_\nu - {q_\nu \over 2}\right) \right] + Q^2\left[ -g_{\mu \nu }- {q_\mu q_\nu \over Q^2}\right] ,\nonumber \\ \end{aligned}$$

(8)

where $k_\mu$ is the four momentum of initial lepton.

If one introduces the deep-inelastic nuclear transition current as $J_N^{\mu }(p_X,s_X,p_N,s_N)$ between the initial nucleon and the final deep inelastic state, X, then the nucleonic tensor can be presented as

$$\begin{aligned} W_N^{\mu \nu }&= {1\over 4\pi m_N }\int \sum \limits _X\sum \limits _{s_X}J^{\mu ,\dagger }(p_X,s_X,p_N,s_N) \nonumber \\&\quad \times J^{\nu }(p_X,s_X,p_N,s_N) (2\pi )^4\delta ^4(q+p_N - p_X) \nonumber \\&\quad \times \delta (p_X^2-M_X^2) {d^4p_X\over (2\pi )^3}, \end{aligned}$$

(9)

where one sums over the all possible finite DIS states, X. For the most general case the nucleonic tensor for unpolarized electron scattering from an unpolarized target can be expressed through two invariant structure functions

$$\begin{aligned} W_N^{\mu \nu }&= - {F_1 \left( x_B, Q^2\right) \over m_N} \left( g^{\mu \nu } - \frac{q^\mu q^\nu }{q^2}\right) \nonumber \\&\quad + \frac{F_2\left( x_B, Q^2\right) }{m_N \left( p_N\cdot q\right) }\left( p_N^\mu - \ \frac{p_N \cdot q}{q^2} q^\mu \right) \left( p_N^\nu -\frac{p_N \cdot q}{q^2} q^\nu \right) , \end{aligned}$$

(10)

where in the reference frame of Eq. (5):

$$\begin{aligned} F_{2}(x,Q^2)&= {m_N (p_N\cdot q)\over (p_N^+)^2} W_N^{++} = {m_N Q^2\over 2x (p_N^+)^2)}W_N^{++} \nonumber \\ F_{1}(x,Q^2)&= {F_{2}(x,Q^2)\over 2x}\left( 1 + {2m_N^2x^2\over Q^2}\right) \nonumber \\&\quad - {m_N\over 2} W_N^{+-}, \end{aligned}$$

(11)

and because of the gauge invariance,

$$\begin{aligned} W_N^{+-} = {2Q\over q^-}W_N^{\perp -}. \end{aligned}$$

(12)

With these structure functions the invariant cross section of the scattering (7) can be presented as follows:

$$\begin{aligned}&{d\sigma \over dQ^2 dx} = {2\pi \alpha ^2\over Q^4}y^2\bigg ( 2F_{1}(x,Q^2) \nonumber \\&\quad + {1\over 2 x y^2} \bigg [(2-y)^2 - y^2\bigg (1 + {4m_N^2 x^2\over Q^2}\bigg )\bigg ]F_2(x,Q^2)\bigg ). \end{aligned}$$

(13)

With respect to Eq. (11) it is worth noting that in the limit of ${2m_N^2x^2\over Q^2}\ll 1$, if $W^{+-} = 0$, one arrives at Callan–Gross relation:

$$\begin{aligned} F_{2}(x,Q^2) = 2 x F_{1}(x,Q^2). \end{aligned}$$

(14)

This, together with Eq. (12) indicates that the Callan–Gross relation is fulfilled when $W^{\perp -} = 0$, which according to Eq. (9) corresponds to the absence of the transverse currents $J^\perp =0$ in the deep-inelastic scattering process in the chosen reference frame (5).

In the partonic model the structure function $F_2(x,Q^2)$ is directly related to the partonic distribution functions [32]:

$$\begin{aligned} F_2(x,Q^2) = \sum \limits _i e_i^2 x f_i(x,Q^2), \end{aligned}$$

(15)

where the $e_i$’s are the charges of interacting quarks relative to the magnitude of electron charge and $f_i(x,Q^2)$ is the PDF of the quark of type i.

From Eqs. (11) and (15) it follows that to calculate the PDFs one needs to calculate the $W^{++}$ component of the nucleonic tensor, $W^{\mu \nu }$.

3.2 Calculation of the scattering amplitude

The transition current, $J^\mu$ in Eq. (9) in our approach is identified with the scattering amplitude $A^\mu$, which corresponds to the process described by the diagram of Fig. 8. To calculate this amplitude we apply the effective light-front diagrammatic rules (summarized in Appendix A) which results in:

$$\begin{aligned}&A^\mu = \sum _{h_V, h_1} \frac{1}{k_V^+} \frac{1}{k_1^+} \bar{u}(k_1',h_1') (ie_1 \gamma ^\mu ) u(k_1,h_1) \nonumber \\&\quad \times {\prod _{i=1}^{3} \bar{u}(k_i, h_i){\varGamma }^{V \rightarrow 3q} \chi _V \over {\mathcal {D}}_2} {\chi ^\dagger _V \chi ^\dagger _R {\varGamma }^{B \rightarrow VR} u(p_N,h_N)\over {\mathcal {D}}_1} \end{aligned}$$

(16)

where ${\varGamma }^{N \rightarrow VR}$ is the effective vertex for the nucleon transition to the valence quark cluster V and residual R systems, the ${\varGamma }^{V \rightarrow 3q}$ vertex describes the transition of the state V into three individual valence quarks. (These vertices are non-perturbative objects and can not be calculated with any finite order diagrammatic approach. In our framework they will be associated with respective wave functions which will be modeled and evaluated by comparing with experiment.) All the h’s denote helicities, the $\chi _V$ and $\chi _R$ denote the respective spin functions for the valence and residual systems.² The light front energy denominators are ${\mathcal {D}}_1$ and ${\mathcal {D}}_2$, which characterize two intermediate states marked by dashed vertical lines in Fig. 8. Here we label the struck quark with $i=1$. The color indices are omitted since they are not involved in electromagnetic interaction. The sequence of $N\rightarrow VR$ followed by $V\rightarrow 3q$ transitions as presented in Fig. 8 is justified based on the assumption that the VR system is characterized by smaller internal momenta than the 3q system. The assumption of small relative momentum in the VR system allows us to approximate the first intermediate state as on-shell. For the 3q system such an approximation is not valid and in principle for the valence quark propagator one will have an on-shell and instantaneous terms due to relation: . However, since we are interested in the calculation of the $A^+$ component of the scattering amplitude which enters in $W^{++}$, the instantaneous term drops out since it couples with the extra $\gamma ^+$ factor in the $\gamma q q$ vertex and $\gamma ^+\gamma ^+ = 0$. Thus, this allows us to represent the second intermediate state by on-shell quark spinors only.

Applying the LF diagrammatic rules (Appendix A) for ${\mathcal {D}}_1$ and ${\mathcal {D}}_2$ denominators one obtains:

$$\begin{aligned} {\mathcal {D}}_1&= P^- - k_R^- - k_V^-\nonumber \\&= \frac{m_N^2}{p_N^+} - \frac{k_{R, \perp }^2 + m_R^2}{k_R^+} - \frac{k_{V, \perp }^2 + m_V^2}{k_V^+}) \nonumber \\&= \frac{1}{p_N^+} \left( m_N^2 - \frac{k_{R, \perp }^2 + m_R^2}{x_R} - \frac{k_{V, \perp }^2 + m_V^2}{x_V} \right) \nonumber \\ {\mathcal {D}}_2&= \frac{1}{k_V^+} \left( m_V^2 - \sum _{i=1}^{3} \frac{ k_{i, \perp }^2 +m_i^2 }{\beta _i} \right) , \end{aligned}$$

(17)

where $k_R$ and $k_V$ are the four momenta of the residual and valence systems respectively, $m_R$ and $m_V$ are their masses. Also, $x_j =\frac{k_j^+}{p_N^+}$ is the light-cone momentum fraction of outgoing particles, $j=1^\prime ,2,3$ and R. and, the summation is over the valence quarks, i, with $\beta _i = \frac{x_i}{x_V}$.

Substituting Eq. (17) into Eq. (16) results in:

$$\begin{aligned} A^\mu&= \sum _{h_1, h_V} \frac{1}{x_V} \frac{1}{\beta _1} \bar{u}(k_1',h_1') (ie_1 \gamma ^\mu ) u(k_1,h_1) \nonumber \\&\quad \times \frac{ \prod _{i=1}^{3} \bar{u}(k_i, h_i){\varGamma }^{V \rightarrow 3q} \chi _V}{m_V^2 - \sum _{i=1}^{3} \frac{ k_{ i, \perp }^2 +m_i^2 }{\beta _i} } \nonumber \\&\quad \times \frac{\bar{\chi _V} {\bar{\chi }}_R {\varGamma }^{B \rightarrow VR} u(p_N,h_N)}{m_N^2 - \frac{k_{V, \perp }^2 + m_V^2}{x_V} - \frac{k_{R,\perp }^2 + m_R^2}{x_R}}. \end{aligned}$$

(18)

To proceed, we introduce light-front wave functions for VR and 3q systems as follows:

$$\begin{aligned}&\psi _{VR} (x_V, {\mathbf {k}}_{R, \perp }, x_R, {\mathbf {k}}_{V, \perp }) \nonumber \\&\quad = \frac{\bar{\chi _V }\bar{\chi _R} {\varGamma }^{B \rightarrow VR} u(p_N,h_N)}{m_N^2 - \frac{k_{V, \perp }^2 + m_V^2}{x_V} - \frac{k_{R,\perp }^2 + m_R^2}{x_R}} \nonumber \\&\psi _{3q} (\{\beta _i, {\mathbf {k}}_{i, \perp }, h_i \}_{i=1}^3) = \frac{\prod _{i=1}^3 \bar{u}(k_i,h_i) {\varGamma }^{V \rightarrow 3q} \chi _V }{m_V^2 - \sum _{i=1}^{3} \frac{ k_{ i, \perp }^2 +m_i^2 }{\beta _i} }, \end{aligned}$$

(19)

where $\{\beta _i, {\mathbf {k}}_{i, \perp }, h_i \}_{i=1}^3$ denotes the LC momenta and helicities of the three valence quarks in the wave function.

Using the above definitions of LF wave functions the scattering amplitude can be presented in the following form:

$$\begin{aligned} A^\mu&= \sum _{h_1, h_V}\bar{u}(k_1',h_1') (ie_1 \gamma ^\mu ) u(k_1,h_1) \nonumber \\&\quad \times \frac{\psi _{VR}(x_V, {\mathbf {k}}_{R, \perp }, x_R, {\mathbf {k}}_{V, \perp })}{x_V} \frac{\psi _{3q}(\{\beta _i, {\mathbf {k}}_{i, \perp }, h_i \}_{i=1}^3) }{\beta _1}. \end{aligned}$$

(20)

3.3 Calculation of the valence PDF

Fig. 9

Hadronic tensor expressed as a cut diagram

The hadronic tensor, $W_N^{\mu \nu }$ of Eq. (9) for the considered scenario of scattering in which the final state consists of three outgoing valence quarks and a residual nucleon system, in the leading order, corresponds to the cut diagram of Fig. 9 for which one obtains:

$$\begin{aligned}&W_N^{\mu \nu } (x, Q^2) = \frac{1}{4\pi m_N} \sum _{ q, h_i} \int \delta (k_R^2- m_R^2) \frac{d^4 k_R}{(2\pi )^3} \nonumber \\&\quad \times \delta (k^{\prime ,2}_1- m_1^2) \frac{d^4 k^\prime _1}{(2\pi )^3} \prod _{i=2}^3 \delta (k_i^2- m_i^2) \frac{d^4 k_i}{(2\pi )^3} \nonumber \\&\quad \times (2\pi )^4 \delta ^{(4)}\left( P+q - k^\prime _1 - \sum _{i=2}^3 k_i -k_R\right) A^{\mu \dagger } A^\nu , \end{aligned}$$

(21)

where the sum $\sum \nolimits _{ q, h_i}$ is over all the spins/flavors of the outgoing valence quarks and the residual system. The Lorentz invariant phase space of all final particles can be presented as follows:

$$\begin{aligned} \delta (k_j^2-m_j^2) d^4k&= \frac{1}{2} \delta (k_j^+k_j^- - k_{j,\perp }^2 -m_j^2) dk_j^- dk_j^+ d^2{\mathbf {k}}_{j, \perp } \nonumber \\&= \frac{ dx_j d^2{\mathbf {k}}_{j, \perp }}{2x_j} |_{k_j^- = \frac{k_{j,\perp }^2 +m_j^2}{x_jp_N^+}}, \end{aligned}$$

(22)

To evaluate the energy–momentum conservation condition, first we express the $\delta ^4()$ function through the LC components:

$$\begin{aligned}&\delta ^{\left( 4\right) }\left( p_N+q - k^\prime _1 - \sum _{i=2}^3 k_i -k_R\right) \nonumber \\&\quad = 2 \delta \left( p_N^+ + q^+ - k^{\prime ,+}_1 - \sum _{i=2}^3 k_i^+ - k_R^+\right) \nonumber \\&\qquad \times \delta \left( p_N^- + q^- - k^{\prime ,-}_1 - \sum _{i=2}^3 k_i^- - k_R^-\right) \nonumber \\&\qquad \times \delta ^{\left( 2\right) }\left( {\mathbf {p}}_{N,\perp } + {\mathbf {q}}_\perp - k^{\prime }_{1,\perp } - \sum _{i=2}^3 {\mathbf {k}}_{i, \perp } - {\mathbf {k}}_{R, \perp }\right) . \end{aligned}$$

(23)

Using the reference frame of Eq. (5) in which $q^+=0$, one observes that $k^+_1 = k^{\prime ,+}_1$ and using the definition of LC momentum fractions $x_j =\frac{k_j^+}{p_N^+}$ one obtains:

$$\begin{aligned}&\delta \left( p_N^+ + q^+ - k^{\prime +} _1- \sum _{i=2}^3 k_i^+ - k_R^+\right) \nonumber \\&\quad = \frac{1}{p_N^+} \delta \left( 1 - \sum \limits _{i=1}^3x_i- x_R\right) . \end{aligned}$$

(24)

For the “–” component we use $q^- = {2p_N \cdot q\over p_N^+}$ (Eq. (5)), the on-mass shell condition for outgoing particles (see Eq. (22)), and conditions of high energy scattering of Eq. (6). In the large $Q^2$ limit for which $k^\prime _{1\, \perp } = k_{1,\perp } + q_\perp \gg k_{1,\perp }$ and $k^{\prime 2}_{1, \perp } \approx Q^2$ one obtains:

$$\begin{aligned}&\delta \left( p_N^- + q^- - k^{\prime -}_{1} - \sum _{i=2}^3 k_i^- - k_R^-\right) \nonumber \\&\quad = \delta \bigg ( \frac{m_N^2}{p_N^+} + \frac{2 p_N \cdot q}{p_N^+} - \frac{k^{\prime 2}_{1,\perp } +m^2}{x^\prime _{1}p_N^+} \nonumber \\&\qquad - \sum _{i=2}^3 \frac{k_{i,\perp }^2 +m^2}{x_{i}p_N^+} -\frac{k_{R,\perp }^2 +m_R^2}{x_{R}p_N^+}\bigg ) \nonumber \\&\quad \approx \delta \left( \frac{2 p_N \cdot q}{p_N^+}- \frac{Q^2}{x_{1}'p_N^+}\right) \nonumber \\&\quad = \frac{x_{1}p_N^+}{2 p_N \cdot q} \delta (x_1 - x_B) |_{x_B = \frac{Q^2}{2 p_N \cdot q}}. \end{aligned}$$

(25)

Furthermore, introducing the initial transverse momentum of the struck quark ${\mathbf {k}}_{1\perp } = {\mathbf {k}}^\prime _{1,\perp } - {\mathbf {q}}_\perp$ for the conservation of transverse momentum one obtains:

$$\begin{aligned}&\delta ^{(2)}\left( {\mathbf {p}}_{N,\perp } + {\mathbf {q}}_\perp - \mathbf{k}^{\prime }_{1,\perp } - \sum _{i=2}^3 {\mathbf {k}}_{i, \perp } - {\mathbf {k}}_{R, \perp }\right) \nonumber \\&\quad = \delta ^{(2)}\left( \sum _{i=1}^3 {\mathbf {k}}_{i, \perp } + {\mathbf {k}}_{R, \perp }\right) . \end{aligned}$$

(26)

The results of Eqs. (24), (25) and (26) can be summarized as:

$$\begin{aligned}&\delta ^{\left( 4\right) }\left( p_N+q - k^\prime _1 - \sum _{i=2}^3 k_i -k_R\right) \nonumber \\&\quad = \frac{x_1}{ p_N \cdot q} \delta \left( 1 - \sum \limits _{i=1}^3x_i- x_R\right) \delta \left( x_1 - x_B\right) \nonumber \\&\qquad \times \delta ^{\left( 2\right) }\left( \sum _{i=1}^3 {\mathbf {k}}_{i, \perp } + {\mathbf {k}}_{R, \perp }\right) . \end{aligned}$$

(27)

Substituting Eq. (27) in Eq. (21) one obtains:

$$\begin{aligned}&W_N^{\mu \nu } (x, Q^2) \nonumber \\&\quad = \frac{1}{4 m_N} \sum _{h_i, q} \int [dx][d^2 {\mathbf {k}}_\perp ] \frac{x_1}{ p_N \cdot q} \delta (x_1 - x_B) A^{\mu \dagger } A^\nu , \end{aligned}$$

(28)

where the integration factors are defined as:

$$\begin{aligned}&[dx] =\delta \left( 1- \sum _{i=1}^3 x_i - x_R\right) \frac{dx_R}{x_R}\prod _{i=1}^{3} \frac{dx_i }{x_i} \end{aligned}$$

(29)

$$\begin{aligned}&[d^2 {\mathbf {k}}_\perp ] = 16 \pi ^3 \delta ^{\left( 2\right) }\left( \sum _{i=1}^3 {\mathbf {k}}_{i,\perp } + {\mathbf {k}}_{R, \perp }\right) \frac{d^2 {\mathbf {k}}_{R, \perp }}{16 \pi ^3 }\prod _{i=1}^{3} \frac{d^2 {\mathbf {k}}_{i,\perp }}{16 \pi ^3}.\nonumber \\ \end{aligned}$$

(30)

Substituting now the “+” component of the scattering amplitude from Eq. (20) to Eq. (28) and using Eq. (11) one obtains for the $F_2$ structure function the following expression:

$$\begin{aligned}&F_2(x_B) = \sum _q \sum _{h_i} \int [dx] [d^2{\mathbf {k}}_\perp ] e_q^2 x_1 \delta (x_1 - x_B)\nonumber \\&\quad \times |\psi _{3q}(\{\beta _i, {\mathbf {k}}_{i, \perp }, h_i \}_{i=1}^3) |^2 |\psi _V (x_V, {\mathbf {k}}_{R, \perp }, x_R, {\mathbf {k}}_{V, \perp })|^2, \end{aligned}$$

(31)

which together with Eq. (15) allows us to obtain the expression for the valence PDF in the form:

$$\begin{aligned}&f_q(x_B)= \sum _{h_i} \int [dx] [d^2{\mathbf {k}}_\perp ] \delta (x_1 - x_B) \nonumber \\&\quad \times |\psi _{3q}(\{\beta _i, {\mathbf {k}}_{i, \perp }, h_i \}_{i=1}^3) |^2 |\psi _V (x_V, {\mathbf {k}}_{R, \perp }, x_R, {\mathbf {k}}_{V, \perp })|^2. \end{aligned}$$

(32)

3.4 Modeling the wave functions

3.4.1 An approach for modeling LF wave functions

There have been many approaches in modeling light-front valence quark wave functions of hadrons (see e.g. [3, 34, 55]). Quantum mechanically, for a bound system the expansion of the potential close to the stability point results in a harmonic oscillator (HO) type interaction potential, thus justifying the use of the HO eigenfunctions as a basis for modeling nonperturbative wave function of hadrons. From a practical point of view, the advantage of using HO wave functions was that they naturally contain the effect of confinement and the Gaussian form was convenient for analytic calculations. The HO approach was used both in constituent (e.g. Ref.  [1]) and current (e.g. Ref.  [34]) quark descriptions of hadrons. For the case of current quark description, one important issue was a proper relativistic generalization of the wave functions. For two valence quark systems such as pions, one straightforward step in the relativistic generalization was introducing the relative momentum of two valence quarks in the center of mass frame according to $k^2 = {m_q^2 + k^2_{\perp }\over x(1-x)} - m_q^2$, which allowed to represent the wave function in the boost-invariant form along the momentum transfer direction $\mathbf{q}$ [3, 34].

In our approach we again use the HO basis for the LF wave functions. In the relativistic generalization of such wave functions we introduce the relative momentum similar as it was defined above for two-valence quark system. However, we also introduce the phase space factor for the residual (spectator) system that allows us to recover the non-relativistic normalization of the wave function in the lab frame where the wave function represents the eigenstate of the HO Hamiltonian. This phase factor appears also if we relate Weinberg (or Bethe–Salpeter) type equations in the nonrelativistic limit to the Lippmann–Schwinger equation for bound systems in the lab frame, where the potential of the interaction is defined (see e.g. [54]). For example, in considering a generic two-body bound system, the HO based LF wave function is represented in the following form:

$$\begin{aligned} \psi _2(\alpha _1, {\mathbf {p}}_{1\perp },\alpha _2, {\mathbf {p}}_{2\perp }) = \sqrt{2(2\pi )^3M_T} {\varPsi }_0(k_{cm}) \sqrt{\alpha _2}, \end{aligned}$$

(33)

where ${\varPsi }_0$ is the ground state wave function of the non-relativistic HO Hamiltonian in the lab frame of the target with mass $M_T$, normalized conventionally as:

$$\begin{aligned} \int |{\varPsi }_0(p)|^2 d^3p = 1. \end{aligned}$$

(34)

In Eq. (33) the residual (spectator) particle “2” is treated as on-shell with $\alpha _1 = 1-\alpha _2$, $\mathbf{p}_{1, \perp } = - \mathbf{p}_{2, \perp }$, and $\alpha _2 = {E^{lab}_{2} + p^{lab}_{2z}\over M_T} = {E^{cm}_{2} + k^{cm}_{z} \over E_{2,cm} + E_{1,cm}}$, where

$$\begin{aligned} k_{cm}^2 = \frac{(s-(m_1-m_2)^2)(s-(m_1+m_2)^2)}{4s}, \end{aligned}$$

(35)

with

$$\begin{aligned} s = (p_1 + p_2)^2 = \frac{p_{1,\perp }^2 +m^2}{\alpha _1} + \frac{p_{2,\perp }^2 +m^2}{\alpha _2}. \end{aligned}$$

(36)

The extra factor of $\sqrt{\alpha _2}$ in Eq. (33) is to account for the phase space of the residual nucleon that allows to relate the relativistic normalization of the LF wave function to the nonrelativistic normalization (34) in the lab frame. Namely starting with the relativistic normalization:

$$\begin{aligned} \int |\psi _2(\alpha _1, {\mathbf {p}}_{1, \perp },\alpha _2, {\mathbf {p}}_{2, \perp })|^2 \frac{d\alpha _2}{\alpha _2}\frac{ d^2 p_{2,\perp }}{16\pi ^3}= 1, \end{aligned}$$

(37)

and using Eq. (33) one obtains:

$$\begin{aligned}&M_T16\pi ^3 \int |{\varPsi }_0(k_{cm})|^2 \alpha _2 \frac{d\alpha _2}{\alpha _2}\frac{ d^2 p_{2,\perp }}{16\pi ^3} \nonumber \\&\quad \approx M_T 16\pi ^3 \int |{\varPsi }_0(k_{cm})|^2 {d p_{2,z}\over M_T} \frac{ d^2 p_{2,\perp }}{16\pi ^3} \nonumber \\&\quad = \int |{\varPsi }(p_2)|^2 d^3p_{2} = 1, \end{aligned}$$

(38)

where in the above derivation we used the nonrelativistic limit for $k_{cm}\rightarrow p_2$ and $\alpha _2 = \frac{E_2 + p_{2,z}}{M_T}\approx \frac{m_2 + p_{2,z}}{M_T}$ with $\rightarrow d\alpha _2 = \frac{dp_{2,z}}{M_T}$.

It is worth mentioning that the above described prescription is not unique to HO wave functions and can be applied to any quantum mechanical wave function.

3.4.2 Wave function of three valence quark system

Fig. 10

Three valence quarks coupled in pairwise HO interaction

We model the three valence quark system as mutually coupled oscillators (Fig. 10) for which the ground state wave function within the above described prescription can be presented in the following form:

$$\begin{aligned}&\psi _{3q} ( \{x_i, {\mathbf {k}}_{i,\perp }\}_{i=1}^{3}) \nonumber \\&\quad = 16\pi ^3 m_N A_V \exp { \left[ - \frac{B_V}{2}( k_{12,cm}^2 + k_{23,cm}^2 + k_{31,cm}^2) \right] } \nonumber \\&\qquad \times \sqrt{x_2 x_3}, \end{aligned}$$

(39)

where $A_V$ and $B_V$ are parameters and $x_i, k_{i,\perp }$, ($i\ne j=1,2,3$) are LC momentum fractions and transverse momenta of each valence quark in the reference frame defined in Eq. (5). The $k_{ij, cm}^2$s, ($i\ne j=1,2,3$) in the exponent of the wave function represent the relative three momenta in the CM system of i, j pairs defined as follows:

$$\begin{aligned} k_{ij, cm}^2 = \frac{(s_{ij}-(m_i-m_j)^2) (s_{ij} - (m_i + m_j)^2)}{4s_{ij}}, \end{aligned}$$

(40)

where the invariant energy of the i, j pair is:

$$\begin{aligned} s_{ij}&= (k_i+k_j)^+ (k_i + k_j)^- - ({\mathbf {k}}_{i,\perp } + {\mathbf {k}}_{j,\perp })^2 \nonumber \\&= (x_i + x_j) \left( \frac{k_{i, \perp }^2 + m_i^2}{x_i} + \frac{k_{j, \perp }^2 + m_j^2}{x_j}\right) \nonumber \\&\quad -({\mathbf {k}}_{i,\perp } + {\mathbf {k}}_{j,\perp })^2. \end{aligned}$$

(41)

Since the momentum fractions defined as $x_i = {k^+_{i}\over p^+_N}$, one can introduce total momentum fraction carried by three-valence quark system as:

$$\begin{aligned} x_V = x_1 + x_2 + x_3. \end{aligned}$$

(42)

Next we introduce the valence quark transverse momenta defined with respect to the total momentum of 3q system $\mathbf{k}_V$:

$$\begin{aligned} {\tilde{\mathbf {k}}}_{i, \perp } = {\mathbf {k}}_{i, \perp } - \frac{x_i}{x_V} {\mathbf {k}}_{V, \perp }, \quad (\text{ i } = 1,2,3), \end{aligned}$$

(43)

where

$$\begin{aligned} {\tilde{k}}_{1,\perp } + {\tilde{k}}_{2,\perp } + {\tilde{k}}_{2,\perp } = 0. \end{aligned}$$

(44)

Using Eqs. (43) and (44) and assuming same masses for valence quarks one obtains:

$$\begin{aligned} s_{12} + s_{23} + s_{31} = \sum _{i=1}^3 x_V\frac{{\tilde{k}}_{i, \perp }^2 + m^2}{x_i} + 3m^2, \end{aligned}$$

(45)

and

$$\begin{aligned} k_{12,cm}^2 + k_{23,cm}^2 + k_{31,cm}^2= \frac{1}{4} \left( \sum _{i=1}^3 x_V\frac{{\tilde{k}}_{i, \perp }^2 + m^2}{x_i} - 9m^2 \right) .\nonumber \\ \end{aligned}$$

(46)

Inserting the above relation into Eq. (39) the wave function of 3q-state can be expressed as follows:

$$\begin{aligned}&\psi _{3q} ( \{ x_i, {\mathbf {k}}_{i, \perp } \}_{i=1}^3) \nonumber \\&\quad = 16\pi ^3 m_N A_V \exp \left[ -\frac{B_V}{8} \left( \sum _{i=1}^3 x_V\frac{{\tilde{k}}_{i, \perp }^2 + m^2}{x_i} - 9m^2 \right) \right] \nonumber \\&\qquad \times \sqrt{x_2 x_3}. \end{aligned}$$

(47)

3.4.3 Wave function of recoil system

Using the same approach (Sect. 3.4.1) we model the wave function of the recoil system in the following form:

$$\begin{aligned} \psi _R (x_R, {\mathbf {p}}_{R, \perp }) = \sqrt{16\pi ^3 m_N} A_R e^{- \frac{B_R}{2} p_R^2} \sqrt{x_R} \end{aligned}$$

(48)

where $x_R$ is the light-cone momentum fraction of the nucleon carried by the recoil system and $p_R^2 = k_{R, \perp }^2 + p_{R, z}^2$ is its three momentum in the CM of the nucleon. As a first approximation that allows to perform calculations analytically, we use a non-relativistic kinematics for the recoil system. In this case in the CM frame of the nucleon one can approximate $x_R \approx {m_R + p_{R,z}\over m_N}$ and substitute

$$\begin{aligned} p_{R,z} = (x_{R}m_N- m_R). \end{aligned}$$

(49)

In numerical evaluations we will treat $B_R$ and $m_R$ as free parameters characterizing the recoil system and evaluate $A_R$ from the normalization condition similar to Eq. (37).

3.5 Completing the calculation of valence PDF

Substituting the wave functions of 3q (Eq. (47)) and the recoil (Eq. (48)) systems in Eq. (32) for the valence quark partonic distribution one obtains:

$$\begin{aligned}&f_q(x_B)= \left( 16 \pi ^3m_N\right) ^3 |A|^2 \int ^{Q^2} 16 \pi ^3\delta (x_1-x_B)\nonumber \\&\delta \left( 1-\sum \limits _{1}^3x_i - x_R\right) \delta ^2\left( \sum \limits _{1}^3k_{i,\perp }+k_{R,\perp }\right) \nonumber \\&\quad \times \exp \bigg [-{B_V\over 4}\sum _{i=1}^3 x_V \frac{{\tilde{k}}_{i,\perp }^2 +m_i^2}{x_i} \nonumber \\&\quad - B_R \left( (m_N x_R - m_R)^2 + k_{R,\perp }^2 \right) \bigg ] x_2x_3x_R\nonumber \\&\quad \times \prod \limits _{i=1}^3 {dx_i\over x_i} {d^2k_{i,\perp }\over 16 \pi ^3}{dx_R\over x_R} {d^2k_{R,\perp }\over 16 \pi ^3}, \end{aligned}$$

(50)

where $A = A_R A_V e^{{9\over 4} B_Vm_q^2}$. Note that the transverse momenta $\mathbf{k}_{i,\perp }$ and ${\tilde{\mathbf{k}}}_{i,\perp }$ are defined in the CM of nucleon and 3q-system respectively and related to each other according to Eq. (43). Using Eq. (43) and relation $\mathbf{k}_{V,\perp } = - \mathbf{k}_{R,\perp }$ and $d^2 {\mathbf {k}}_{i, \perp } = d^2 {\tilde{\mathbf {k}}}_{i, \perp }$ allows us to factorize integrations by ${\tilde{k}}_{i,\perp }$ and $k_{R,\perp }$ resulting in:

$$\begin{aligned}&f_q(x_B, Q^2) = m_N^3 |A|^2 \int \prod \limits _{i=1}^3 {dx_i\over x_i} \delta (x_B-x_1) \nonumber \\&\quad \times \exp \left[ -\frac{B_V}{4} \sum _{i=1}^3 x_V\frac{ m_i^2}{x_i} - B_R m_N^2\left( x_R - \frac{m_R}{m_N}\right) ^2 \right] \nonumber \\&\quad \times x_2 x_3 \int ^{Q^2} \prod \limits _{i=1}^3 d^2{\tilde{k}}_{i,\perp } \delta ^2\left( \sum \limits _{i=1}^3 {\tilde{\mathbf{k}}}_{i,\perp } \right) \nonumber \\&\quad \times \exp { \left[ -\frac{B_V}{4} \sum _{i=1}^3 x_V \frac{{\tilde{k}}_{i,\perp }^2}{x_i}\right] } \int ^{Q^2} d^2k_{R,\perp } \exp {\left[ - B_R k_{R,\perp }^2 \right] }. \end{aligned}$$

(51)

Furthermore we evaluate the following integrals in the above equation (see Appendix B):

$$\begin{aligned}&\int ^{Q^2_R} \exp {\left[ - B_R k_{R,\perp }^2\right] } d^2p_{R,\perp } = {\pi \over B_R}(1-e^{-B_R Q^2_{R,Max}});\nonumber \\&\int ^{Q^2} \prod \limits _{i=1}^3 d^2{\tilde{k}}_{i,\perp } \delta ^2\left( \sum \limits _{i=1}^3 {\tilde{\mathbf{k}}}_{i,\perp } \right) \exp { \left[ -\frac{B_V}{4} \sum _{i=1}^3 x_V \frac{{\tilde{k}}_{i,\perp }^2}{x_i}\right] } \nonumber \\&\quad = {16\pi ^2\over B_V^2} {x_1x_2x_3\over x_V^3}\nonumber \\&\qquad \times (1-e^{-a_{cm}Q^2_{cm,Max}})(1-e^{-a_{rel}Q^2_{rel,Max}}), \end{aligned}$$

(52)

where $a_{cm} = {B_V x_V\over 4} {x_V\over x_3(x_1+x_2)}$ and $a_{rel} = {B_V x_V\over 4} {x_1+x_2\over x_1x_2}$. Here $Q^2_{R,Max}$ is the maximal transverse momentum square of the recoil system. Also in the derivation (see Appendix B) we defined ${\tilde{\mathbf{k}}}_{12,\perp }^{cm} = {\tilde{\mathbf{k}}}_{1,\perp } + {\tilde{\mathbf{k}}}_{2,\perp }$ and ${\tilde{\mathbf{k}}}_{12,\perp }^{rel} = {x_2{\tilde{\mathbf{k}}}_{1,\perp } - x_1 {\tilde{\mathbf{k}}}_{2,\perp }\over x_1 + x_2}$ and $Q^2_{cm,Max}$ and $Q^2_{rel,Max}$ represent the maximal transverse momentum squares of the center of mass and relative transverse momenta of “1” and “2” quark system. Inserting the expressions of Eq. (52) into Eq. (51) and taking the integral over $d x_1$ using $\delta (x_B-x_1)$ (see Appendix B for details) one arrives at:

$$\begin{aligned}&f_q(x_B, Q^2) = {\mathcal {N}} \displaystyle \int _0^{1-x_B}dx_2 \displaystyle \int _0^{1-x_B - x_2} dx_3 \nonumber \\&\quad \times \exp \left[ -\frac{B_Vx_V}{4} \sum _{i=1}^3 \frac{ m_i^2}{x_i} - B_R M_N^2\left( x_V - \left( 1- \frac{M_R}{M_N}\right) \right) ^2 \right] \nonumber \\&\quad \times \frac{x_2 x_3}{x_V^3} \left( 1-e^{-a_{cm}Q_{cm}^{max 2}}\right) \left( 1-e^{-a_{rel}Q_{rel}^{max 2}}\right) \nonumber \\&\quad \times \left( 1-e^{-B_R Q^2}\right) , \end{aligned}$$

(53)

where $x_1= x_B$ and $x_V= x_B + x_2 + x_3$ and the normalization constant

$$\begin{aligned} {\mathcal {N}} = {16\pi ^3 A_V^2 A_R^2 m_N^3\over B_R B_V^2} e^{{9\over 4} B_V m_q^2}. \end{aligned}$$

(54)

The above equation simplifies further when $Q^2$ is large so that terms with it in the exponential are negligible in Eq. (53):

$$\begin{aligned}&f_q(x_B, Q^2) = {\mathcal {N}} \int _0^{1-x_B}dx_2 \int _0^{1-x_B - x_2} dx_3 \frac{x_2 x_3}{x_V^3}\nonumber \\&\quad \times \exp \left[ -\frac{B_Vx_V}{4} \sum _{i=1}^3 \frac{ m_i^2}{x_i} - B_R M_N^2\left( x_V - \left( 1- \frac{m_R}{m_N}\right) \right) ^2 \right] . \end{aligned}$$

(55)

We further simplify this expression by considering the massless limit of valence quarks and changing the integration variable from $x_3$ to $x_V$. This allows us to evaluate the $dx_2$ integration analytically resulting in the final expression for the valence quark distribution in the following form:

$$\begin{aligned}&f_q(x_B, Q^2) = {{{{\mathcal {N}}}}\over 6} \int \limits _{x_B}^1 dx_V {(x_V-x_B)^3\over x_V^3}\nonumber \\&\quad \times \exp \left[ -B_R m_N^2\left( x_V - \left( 1- \frac{m_R}{m_N}\right) \right) ^2 \right] . \end{aligned}$$

(56)

Note that as it follows from the above equation, in the massless limit, the parameter $B_V$ enters in the structure function only through the normalization factor ${{\mathcal {N}}}$.

3.6 Qualitative features of the model

One important feature of the model can be observed if we consider Eq. (56) at $x_B\sim x_p$, when $x_V\sim 1$. In this case the qualitative behavior of $f_q(x_B, Q^2)$ can be obtained by evaluating Eq. (56) at the maximum of the exponent; $x_V = 1- \frac{m_R}{m_N}$. This results in $f_q(x_B,Q^2)\sim (1- x_B - {m_R\over m_N})^3$ and for the valence quark structure function one obtains:

$$\begin{aligned} h(x_b,t) = x_B f_q(x_B,Q^2)\sim x_B\left( 1- x_B - {m_R\over m_N}\right) ^3. \end{aligned}$$

(57)

From the above relation it follows that h(x, t) function peaks at

$$\begin{aligned} x_{p} \approx {1\over 4}\left( 1-{m_R\over m_N}\right) . \end{aligned}$$

(58)

Thus one observes that within the considered model the peak of the h(x, t) distribution for valence quarks (see Fig. 3b) is related to the mass of the residual system in the nucleon.

Now if we use the characteristic value for the peak position $x_{p} \simeq 0.2$ for valence quarks at starting $Q^2=4$ GeV$^2$, from Eq. (58) one observes that it will correspond to $m_R\approx m_\pi$. Thus the present model mimics a pion-cloud like picture for the residual system.

Equation (58) results in an interesting parameter free observation that for any realistic PDF for the all range of $Q^2$, where they are defined, the position of the peak of valence quark distribution should not exceed ${1\over 4}$: that is $x_p\lesssim {1\over 4}$.³ Since QCD evolution shifts $x_p$ towards the lower values (Fig. 4) it is interesting to check ${1\over 4}$ limit for starting $Q_0$ from which input PDFs evolved to higher $Q^2$ values. For this, we considered all current ‘standard’ PDFs for which starting $Q_0 > rsim 1$ GeV$^2$. Here by ‘standard’ we refer to all PDFs for which the input PDF at starting $Q_0$ corresponds to a physical nonperturbative parton distribution. (This to be distinguished from ‘dynamical’ PDFs in which $Q_0$ (which can be less than 1 GeV$^2$) is introduced as a free parameter with input PDFs fixed by the best fit to the experimental data. As a result PDFs at starting $Q_0$ do not correspond to realistic partonic distributions. For details see Ref.  [57].)

Fig. 11

The position of the peak of $x q_v$ distribution for valence u-quarks at different starting $Q_0$. References  [36,37,38, 41,42,43, 58, 59]

Considering the above described ‘standard’ PDFs in Fig. 11 we observe that all the different parameterizations produce peak position for valence u-quarks less than ${1\over 4}$ for the commonly used charm quark mass scale $Q_0 = m_C = 1.3$ GeV, which is a natural starting point if one introduces charm component through QCD evolution [60, 61]. It is also intriguing that those ‘standard’ PDFs which have smaller starting scale at $Q_0= 1$ GeV indicate an $x_p$ approaching to ${1\over 4}$. Note that for all these sets $x_p^{d_V} <x_p^{u_V}$, therefore the same limits holds also for the down valence quarks.

It is interesting that in Ref.  [10] it was observed that within diquark model the magnitude of $x_p$ is related to the mass of the recoil diquark in the form:

$$\begin{aligned} x_p \approx \left( 1- {m_d\over m_N}\right) , \end{aligned}$$

(59)

where $m_d$ is the mass of the diquark. This and Eq. (59) however are different, reflecting different dynamics of valence quark generation in the nucleon.

It is worth mentioning that with an increase of $Q^2$, because of the radiation of valence quarks the magnitude of $x_V$ diminishes and due to the ${(x_V-x_B)^3\over x_V^3}$ factor in Eq. (56) the $x_p~-~m_R$ relation become more complicated than Eq. (58).

Next we evaluate valence quark PDFs at $x_B\rightarrow 1$ limit. For this we substitute $x_B = 1 - \epsilon$ in Eq. (56) and in the $\epsilon \rightarrow 0$ limit evaluate the integral which results in

$$\begin{aligned} f_q(x_B, Q^2)\mid _{x_B\rightarrow 1} = {{{{\mathcal {N}}}}\over 24}e^{ -B_R m_R^2} (1-x_B)^4. \end{aligned}$$

(60)

This result indicates that the considered mean field model predicts a fall off of the valence quark PDF faster than the one predicted within perturbative QCD [34]; $(1-x_B)^3$. Thus the observation of the $(1-x_B)^4$ behavior at $x_B\rightarrow 1$ limit will indicate the dominance of the mean field dynamics.

4 Numerical estimates of valence quark distributions

4.1 Strategy of choosing the parameters of the model

We use Eq. (56) as a baseline formula for fitting to the phenomenological PDFs evaluated in leading order. In general the model has five parameters $A_V$ $B_V$, $A_R$, $B_R$ and $m_R$.

For the valence quarks we assume that characteristic separations in the 3q system in the impact parameter space is $\langle b^2_{i,j}\rangle \sim$ (0.3Fm)$^2$. This allows us, based on Eq. (47), to evaluate $B_V = 4 \langle b^2_{i,j}\rangle {x_i\over x_V} \approx {4\over 3} \langle b^2_{i,j}\rangle \approx 3.08$ GeV$^{-2}$.

For the recoil system we relate $A_R = \left( {B_R\over \pi }\right) ^{3\over 4}$. Finally, the remaining parameter $A_V$ is fixed through the normalization factor, ${{{\mathcal {N}}}}$ and Eq. (54), yielding:

$$\begin{aligned} A_V =\sqrt{B_R{{{\mathcal {N}}}}\over 16 \pi ^3 m_N ^3} {B_V\over A_R}. \end{aligned}$$

(61)

In the fitting procedure the parameters ${{{\mathcal {N}}}}$, $m_R$ and $B_R$ are evaluated by the fitting Eq. (56) using the maximum likelihood method with the constraint to match the height and position in leading order approximation.

Fig. 12

The h(x) distribution at $Q^2=4$ GeV$^2$. Left panel for valence d-quark distribution. Right panel for valence u-quark distributions

As a starting $Q^2$ we choose $Q^2=4$ GeV$^2$ such that it provides large enough produced mass $W\ge 2$ GeV in the range of up to $x \sim 0.65$ for which we expect that mean-field dynamics have non-negligible contribution. This $Q^2$ is also large enough that it resolves the 3q valence cluster in the nucleon, for which one requires $Q^2\gg {1\over B_V}$. In Fig. 12 we show the current status of the valence PDFs in LO for d- and u-quarks evaluated by using modern PDF parameterizations [36, 38, 41,42,43] at $Q^2=4$ GeV$^2$. As the figure shows the height and the peak-position of the distributions vary between different parameterization on the level of 15–20%. Thus by fitting to these distributions we evaluate the range of parameters as they change from one to other parameterizations. We fit using maximum likelihood approximation without estimating the errors of parameters for specific PDF parameterization. In such way our estimates can be considered qualitative (see also discussion in Sect. 5).

4.2 Fitting results and estimation of the magnitude of high momentum component of the valence quark distributions

The results for parameters ${{\mathcal {N}}}$, $B_R$ and $m_R$ for considered PDF parameterizations are given in Table 1. As these parameters show despite their variations, there are common features shared by different parameterizations. The one is that the residual mass for the case of valence u-quarks ($0.04 \le m_R^u \le 0.07$ GeV) is less than the corresponding mass for the d-quark ($0.17 \le m_R^d \le 0.32$ GeV). Additionally, $m_R^u\le m_{pion}$ and $m_R^d > m_{pion}$. These evaluations predict that in the case of scattering from valence d-quark in the proton there is a larger probability to detect recoil pions than in the case of scattering from the u-quark. Also, the maximum of the integrand in Eq. (56) is at $x_V = 1 - \frac{m_R}{m_N}$ and thus we expect $x_R =\frac{m_R}{m_N} \sim 0.05 - 0.34$. If we assume that the residual system is a pion, then the quark in the pion should have momentum fraction $x = \frac{x_R}{2} \sim 0.025 - 0.17$. This is consistent with the expected region where pion contributions via the Sullivan process should be significant in DIS probing of the nucleon [44]. Our systematic underestimation of the lower x region, as seen in Fig. 13, could be in part from neglecting the contribution from the scattering off the quarks in residual system.

Table 1 Fitting parameters for valence d- and u-quarks

Another characteristics of the parameters is that the slope factor for the residual system in the case of scattering from valence u-quarks ($42 \le B_R^u \le 65$ GeV$^{-2}$) is systematically larger than for the case of d-quarks ($29 \le B_R^d \le 50$ GeV$^{-2}$) (not included is the CJ15 parametrization for the d-quarks which yields $B_R^d \le 6$ GeV$^{-2}$.). This indicates a more compact state for the case of scattering from d quarks which can be checked by measuring the attenuation of the residual system in the nucleus depending whether scattering took place on valence d- or u-quark.

In Fig. 13 we present a comparison of the model with the CT14nnlo PDFs at $Q^2=4$ GeV$^2$.⁴ As the figure shows, a better description is achieved for the valence d-quark than for the u-quark distribution. This indicates that larger high momentum component is needed for valence u-quarks than for the d-quark.

To evaluate the expected contribution from the high momentum component of valence quark distribution, we used the fitting parameters to evaluate normalization and momentum sum rules for both valence d- and u-quarks. For the normalization one obtains for the d-quark $N_d = 0.64$ and for the u-quark ${\mathcal {N}}_u = 1.37$. It is interesting that these normalizations result in a ${{\mathcal {N}}_d\over {\mathcal {N}}_u}\sim 0.5$. These results also indicates that one expects that Regge mechanism at $x<x_p$ together with qq-correlations to contribute $\sim$ 36% and $\sim$32% of total normalizations for d- and u-quarks respectively.

Fig. 13

Comparison of the results of the fitting of mean-field model to the $x q_V(x)$ distribution for d- (a) and u- (b) quarks. Data points corresponds to the CT14nnlo parameterization for central PDF sets and the shaded areas shows the uncertainties of PDFs estimated through the PDF error sets (for details see Ref.  [38])

The evaluation of the momentum sum rule ($P_q = \int x q_V(x) dx$) yields $P_{d} = 0.095$ and $P_u = 0.24$ which should be compared with $P_{d} = 0.1$ and $P_{u}= 0.268$ obtained from the CT14llo distribution. Expecting that the Regge mechanism, which dominates at $x<x_p$, has a smaller contribution to the momentum sum rule the above estimates allow us to evaluate better the contribution from qq correlations. As the mentioned numbers show, the expected contribution from qq correlations in the momentum sum rule for the d-quark is $\sim$ 5% and for the u-quark is $\sim$ 10%. These evaluations however can be somewhat of an underestimation, since in fitting the model to the PDFs we assumed that all the strength of PDF at $x=x_p$ is due to the mean-field mechanism. The larger contribution for quark-correlation mechanism may require a renormalization of the mean field distribution. This problem can be addressed after modeling high-x tail of valence PDFs and combining it with the current mean-field model. Overall our current result indicates an expected larger contribution of high momentum tail to the valence u-quark distribution compared to that of d-quarks.

Another interesting observable is the ratio of the d- to u-quarks at $x\rightarrow 1$ limit. As it follows from Eq. (60) the model predicts

$$\begin{aligned} {f_d(x)\over f_u(x)} \mid _{x_B\rightarrow 1} = {{\mathcal {N}}_d\over {\mathcal {N}}_u}e^{\left( B_R^u(m_R^u)^2 -B_R^d (m_R^d)^2\right) }. \end{aligned}$$

(62)

The interesting feature of the model is that the properties of valence PDFs (the mass of the residual system, the slope factor and normalization constant) at the vicinity of the peak define the magnitude of the d/u ratio at x=1. The mechanism that results in the decrease of ${d_V\over u_V}$ ratio with an increase of x, is due to the inequality of Eq. (3) because of which the recoil wave function for the case of the d-quark decreases faster than that of the u-quark. Thus the value of the ${d_V\over u_V}$ ratio at $x\rightarrow 1$ in the present model is related to the difference in the peak positions for h(x, t) distribution for d- and u-quarks.

Using parameters of Table 1 one evaluates

$$\begin{aligned} 0.06 \le {d_V\over u_V}\mid _{x\rightarrow 1} \le 0.1 (0.14) \end{aligned}$$

(63)

where (0.14) is the estimate for CJ15 parameterization. This is the prediction of the current mean-field model if quark–quark correlations will be found to be negligible at large x.

In Fig. 14 we present the x dependence of the ${d_V\over u_V}$ ratio for the same parameters used in Fig. 13.

Fig. 14

The x dependence of ${d_V\over u_V}$ ratio in valence quark region, solid curve. The calculation is compared with the one following from CT14nnlo parametrization. The area takes into account the spread of the PDFs for d- and u-quark distributions [38]

If, however, the current estimates of d/u ratio from DIS processes off $A\!\!=\!\!3$ nuclei [62, 63] is confirmed, it will indicate the need of qq short range interactions to enhance the d/u ratio at $x=1$, to the values, $\sim 0.2$.

5 Summary and outlook

We have developed a model which describes the mean-field dynamics of valence quarks in the nucleon by separating the nucleon into a valence 3q-cluster and residual system. Within this model based on effective light-front diagrammatic approach we calculated the valence quark distributions in the nucleon in leading order approximation. The parameters entering the model are ones that describe the wave functions of nonperturbative 3q and residual systems.

5.1 Summary of results and predictions of the model

An important result of the model is the prediction of a new relation for the position of the peak in the $x q_V(x)$ distribution and the mass of the residual system according to Eq. (58) for moderately large $Q^2$. This relation naturally describes the fact that the valence d-quarks peak at lower $x_p$ than u-quarks following from the expectations, that in the proton, the residual system for the d-quark in general has larger mass than for the u-quarks (Eq. (2)).

The model predicts $(1-x_B)^4$ behavior for valence quark PDFs at $x_B\rightarrow 1$, which is stronger than one predicted in pQCD, $(1-x_B)^3$. This indicates the need of inclusion of 2q and 3q correlations to describe high x tail of valence PDFs.

To evaluate the parameters of the model we fitted the calculated valence quark distributions to the several leading order phenomenological PDFs with starting $Q^2=4$ GeV$^2$. The obtained parameters depend on the particular choice of the phenomenological PDF. However they demonstrate the following common features

• The mass of the residual system for probed valence u-quarks is systematically lower than that of the valence d-quarks and it is less than the pion mass

• The slope factor of the residual system in the proton is larger for the u-quarks than for the d-quark indicating that d-quark-residual system is more compact

• The magnitude of the parameters indicate the need for larger (by factor of two) high momentum component for the valence u-quarks than for the d-quark.

The above features allow us to make specific predictions that can be checked experimentally. One is that there should be a correlation between the multiplicity of pion productions in the target fragmentation region with the flavor of the knock-out quark in the forward direction. Specifically one expects more pions to be produced in the target fragmentation region of the proton if the forward scattering takes place on the valence d-quark. Another prediction of the model is that for asymmetric nuclei, attenuation of the residual system is correlated with the flavor of struck quark in the forward direction. This prediction is based on the observation that the recoil system is more compact (larger $B_R$) in the case of interaction with the valence d-quark in the proton.

5.2 Possible improvements and fitting to beyond the leading order PDFs

To be able to calculate the nucleon structure function analytically we treated the kinematics of residual system non-relativistically (Eq. (49)). However, the fitting results of $m_R$ indicate the need for a relativistic treatment of the residual system. This will be done in the process of fitting the model parameters to the next to leading order PDFs. For the latter case since all modern NLO PDFs adopt ${\overline{MS}}$ factorization scheme the relation between structure function $F_2$ and PDF, $f_i(x,t)$ is a convolution integral which covers region above the probed Bjorken x kinematics. Thus, fitting to the NLO PDFs will require the modeling of the contribution from qq correlations presented in Fig. 1b, c. These calculations are in progress and will be published in follow-up paper.

In the current model we described the residual system by a single mass, $m_R$ assuming the latter is independent of the momentum fraction of valence quark struck from the cluster. We also neglected the possible role of higher residual masses with the isospin different from that of the proton, which will correspond to different 3q cluster wave function. The model could be improved further with consideration of a more elaborate structure for the residual system. However this will require a comprehensive comparison of the full calculation with different experimental data.

Finally, the reliability of the obtained parameters for the 3q-cluster and residual system wave functions can be checked by calculating the mean-field effects at $x\sim 0.1-0.4$, in many different quantities in QCD, such as form-factors, generalized partonic distribution functions or transverse momentum distributions for different DIS processes.

With the addition of the two- and three-quark short-range interactions to the present model we hope to have a framework that will allow us to study the dynamics of the valence quarks in the entire range of $x>0.1$. The parameters obtained in this work for the wave functions of the mean field 3q-cluster (Eqs. (47)) and residual system (Eq. (48)) will allow us to use them in the calculation of the correlation diagrams of Fig. 1b, c. Note that in this work we did not use polarized PDF’s to fix the parameters of the model since one of the assumptions was that the mean-field dynamics adheres SU(6) symmetry. However once two- and three-quark short-range interactions are added they will break the SU(6) symmetry because of the vector exchange nature of interactions. Then by fitting calculated distributions to polarized PDFs we will be able to calibrate the contribution from short range interaction at large x.

Appendices

Appendix A: Light-front diagrammatic rules

The light-front diagrammatic rules are given for the processes ordered in light-cone time, $\tau = x^+ = t + z$ with the possibility of an instantaneous interaction.

• 1. Draw topologically distinct $\tau$-ordered diagrams identifying effective interaction vertices. In addition to the usual advanced and retarded propagation between two events one needs to include a possibility in which the two events connected by an internal fermion or photon lines that interact at the same LC $\tau -$time, commonly referred as an instantaneous term.

• 2. For each particle line assign four-momentum $p^\mu$ and spin s (or helicity, $\lambda$) with the on-mass-shell conditions $p^2=m^2$.

• 3. Incoming fermions/anti-fermions get $u(p,h) / \bar{v}(p,h)$ . Outgoing fermions/anti-fermions have $\bar{u}(p,h) / v(p,h)$.

• 4. Incoming photons or gluons with polarization $\lambda$ are assigned with polarization vector $\epsilon ^\mu _\lambda (k)$

• 5. For each internal particle line assign a factor $\frac{\theta (p^+)}{p^+}$

• 6. Intermediate states get a light front energy denominator of:

  $$\begin{aligned} {1\over {\mathcal {D}}} = \frac{1}{\sum _{init.}p^- - \sum _{interm.}p^- + i\epsilon } \end{aligned}$$

  where sums go over the initial and intermediate states.

• 7. Assign to the each vertex in the diagram the effective transition factor ${\varGamma }$. For bare particle elementary interactions these factors will correspond to the fundamental vertices with coupling constants.

• 8. For the transition of the composite particle, A with momentum $p_A$ and helicity $h_A$ to its n-constituents one introduces LF wave function in the form:

  $$\begin{aligned} \psi (\{ x_i,k_{i,\perp },h_i\}_i^n) = {\prod \limits _{i=1}^{n} \chi _{fi}(x_i,k_{i,\perp },h_i){\varGamma } \chi _A(p_A,h_A) \over p_{A}^+ {\mathcal {D}}},\nonumber \\ \end{aligned}$$

  (A.1)

  where $x_i$ and $k_{i,\perp }$ are LC momentum fractions and transverse momenta of constituent particles with helicity $h_i$. The spin wave functions of outgoing particles are described by $\chi _{fi}(x_i,k_{i,\perp },h_i)$ and ${\varGamma }$ is the effective vertex of the transition of particle A to n-constituents.

• 9. Each closed loop is integrated by the LC space factor

  $$\begin{aligned} \int \frac{dp^+d^2 \mathbf {p_\perp } }{16\pi ^3}. \end{aligned}$$

Appendix B: Integrations in the transverse momentum space

Now, if we introduce the relative transverse momentum ${\mathbf {k}}_{i, \perp } = {\tilde{\mathbf {k}}}_{i, \perp } - \frac{x_i}{x_V} {\mathbf {k}}_{V, \perp }$ and using the fact that

$$\begin{aligned} \delta ^{\left( 2\right) } \left( \sum _{i=1}^3 {\mathbf {k}}_{i, \perp } + {\mathbf {k}}_{R, \perp }\right) = \delta ^{\left( 2\right) } \left( \sum _{i=1}^3 {\tilde{\mathbf {k}}}_{i, \perp }\right) , \end{aligned}$$

(B.1)

one factorizes the $d^2{\tilde{k}}_{R,\perp }$ and $d^2 k_{i,\perp }$ integrations resulting in a following expression for the transverse integrals of Eq. (50) :

$$\begin{aligned}&\int ^{Q^2} [d^2{\mathbf {k}}_{\perp }] \exp \left[ -\frac{B_Vx_V}{4} \sum _{i=1}^3 \frac{k_{i,\perp }^2}{x_i} - B_R k_{R,\perp }^2 \right] \nonumber \\&\quad = \int ^{Q^2} \prod _{i=1}^3 \frac{d^2 {\tilde{\mathbf {k}}}_{i, \perp }}{16 \pi ^3} 16\pi ^3 \delta ^{(2)} \left( \sum _{i=1}^3 {\tilde{\mathbf {k}}}_{i, \perp } \right) \nonumber \\&\qquad \times \exp \left[ -\frac{B_V x_V}{4} \sum _{i=1}^3 \frac{{\tilde{k}}_{i, \perp }^2}{x_i}\right] \int ^{Q^2} \frac{d^2 {\mathbf {k}}_{R, \perp }}{16 \pi ^3} e^{-B_R k_{R, \perp }^2}. \end{aligned}$$

(B.2)

Furthermore, we take the $d^2{\tilde{\mathbf {k}}}_{3,\perp }$ integral using the $\delta ^2()$ function and evaluate the $d^2 {\mathbf {k}}_{R, \perp }$ integral analytically, arriving at:

$$\begin{aligned}&\frac{1}{(16\pi ^3)^2} \int ^{Q^2} d^2 {\tilde{\mathbf {k}}}_{1, \perp } d^2 {\tilde{\mathbf {k}}}_{2, \perp } \nonumber \\&\qquad \times \exp {\left[ -\frac{B_V x_V}{4} \left( \frac{{\tilde{k}}_{1, \perp }^2}{x_1} + \frac{{\tilde{k}}_{2, \perp }^2}{x_2} + \frac{({\tilde{\mathbf {k}}}_{1,\perp } + {\tilde{\mathbf {k}}}_{2,\perp })^2}{x_3}\right) \right] } \nonumber \\&\qquad \times \frac{\pi }{16\pi ^3} \int _0^{Q^2} dk_{R,\perp }^2 e^{-B_R k_{R, \perp }^2} \nonumber \\&\quad = \frac{1}{(16\pi ^3)^2} \int ^{Q^2} d^2 {\tilde{\mathbf {k}}}_{1, \perp } d^2 {\tilde{\mathbf {k}}}_{2, \perp } \nonumber \\&\qquad \times \exp { \left[ -\frac{B_V x_V}{4} \left( \frac{{\tilde{k}}_{1, \perp }^2}{x_1} + \frac{{\tilde{k}}_{2, \perp }^2}{x_2} + \frac{({\tilde{\mathbf {k}}}_{1,\perp } + {\tilde{\mathbf {k}}}_{2,\perp })^2}{x_3}\right) \right] } \nonumber \\&\qquad \times \frac{\pi (1-e^{-B_R Q^2})}{16\pi ^3 B_R}. \end{aligned}$$

(B.3)

Next we introduce CM and relative transverse momenta for “1” and “2” quarks:

$$\begin{aligned} {\mathbf {k}}_{12}^{cm}= & {} {\tilde{\mathbf {k}}}_{1, \perp } + {\tilde{\mathbf {k}}}_{2, \perp }, \quad {\tilde{\mathbf {k}}}_{12}^{rel} = \frac{x_2{\tilde{\mathbf {k}}}_{1, \perp } - x_1{\tilde{\mathbf {k}}}_{2, \perp }}{x_1 +x_2}, \end{aligned}$$

(B.4)

which allows to decouple the momentum dependent term in Eq. (B.3) as follows:

$$\begin{aligned}&\frac{k_{1, \perp }^2}{x_1} + \frac{k_{2, \perp }^2}{x_2} + \frac{({\mathbf {k}}_{1,\perp } + {\mathbf {k}}_{2,\perp })^2}{x_3} \nonumber \\&\quad = (k_{12}^{rel})^2 \left( \frac{1}{x_1} + \frac{1}{x_2} \right) + (k_{12}^{cm})^2 \left( \frac{x_V}{x_3(x_1 + x_2)} \right) . \end{aligned}$$

(B.5)

Using this and the relation $d^2{\tilde{\mathbf {k}}}_{1,\perp } d^2{\tilde{\mathbf {k}}}_{2,\perp } = d^2{\mathbf {k}}_{12}^{cm} d^2{\mathbf {k}}_{12}^{rel}$ for Eq. (B.2) one obtains:

$$\begin{aligned}&\int ^{Q^2} [d^2{\mathbf {k}}_{\perp }] \exp \left[ -\frac{B_Vx_V}{4} \sum _{i=1}^3 \frac{k_{i,\perp }^2}{x_i} - B_R k_{R,\perp }^2 \right] \nonumber \\&\quad = \int ^{Q^2} d^2{\mathbf {k}}_{12}^{cm} d^2{\mathbf {k}}_{12}^{rel} \nonumber \\&\qquad \times \exp \bigg [-\frac{B_V x_V}{4} \bigg ( (k_{12}^{rel})^2 \left( \frac{1}{x_1} + \frac{1}{x_2} \right) \nonumber \\&\qquad + (k_{12}^{cm})^2 \frac{x_V}{x_3(x_1 + x_2)} \bigg ) \bigg ] \nonumber \\&\qquad \times \frac{1}{(16\pi ^3)^2} \frac{\pi (1-e^{-B_R Q^2})}{(16\pi ^3) B_R}. \end{aligned}$$

(B.6)

Introducing $a_{rel} = \frac{B_V x_V }{4} \left( \frac{1}{x_1} + \frac{1}{x_2} \right)$ as well as $a_{cm} = \frac{B_V x_V}{4} \frac{x_V}{x_3(x_1 + x_2)}$, the remaining part of the integration in Eq. (B.6) can be evaluated analytically resulting in:

$$\begin{aligned}&\int ^{Q^2} [d^2{\mathbf {k}}_{\perp }] \exp \big [ -B_V \sum _{i=1}^3 x_V \frac{k_{i,\perp }^2}{x_i} - B_R k_{R,\perp }^2 \big ] \nonumber \\&\quad = \frac{1}{\left( 16\pi ^3\right) ^2} \frac{\pi \left( 1-e^{-a_{cm}Q_{cm}^{max 2}}\right) }{16 \pi ^3 a_{cm}} \frac{\pi \left( 1-e^{-a_{rel}Q_{rel}^{max 2}}\right) }{16 \pi ^3 a_{rel}} \nonumber \\&\qquad \times \frac{ \pi \left( 1-e^{-B_R Q^2}\right) }{16\pi ^3 B_R} \nonumber \\&\quad = \frac{x_1 x_2 x_3 }{\left( 16\pi ^3\right) ^4 x_V^3 B_V^2 B_R} \left( 1-e^{-a_{cm}Q_{cm}^{max 2}}\right) \left( 1-e^{-a_{rel}Q_{rel}^{max 2}}\right) \nonumber \\&\qquad \times \left( 1-e^{-B_R Q^2}\right) , \end{aligned}$$

(B.7)

where $Q_{cm}^{max}$ and $Q_{rel}^{max}$ denote the max values of $k_{12}^{cm}$ and $k_{12}^{rel}$, respectively. Substituting the above expression into Eq. (50) one obtains:

$$\begin{aligned}&f_q(x_B, Q^2) = {\mathcal {N}} \int [dx]\delta (x_B-x_1)\nonumber \\&\quad \times \exp \left[ -\frac{B_Vx_V}{4} \sum _{i=1}^3 \frac{ m_i^2}{x_i} - B_R M_N^2(x_R - \frac{M_R}{M_N})^2 \right] \nonumber \\&\quad \times \frac{x_1 x_2^2 x_3^2 x_R}{x_V^3} \left( 1-e^{-a_{cm}Q_{cm}^{max 2}}\right) \left( 1-e^{-a_{rel}Q_{rel}^{max 2}}\right) \nonumber \\&\quad \times \left( 1-e^{-B_R Q^2}\right) , \end{aligned}$$

(B.8)

where ${\mathcal {N}}= {16\pi ^3 A_V^2 A_R^2 m_N^3\over B_R B_V^2} e^{{9\over 4} B_V m_q^2}$. The two remaining Dirac delta functions can be used to evaluate the $d x_1$ and $d x_R$ integrals.