Correlations in Double Parton Distributions: Perturbative and Non-Perturbative effects

The correct description of Double Parton Scattering (DPS), which represents a background in several channels for the search of new Physics at the LHC, requires the knowledge of double parton distribution functions (dPDFs). These quantities represent also a novel tool for the study of the three-dimensional nucleon structure, complementary to the possibilities offered by electromagnetic probes. In this paper we analyze dPDFs using Poincar\'e covariant predictions obtained by using a Light-Front constituent quark model proposed in a recent paper, and QCD evolution. We study to what extent factorized expressions for dPDFs, which neglect, at least in part, two-parton correlations, can be used. We show that they fail in reproducing the calculated dPDFs, in particular in the valence region. Actually measurable processes at existing facilities occur at low longitudinal momenta of the interacting partons; to have contact with these processes we have analyzed correlations between pairs of partons of different kind, finding that, in some cases, they are strongly suppressed at low longitudinal momenta, while for other distributions they can be sizeable. For example, the effect of gluon-gluon correlations can be as large as 20 $\%$. We have shown that these behaviors can be understood in terms of a delicate interference of non-perturbative correlations, generated by the dynamics of the model, and perturbative ones, generated by the model independent evolution procedure. Our analysis shows that at LHC kinematics two-parton correlations can be relevant in DPS, and therefore we address the possibility to study them experimentally.


Introduction
Multi Parton Interactions (MPI) occur when more than one parton scattering takes place in one hadron-hadron collision. They have been defined long time ago [1], have been recently rediscovered and are presently attracting remarkable attention, thanks to the activity of Large Hadron Collider (LHC), where specific signatures are expected to be observed (see Refs. [2,3,4,5,6] for recent reports).
In particular, the cross section for hard double parton scattering (DPS), the simplest MPI process, depends on non-perturbative objetcs, the double parton distribution functions (dPDFs), describing the number density of two partons located at a given transverse separation in coordinate space and with given longitudinal momentum fractions. dPDFs encode, for example, the novel information on the probability that partons which are close to each other are faster, or slower, than those which are far from each other. They are therefore naturally related to parton correlations, as noticed several years ago [7], and represent a novel tool to access the three-dimensional (3D) nucleon structure, presently studied using electromagnetic probes [8,9]. The correlations in DPS are presently deeply investigated (see, e.g., [3,10,11]).
In addition to this non perturbative information, the knowledge of dPDFs, DPS and MPI in general could be very useful to constrain the background to the search of new Physics at the LHC, making their study very timely. No data are presently available for dPDFs and their calculation using non perturbative methods is cumbersome. A few model calculations, able in principle to grasp the most relevant features of dPDFs, have been therefore performed [12,13,14,15,16]. In particular, in Ref. [14], a Light-Front (LF) Poincaré covariant approach, reproducing the essential sum rules of dPDFs without ad hoc assumptions and containing natural two-parton correlations, has been described. We note in passing that, although it has not yet been possible to extract dPDFs from data, the so called "effective cross section", σ ef f , the ratio of the product of two single parton scattering cross sections to the DPS cross section with the same final states, has been extracted, in a model dependent way, in several experiments [17,18,19,20,21,22]. Despite of large error bars, the present experimental scenario is consistent with the idea that σ ef f is constant w.r.t. the center-of-mass energy of the collision. In Ref. [23] we have presented a predictive study of σ ef f , making use of the LF quark model approach to dPDFs developed in Ref. [14]. It was found that the order of magnitude of the measured σ ef f is correctly reproduced by the model and, more interestingly, in the valence region, a clear dependence is predicted on the longitudinal momentum fractions of the proton carried by the two partons. If measured, this feature could represent a first access to the observation of 2-partons correlations in the proton.
Beyond these intriguing results, already found in the valence region, one should check if similar possibilities survive at LHC kinematics, dominated by low-x partons, at very high energy scales. In this paper, using our model predictions, we plan therefore: i) to test the validity of factorization assumptions, which basically neglect at least part of the correlations between the partons, often used in dPDFs studies, at the scale of the model and after evolution to experimental energy scales; ii) to test if correlations in longitudinal and transverse momenta survive the evolution procedure; iii) to develop an extension of our approach to include, at the low energy scale of the model, sea quarks and gluon degrees of freedom; iv) to study the importance of 2-body correlations between different kinds of partons (valence quarks, sea quarks and gluons) at values of longitudinal momenta and energy scales close to the experimental ones, to establish the possibility to observe them at the LHC.
The paper is structured as follows. The first section is dedicated to present a short summary of the formalism and the results obtained in Ref. [14]. The second section is dedicated to compare our results, where correlations are naturally produced by the dynamics of the model, with a few factorized forms of dPDFs. In the third section, we study how QCD evolution to high momentum scales affects the results of the model. In the following section we describe a strategy to introduce sea quarks and gluons at the low momentum scale of the model. In section five we quantify, within our scheme, how large are the correlation effects between different kind of partons at very low values of longitudinal momentum fractions and very high energy scales. This is very important to address measurable signatures of two-partons correlations. We end by drawing some conclusions of our study.

Calculating Double Parton Distribution Functions
Recently dPDFs have been explicitly calculated by us within a Light-Front (LF) approach [14]. The method is fully covariant and is based on a fixed-number Light-Front SU (6)-symmetric Hamiltonian making use of an Hypercentral potential introduced in Ref. [24] as a generalization of a non-relativistic constituent quark model proposed in Ref. [25]. The approach is particularly suitable for the description of Deep Inelastic Scattering (DIS) processes which find their natural environment in a LF -description.
The numerous applications to a large varieties of DIS observables like polarized [24] and unpolarized [26,27,28] structure functions, spin and angular momentum distributions [29,30], helicity-independent and dependent GPDs [31], demonstrate the reliability and flexibility of the approach.

The Light-Front formulation
Let us briefly summarize the main steps for the LF-evaluation of the dPDFs. In terms of the Light-Cone (LC) quantized fields q i for a quark of flavor i, helicity λ in an unpolarized proton, the dPDFs in momentum space, often called " 2 GP Ds" in the literature [32,33], read (see, e.g., [12,13]) andÔ In the above equations, both the light-like four vectorn = (1, 0, 0, −1) and the rest frame state of the nucleon with helicity λ, P = 0, λ , have been introduced. The "±" components of a four-vector b are defined according to b ± = b 0 ± b z and x i = k + i /P + is the fraction of the system momentum carried by the parton "i", while the notationb = (b + , b ⊥ ) is used for light-cone vectors. The LC free quark fields are defined as where the operator a ĩ k,r destroys a quark of flavor i, helicity r and LC momentumk. The spinors are indicated by u LF (k, r) (we adhere to the definitions and notations of Ref. [34]). The proton state P = 0, λ can be expanded in its Fock components retaining only the first (valence) contribution (the short-hand notation ({α i }) is adopted, here and in the following, for (α 1 , in terms of the LF one-quark states of isospin τ i , |k i , λ f i , τ i . At variance the same proton state can be described in terms of canonical, The two descriptions are related by Melosh rotations [35]. Following our previous developments (e.g. Refs. [24,31]) the considerations made for free canonical states can be generalized to interacting quarks in a proton, by means of a suitable representation of the Poincaré operators, namely the Bakamjian-Thomas construction [36]. The extension to interacting systems requires, in fact, a dynamical representation of the Poincaré group. One way to achieve this result is to add an interaction V to the free mass operator M 0 to obtain the mass operator M = M 0 +V . Since the LF boosts we use are interaction independent, all the other definitions remain unaffected. All required commutation relations are satisfied if the mass operator commutes with the total spin and with the kinematic generators. In practice, the conditions are realized if: i) V is independent on the total momentumP; ii) V is invariant under ordinary rotations.
Summarizing: in the LF formulation of the quark dynamics, the intrinsic momenta of the quarks (k i ) can be obtained from the corresponding momenta (p i ) in a generic frame through a LF boost (K i = u(P ) · p i , P ≡ 3 i 1 p i ) such that the Wigner rotations reduce to the identity. The spin and spatial degrees of freedom are described by the wave function where µ i refers to the eigenvalue of the LF spin, so that the spin part of the wave function is transformed by the tensor product of three independent Melosh rotations, namely R † = 3 ⊗i=1 R † (k i , m i ). The internal wave function is an eigenstate of the baryon mass operator i + m 2 i and where the interaction term V must be independent on the total momentum P and invariant under rotations. The nucleon state is then characterized by isospin (and its third component), parity, Light-Front (non-interacting) angular momentum operators with well defined projection along the quantization axis.
The relativistic mass equation chosen is built according to such a dynamical construction [24]. Thanks to the correct kinematical conditions on the longitudinal momentum fraction carried by the quark as described by the LF-approach, dPDFs vanish in the forbidden kinematical region, x 1 +x 2 > 1. (see Ref. [14] for further details).

Light-Front results at the low scale of the model
Reducing Eq. (1) to the first (valence) Fock components and specializing the result to the u quarks as an example, one has (λ 1 , with The spin projector values are determined by the Melosh rotationsD i : to be calculated using the canonical spin-isospin states corresponding to the SU(6) symmetric matrix elements.
In particular the combinations will describe two unpolarized u-valence quarks, and two (longitudinally) polarized u-valence quarks. These two distributions only contribute to the total cross section of events involving unpolarized proton targets. In Fig. 1 the numerical results of the Eqs. (8)- (11) for two unpolarized u-valence quarks. The dPDFs vanish in the region x 1 + x 2 > 1 and the correlations in x 1 , x 2 are dictated by the LF-quark dynamics, which governs also the dependence in k ⊥ , clearly seen in the right panel of the same figure.
Since the Fock expansion of the proton state has been restricted to the three valence quarks (cfr. Eq. (5)), it is natural that the full momentum is carried by those quarks and the resulting appropriate energy scale remains quite low (the so-called hadronic scale, µ 2 0 ≈ 0.1 GeV 2 [37]) as indicated by analogous Leading − Order calculations (see, e.g., Refs. [24,27,28]). That scale is clearly indicated in the resulting expressions (11) and (12) and in both panels of Fig. 1.

2.3
Factorization and approximations at the low scale of the model In the present Section we will compare our approach to a number of strategies used in the literature to calculate dPDFs, strategies that we call, in general, "factorization schemes". Differences and analogies will help in understanding the role of correlations and their dependence on the evolution scale.

Phenomenological factorizations
As a first illustrative example we can restrict the discussion to the approach proposed by Diehl, Kasements and Keane in Ref. [38], which has motivated in part the discussion presented in this section.
In fact the interest of those authors is on the influence of the evolution scale on correlation effects, precisely one of the goals of the present work. The model they propose refers to the description of the dPDFs at the starting scale, where they assume independent partons. In that case, in fact, the dPDFs in coordinate space can be simply written as a convolution of f a,b (x, b) functions, which are impact parameter dependent generalized parton distributions (see, e.g, Ref. [3]): with a, b denoting parton species. This idea has been firstly presented in Refs. [32,33]. The authors of Ref. [38] assume a Gaussian b dependence with an x-dependent width, namely  where f a (x) denotes the usual parton densities (taken from the LO set of the MSTW 2008 analysis [39]), while Eq. (14) is assumed to be valid at the starting scale Q 2 0 = 2 GeV 2 . Diehl et al stress that the approach is tailored for the region x 1 , x 2 < 0.1 and its parameters are specified for gluons and for the sum, q + = q +q, and difference, q − = q −q, of quark and antiquark distributions. The expressions for h a (x) are found in Ref. [38] and not reported here; the parameters which are necessary to define h a (x) are fixed so that the resulting parton densities are in tentative agreement with phenomenology.
The final expression for the unpolarized dPDFs Eq. (13) reads: and, as a consequence, one has, for the Fourier transform, The term is assumed, at the same scale, to introduce correlations between x 1 and x 2 , in fact Eq. (17) does not factorize into separate contributions from each of the two partons, a and b (the values of the parameters in Eq. (17) can be found in Ref. [38]). The combinations u − and u + are taken as representatives of the quark sector.

Factorization by means of Generalized Parton Distributions
In Ref. [3], a systematic study of relations between single parton and double parton distributions has been performed. To reduce F ab to single-particle distributions the authors find it more convenient to work in the transverse-momentum k ⊥ space, rather than transverse distance y representation and the result reads: where M p is the proton mass and H q (x, ξ, t) and E q (x, ξ, t) are Generalized Parton Distributions (GPDs) (see, e.g., [40] and references therein). H q generalize the unpolarized quark densities q(x) while E q is related to unpolarized quarks in a transversely polarized proton. The first term in Eq. (18) depends on H q only and it corresponds to the simplest approximation of the two-parton distribution as a product of single-parton distributions (cfr. Eq. (16)). In Ref. [3] one can read: "Although the relation between multiparton distributions and GPDs is an approximation whose accuracy is not easy to estimate (and although our current knowledge of GPDs is far less advanced than that of ordinary parton densities) this relation provides opportunities to obtain information about multiple interactions that is hard to get by other means".
The question on the accuracy is particularly relevant in view of possible experimental studies of multi-parton effects, and the LF-approach we are presenting can shed some light on the approximation (18), including the role played by the E correction term. Since GPDs have been studied, precisely within the same LF-approach, by Pasquini, Boffi and Traini [31] some years ago, one can check directly the accuracy of Eq. (18) 1 . Because of the natural normalization of the expression Eq. (11): and the normalization of the H u V GPDs the comparison holds for       A detailed comparison is shown in Fig. 3 and Fig. 4 for two specific values of x 2 , namely x 2 = 0.1 and x 2 = 0.3. In these two cases the comparison is not restricted to k ⊥ = 0 only (see Fig. 3), but it extends to the kinematical region up to k 2 ⊥ = 0.5 GeV 2 (see Fig. 4). Once again no systematic agreement is found. The only weak improvement, for k ⊥ > 0, is due to the E u V dependent correction term.
3 Scale impact on correlations 3.1 Analysis of the approach of Ref. [38] Let us first analyze the scale dependence of the correlations introduced at Q 2 0 within the assumptions Eqs. (15) and (16), as proposed by Diehl et al. in Ref. [38]. To this end, we study the QCD-evolution of the dPDFs. In particular one could ask oneself to which extent the Gaussian y-dependence (or k ⊥ -dependence) of the starting scale is preserved under evolution. Quantities particularly suitable to this end are the ratios which, at Q 2 0 and x 1 = x 2 = constant, are just straight lines as functions of y 2 or k 2 ⊥ . Perturbative evolution of the dPDFs is summarized and discussed in Appendix A; however let us anticipate the results in this example, proposed in Ref. [38]. As it is done also in Ref. [38], only the homogeneus part of dPDFs evolution is implemented, for the moment being, in our scheme. According to some studies, the inhomogeneus part could play some role in this phenomenology [33,41]; its analysis is beyond the scope of the present paper.
In Fig. 5 we show the ratios Eqs. (22) for different quark and gluon combinations: F u + u + , F u − u − and F gg at x 1 = x 2 = 0.1 at different scales, namely the starting scale Q 2 0 = 2 GeV 2 and Q 2 = 10 4 GeV 2 . One can check that, for quarks, the shape remains approximately Gaussian (a straight line)  Figure 5: Effects of evolution on correlations according to the scheme of Ref. [38]. Upper panel: at fixed values of Q 2 and following the assumptions of Ref. [38]. Middle panel: As in the upper panel for ln[F u + u + ( y 2 )/F u + u + (0)]. Lower panel: As previous panels, for ln[F gg ( y 2 )/F gg (0)]. up to scales as high as Q 2 = 10 4 GeV 2 even if the slope changes rather strongly. For gluons, also the Gaussian property is not preserved.  Figure 6: Effects of evolution on correlations for the LF-Hypercentral approach.
In particular the upper panel of Fig. 5 shows the valence components of . For those distributions, only the non-singlet evolution is relevant. In the case of

Scale dependence within the LF-approach
In the LF-approach the Fock decomposition of the proton state at the lowest scale µ 2 0 includes valence quarks only and one remains with the following reductions with u V u V (x 1 , x 2 , k ⊥ , µ 2 0 ) given by Eqs. (8) and (11). The distributions are now function of k ⊥ since in the LF-approach they are defined in momentum space. The relation with y is a simple Fourier transform (cfr. Eqs. (15) and (16)); however their functional forms, entirely determined by the dynamical structure of the LF-wavefunctions, are far from being Gaussian.
The distributions at the starting point are strongly simplified as indicated by Eqs. (24), but they evolve in a complicated way as combination of non-singlet Fig. 6, the results are shown at fixed x 2 = x 1 = 0.1, Q 2 = 10 4 GeV 2 , as function of k 2 ⊥ . The form is clearly non-Gaussian, since non-Gaussian is its functional form at µ 2 0 , fact which is related to the dynamics of the LF, not on the value of µ 2 0 . The complete results are shown in Fig. 6, following the same notations and criteria of Fig. 5.
Comparing the results of the two set of Figures (5 and 6), it is evident that the evolution effects are similar in the two different cases, but it is also evident that the Gaussian ansatz is rather arbitrary and not supported by LF dynamics.

4
Adding sea quarks and gluons at a low energy scale In the previous sections, our Light-Front approach has been focused on the study of valence degrees of freedom at low momentum scale. Other partons, and their correlation effects, emerge from radiated gluons in the perturbative QCD-evolution of the dPDFs. In the present Section we enlarge the perspective studying how sea quarks and gluons can be included at a low-momentum scale and within the same LF framework. An example (e.g. Refs. [27,28]) is given by inclusive DIS, where the (nonperturbative) meson degrees of freedom can be introduced by means of a description of the meson cloud and the scattering of the virtual photon off the constituents of the mesons (Sullivan process). Analogous approach can be applied to the explicit evaluation of meson cloud effects on GPDs (e.g. Refs. [31,42]).
Hereafter we will propose a simplified approach in which the effects of the valence degrees of freedom (producing the largest part of the dPDFs at low-momentum scale) are calculated using Eqs. (8) and (11), while the non-perturbative sea and gluons components are evaluated by means of a factorized approximation of the kind discussed in Sect. 2.3. In order to minimize the hypothesis on factorization let us start discussing the limiting case k ⊥ = 0. Let us first illustrate, as an example, the uu dPDFs: The pure valence (and dominant) term, the expression (25) in the above equation, can be evaluated in a direct way within the LF-approach described in the previous sections. In order to calculate the residual terms, Eqs. (26) and (27), one can assume factorized forms (see e.g. Ref. [3]). The complete (approximate) expression for F uu becomes: Few comments are in order: i) the contribution Eq. (30) is the term due to valence quarks, it is not approximated by a factorized procedure and it is based on the calculated expressions Eqs. (8) and (11); ii) the residual contributions imply the knowledge of the singlet componentū(x, Q 2 0 ) and fulfill the correct kinematical conditions for x 1 + x 2 > 1, owing to the constraints introduced by the phenomenological function (1 − x 1 − x 2 ) n θ(1 − x 1 − x 2 ). The exponent n has to be fixed phenomenologically, as seen in Sect. 2.3.1 in the case of the model of Ref. [38] and will be discussed in the next Sections for the LF-approach; iii)ū(x, Q 2 0 ) = u sea (x, Q 2 0 ) has, at the low momentum scale Q 2 0 , a non-perturbative origin, basically due to the meson cloud surrounding the nucleon; iv) the scale Q 2 0 is not to be identified with µ 2 0 , i.e. the scale of the bare nucleon, where only the three valence quarks contribute.
In the following sections we will discuss a straight-forward (phenomenological) way of introducing meson and gluon degrees of freedom at the low-momentum non-perturbative scale. QCD evolution will be used to reach the high energy scale of the LHC experiments.

4.1
Factorization procedures within the LF-approach at k ⊥ = 0 The advantage of the approach we are discussing is based on a complete calculation of correlation effects within the LF-dynamics, in the restricted space of valence degrees of freedom. At the same time it allows to discuss the role of the factorization procedure and its validity, comparing our approach with phenomenological factorized models. This comparison aims to identify the coherence and self-consistency of the factorization schemes. In the following we will

Fixing the factorization form
The optimization of factorization procedures for dPDFs is not a simple issue. The most relevant constraints are related to momentum and quark number sum rules [2]. Our LF-approach, on the contrary, fulfills such sum rules by construction and therefore one does not need to implement phenomenological assumptions required to build factorized dPDFs.
As an example the resulting valence dPDF u V u V (x 1 , x 2 , k ⊥ = 0, µ 2 0 ), as well as the single PDFs (sPDFs) calculated within the same LF dynamical approach, fulfill the momentum and quark number sum rules. One can take advantage from such fundamental properties to fix the order of magnitude of the phenomenological exponent n in Eq. (31), trying to combine the knowledge of sPDFs and dPDFs in the following (factorized) relation: The restricted validity of the factorization approach has been already discussed in Sect. 2.3, therefore one cannot expect Eq. (32) to be satisfied with a high degree of accuracy. We expect, however, indications for the value of the exponent n to be used for building the additional sea and gluon contributions to the dPDFs. The value n = 2 has been discussed in the past as a good choice (see, e.g., Ref. [43] and references therein). More recent arguments (see, e.g., Refs. [38], and [2]) are in favor of more sophisticated parametrizations. Given the restricted use we are going to make of the factorization assumption, we prefer to remain within the straight-forward formulation Eq. (32). Our numerical analysis confirms a limited validity of the factorization and, at the same time, suggests n ≈ 0.2 (more precisely, values within the range 0.1 < n < 0.5; n = 0.2 is our optimal choice).

Sea and gluon contribution according to Ref. [39]
The advantages of the approximation Eq. x 1 x 1 can be approximated within a clear and self constrained factorized approach able to select the form of the factorization as discussed in the previous subsection.
In the following we discuss the introduction of sea and gluon degrees of freedom by means of one of the most used phenomenological parametrization of sPDF, the LO MSTW2008 parametrization (see Table 4 of Ref. [39]). The parametrization is valid at Q 2 0 = 1.0 GeV 2 . The fact that we are proposing LO parametrization is specifically due to the evolution properties of the dPDFs, known at LO only.
At the scale Q 2 0 , the sPDF MSTW2008 parametrization is characterized by the presence of partons like u V , d V ,ū,d, s,s and gluons. The total momentum is shared among such degrees of freedom and one has: Since at the scale µ 2 0 the system is determined by the valence degrees of freedom only, one has: dxx u V (x, µ 2 0 ) + d V (x, µ 2 0 ) = 1; Sea(x, µ 2 0 ) = 0; g(x, µ 2 0 ) = 0. As a consequence, to use Eq. (32) at the scale Q 2 0 , the dPDF u V u V (x 1 , x 2 , k ⊥ = 0, µ 2 0 ) at the scale µ 2 0 has to be evolved to Q 2 0 . The evolution is performed by means of the Non-Singlet reduction of the QCD evolution as described in Ref. [14] and summarized in the appendix A, where the complete Mellin procedure we are proposing for both Singlet and Non-Singlet sectors is illustrated in some detail. The result is shown in Fig. 7 for four selected values of x 2 . x 1 x 1 Fig. 7 is compared with x 1 x 2 F uu (x 1 , x 2 , k ⊥ = 0, Q 2 0 ) which contains the additional sea contributions fromū (u = u V +ū, see text) (continuous lines). x 1 In Fig. 8 and Fig. 9 we show the complete set of dPDFs involving the u-quark at the starting scale Q 2 0 . In particular, the combination Fig. 8 giving explicit evidence to the contribution due toū quarks. The largest effects of the sea component are clearly evident for smallest values of the momentum fraction. In Fig. 9 we show the dPDFs (F ug + F gu )/2 containing valence, sea and gluon contributions. The order of magnitude of those components is comparable with the valence part F u V u V at the scale Q 2 0 (cfr. Fig. 7).

4.2
Factorization procedures within the LF-approach at k ⊥ > 0 The dPDFs in a pure valence scenario have been discussed in previous sections and the dependence on k ⊥ has been explicitly investigated (cfr. for example, Figs. 1, 4, and 6). As a result they do not admit simple factorized forms. However, as a first attempt to go beyond the valence scenario at k ⊥ = 0, we could add the other degrees of freedom using factorized expressions. To this aim, the knowledge of the exact LF valence component at k ⊥ = 0 helps to define the additional, factorized contributions. In this section we find a reasonable factorized approximation to the exact valence LF dPDFs. In this way we fix the parameters which will be used for the non-valence degrees of freedom In practice, we want to generalize to k ⊥ > 0 Eqs. (31) and (32), valid at k ⊥ = 0. For instance, Eq. (32) becomes and simple choices are (cfr. Eq. (16)), Within scenario A of Eq. (37), no correlations between x 1 , x 2 and k ⊥ have been introduced: φ A depends on k ⊥ and it does not depend on x 1 , x 2 ; in scenario B of Eq. (38) the exponent depends on x 1 , x 2 , similarly to Eq. (16). The knowledge of u V u V (x 1 , x 2 , k ⊥ , µ 2 0 ) LF from Eq. (11) and of the sPDF u V (x 1 , µ 2 0 ) LF with the additional information n = 0.2 from the analysis of Sect. 4.1.1, can be used to optimize the fit Eq. (36). The results of the optimization procedure are shown in Fig. 10 for selected and extreme examples. Our recommended values are b A = b B = 0.6 GeV −1 and the quality of the fit is, once again, quite poor with a slight preference for the full correlated approximation of scenario B (Eq.(38)). The approximation is crude for the valence-valence correlations, but it is sound and the next subsection will be devoted to the implementation of additional degrees of freedom on the basis offered by the factorization Eq. (36) and scenario B.

4.2.1
Sea and gluon contribution at k ⊥ > 0 The valence-valence dPDFs do not allow for a simple factorization. The approximation proposed in Eq. (36) is therefore quite poor in that case, as Fig. 10 explicitly shows. However, we do not need to approximate valence-valence correlations; we can resort to the exact calculation also in the case of k ⊥ > 0 (cfr. Fig. 1 and Eqs. (8), (11)) and, using the best factorization scheme (n = 0.2 and  37), (38)) to introduce the additional degrees of freedom at the scale Q 2 0 and k ⊥ > 0. We have to follow once again the steps i)-iv) of Sect. 4 and the procedure described in Sect. 4.1.2. Numerically they are more challenging and, in a sense, incomplete since the evolution in k ⊥ is still an open problem [3]. We simply evolve at fixed k ⊥ applying the scale evolution of Appendix A. As an example we show in Fig. 11 the distribution F uu (x 1 , x 2 , k ⊥ , Q 2 0 ) at the scale Q 2 0 which generalizes Eq. (31) and Fig. 8, and whose resulting expression reads (from now on, we will discuss the scenario B of Eq. (38) only): is obtained evolving at the scale of the MSTW parametrization, Q 2 0 = 1.0 GeV 2 , the LF result (8), (11) at fixed k ⊥ . Besides, u V (x, Q 2 0 ) is the PDF obtained at the same scale within the LF approach, andū(x, Q 2 0 ) is taken from the LO MSTW parametrization [39]. By means of the perturbative evolution developed in Appendix A one can now evolve the distribution calculated at low-momentum scale to a typical experimental scale. We evolve to Q 2 = 250 GeV 2 , a scale relevant to study properties of dPDFs, as shown by experiments [17,18,19,20,21,22] and by a quite recent theoretical study within the LF approach [23]. In two series of figures, (Figs. 12 and  13), we compare the results obtained evolving directly from the lowest scale µ 2 0 where only valencevalence dPDFs are present (cfr. Fig. 1), with those obtained with the evolution from the scale of the MSTW parametrization Q 2 0 = 1 GeV 2 , where also gluon and sea dPDFs contribute (cfr. Fig. 11). The presence of the additional Singlet components is quite relevant, in particular for those components containing sea and gluon degrees of freedom, as it appears clearly from the comparison of the set of   collisions at very high energy. On the other hand the evolution obtained from the lowest momentum scale has the merit of being directly connected with quark dynamics and correlations are generated in a transparent way. A detailed study of the interrelations between non-perturbative correlations, generated by the dynamics of the model, and perturbative ones, generated by QCD evolution, is performed, at low-x, in the next sections.

5
Perturbative and non-perturbative two-parton correlations at low-x In this section we present results obtained within our LF scheme, aimed at establishing what kind of error one can do if two-parton correlations are neglected in treating dPDFs, for example in analyzing collider data. In previous papers of ours [14,23] we have already emphasized that this error can be rather sizeable when x 1 , x 2 lie in the valence region. We want now to analyze the low x scenario, reaching x values as low as 10 −2 , using a full, non-singlet and singlet, LO QCD evolution to the very high Q 2 scales typical of pp scattering at the LHC. As in Ref. [38], for the moment being, only the homogeneus part of the evolution of dPDFs is performed. As already said, the Q 2 evolution of the k ⊥ dependence has not been investigated yet and is still a missing item in this phenomenology.

5.1
Characterizing the two-parton correlations at low-x To study the relevance of two-parton correlations at low-x, we found very helpful to show ratios of dPDFs to products of PDFs; in the case of gluon distributions, for example at x 2 = 0.01, this ratio reads where Q 2 is a phenomenologically relevant scale, chosen in the following to be Q 2 = 250 and 10 4 GeV 2 . These scales are reached by performing QCD evolution of the results obtained within our LF scheme for both sPDFs and sdPDFs, starting from the hadronic scale µ 2 0 , where only valence degrees of freedom are present. It is clear that this ratio would be just 1 if it were possible to approximate the dPDF with the product of two sPDFs. The difference from 1 of the ratio is a measure of the error which is done by using that approximation, which amounts to disregard any kind of two-parton correlations.
In general the ratio can be written including other kind of partons; in the following, we will analyze the selected combinations The symmetrization is mandatory from the point of view of the experimental measurements, which cannot distinguish the two combinations. Obviously u V is a N on − Singlet-index, as well as g is a Singlet-index, while the sea indexes have no fixed flavor-symmetries; the different distributions evolve following the corresponding equations, as discussed in Appendix A.
Results of the ratio Eq. (43) for the flavor combinations gg, u V u V , u Vū ,ūū, u V g are shown in Figs. 14 -16.
All the ratios have two common qualitative features: i) results at Q 2 = 250 GeV 2 do not really differ from those at Q 2 = 10 4 GeV 2 ; the role of correlations does not depend therefore on the different high momentum scale which is chosen; ii) in all flavor combinations, when at least one of the momentum fractions of the two partons is in the valence region, correlations are strong and the error which is done in approximating a dPDF with a product of sPDFs is huge.  When both the momentum fractions of the partons are small, the situation is more involved. In the valence-valence sector, one finds negligible correlations and the ratio is basically 1 (cf. Fig. 14,  right panel). This fact, in the Non-Singlet (NS) sector, had been already found and discussed in Ref. [14]. In all other cases, where singlet evolution is playing a role, even at values of x 1 , x 2 as low as 10 −2 , correlations are found to produce sizable deviations of the ratios from 1. The maximum effect is found in the gluon-gluon case (cf. Fig. 14, left panel), when it reaches 20 %. One should realize that, if two-parton correlations were present at the LHC scale, one could access through DPS studies novel information on the proton structure. Our evolved model results show that if one were able to measure dPDFs at a 20 % accuracy, a specific dynamical information would be reachable. The different behavior of the valence-valence sector from the others, as well as the fact that the gluongluon sector experiences the biggets effect, are interesting features of our results and deserve to be understood through a further investigation. This is carried on in the next section.

Perturbative versus Non-Perturbative Two-Parton Correlations
In this subsection we will find that the results described in the previous one can be understood by disentangling perturbative and non-perturbative effects.
To this aim, let us consider again the ratio Eq. (43) In the previous subsection 5.1, results are obtained evolving the numerator from µ 2 0 to Q 2 , considering at the lowest scale the dPDFs predicted by our LF-model. The denominator is obtained evolving to Q 2 the analogous sPDFs of the same LF-model.
A first consideration is in order: if the denominator, given by the product of single PDFs, had been evolved by means of dPDF-evolution criteria, we would have obtained a simplified approximation of the dPDFs at Q 2 , including perturbative correlations only.
Let us define the following quantity In fact, F ab (x 1 , x 2 = 0.01, k ⊥ = 0, Q 2 ) P erturbative contains those correlations which come from dPDF perturbative evolution only. At this point, we could consider three different ratios: i) the ratio ab , Eq. (45); ii) the ratio P erturbative ab : ratio P erturbative which would be strictly 1 if only perturbative correlations were included in the numerator.
The three ratios are very useful to disentangle the effects of perturbative versus non-perturbative double-parton correlations; of course the ratios (43) or (45) are the most complete, including both kind of correlations in a consistent way.
In Figs. 17, 18 and 19, the results for the three ratios are compared at the scale Q 2 = 250 GeV 2 , at x 2 = 0.01, as functions of x 1 .
The ratio gg , shown in Fig. 17 (left panel), is particularly emblematic. The full ratio gg of Eq. (45) (dashed line), clearly influenced by both perturbative (dot-dashed line) and non-perturbative (continuous line) effects, is compared with those where perturbative and non-perturbative correlations are disentangled, contributing to the behavior of gluon − gluon dPDFs at low values of x 1 and x 2 . The same comments hold for dPDFs corresponding to the other partons. An interesting feature of these results, clearly read in Figs. 17, 18 and 19, is that in few cases the perturbative and non-perturbative components tend to cancel (in the case of V alence − V alence illustrated in Fig. 17 (right panel), or u −ū, as it borns out from Fig. 18 (right panel)). In the case of the gluon-gluon sector, the effect tends instead to sum coherently: this explains the persistence of correlations in this sector, even at high Q 2 and low x, observed in the previous subsection.
In closing this Section, we conclude that correlations in dPDFs, for some flavor combinations, are present also at low x 1 and x 2 , even at the large energy scale of LHC experiments. This arises because perturbative and non-perturbative effects sum coherently. These conclusions are not artifacts of the specific LF model used. They hold qualitatively also in ratios obtained starting the evolution from Q 2 0 = 1 GeV 2 >> µ 2 o , using as non-perturbative input the semi-factorized model of Section IV. In order to illustrate this important point, two more plots are included (Fig. (20)). In the first one, the valence-valence ratio is shown, in the other the gluon-gluon one. These examples are illustrative, indeed, of two specific aspects: i) the valence-valence ratios should not depend on the starting point because they converge at the same values at the common hadronic scale µ 2 0 . The small differences which appear in the figures are therefore a clear estimate of the errors introduced by our numerical evolution and one can appreciate the precision of our approach; ii) the second figure, showing the gluon-gluon ratio, is included because the glue is the dominant component at low-x and it contributes in a negligible way to the valence region. The correlations induced at low-x still contain a specific sign of the correlations introduced in the valence sector and this is due to the presence of the valence component in the quark-singlet sector in the evolution procedure. The strength of the correlation seems to become smaller but they are still sizable.

Conclusions
Double Parton Scattering (DPS) represents a background in several channels for the search of new Physics at the LHC. Its correct description depends on our ability of modelling double parton distribution functions (dPDFs). The knowledge of these quantities would represent also a novel tool for the study of the three-dimensional nucleon structure, complementary to possibilities offered by electromagnetic interactions, in the framework of Generalized Parton Distribution functions. In this paper we have analyzed dPDFs, using Poincaré covariant predictions obtained, at a low energy scale, within a Light-Front model proposed by us in a recent paper, evolved using QCD evolution to experimentally relevant scales. We checked to what extent factorized expressions of dPDFs, in terms of products or convolutions of one-body densities, can be used, neglecting, at least in part, two-parton correlations. Our tests were performed using our model predictions starting from a scale where only quark degrees of freedom are relevant, or from higher scales, modeling sea quark and gluon contributions. Our model study demonstrates that factorization procedures strongly fail in reproducing the calculated dPDFs in the valence region, where measurements of DPS could really allow to access two-parton correlations. Besides, a gaussian behavior for the transverse distance in coordinate or momentum space seems rather arbitrary. Anyway, to have contact with measurable processes at existing facilities, everything has to be pushed to very low values of the longitudinal momenta of the interacting partons. This study has been carried on systematically and represents the most interesting part of our investigations. Correlations between pairs of partons of different kind have been considered, finding that, in some cases, their effect tends to be washed out at low-x, as it happens for the valence, flavor non-singlet distributions, while they can affect other distributions in a sizable way, as in the gluon sector, when they can be as large as 20 %. We have shown that this different behavior can be understood in terms of a delicate interference of non-perturbative correlations, generated by the dynamics of the model, and perturbative ones, generated by the model independent evolution procedure. Our analysis shows that at LHC two-parton correlations can be relevant in DPS, opening a possibility to observe them for the first time. Our model dPDFs have now to be used to predict cross sections in specific channels where DPS is known to give an important contribution, such as, for example, the production of two W bosons with the same sign. Our research is now addressing this final goal.

Acknowledgements
This work was supported in part by the Mineco under contract FPA2013-47443-C2-1-P, by GVA-PROMETEOII/2014/066 and SEV-2014-0398. S.S. thanks the Department of Theoretical Physics of the University of Valencia for warm hospitality and support. M.T. and V.V. thank the INFN, sezione di Perugia and the Department of Physics and Geology of the University of Perugia for warm hospitality and support. Several discussions with F.A. Ceccopieri are gratefully acknowledged.

A
Appendix: Perturbative evolution of dPDFs in Mellin space Following Diehl and Kasemets [44] one has to admit that "a consistent formulation of factorization for double parton scattering does not yet exist, so that it remains unclear how dPDFs should best be defined (and how they evolve)". However some phenomenological aspects of QCD-evolution are known since long time (e.g Refs. [45,46]) and have been recently retaken [47,48,49] developing numerical codes able to solve the evolution equations. In addition also theoretical progresses have been reported (for example the demonstration that the exchange of Glauber gluons cancels for the considered observable, a step forward in the proof of QCD factorization for DPS [50]). In the following we develop a systematic numerical approach to the evolution of dPDFs, in Mellin space instead of coordinate space, restricting ourselves to the, so called, homogenous equation, a restriction we share with numerical solutions in coordinate space as applied in several contributions by Diehl and other coauthors (see Ref. [51] and reference therein).
If we assume equal renormalization scales Q 1 and Q 2 for the two partons (i.e. Q 1 = Q 2 = Q), the LO evolution equation for the unpolarized double parton distributions F j 1 j 2 (x 1 , x 2 ; Q 2 ) then reads (see Ref. [52]) The convolution integrals appearing in Eq. (49) have the same structure of the integrals appearing in the evolution of the single parton distributions, namely the renormalization group equation (RGE). In order to solve evolution equations, one can perform a M ellin-transformation of Eqs. (49), in particular for the first two terms where and the θ(1 − x 1 − x 2 ) appearing in the definition of the moments Eq. (51) is a direct consequence of the limit of integration in Eq. (49) and the momentum conservation. P ij are the evolution kernels or splitting functions. They are calculated perturbatively as a series expansion in a s (Q 2 ) = α s (Q 2 )/(4π): and m = 0 indicates the Leading-Order contribution.
(Expressions for P (0) ij in the context of dPDFs can be found (e.g.) in Appendix A of Ref. [44]).

A.1 dPDF (flavor) decomposition and evolution
In order to solve Eqs. (50) one has to combine the flavor indices in a way consistent with evolution, in particular one has to identify combinations evolving as Singlet and Non-Singlet. The combinations depend on the order of the evolution. At LO and NLO a useful transformation is the following Specifically, in the case of dPDFs F ij , the same argument holds for indices i, j combined in such a way to produce T 3 , T 8 and V i structures. Consequently, in addition to will evolve following the simple Non-Singlet rules. Just to give an example, we will discuss, in the next Section, the evolution of the dPDF Neglecting the inhomogeneous term, the solution of Eq. (50), for the Mellin moments of combination Eq. (57) is: (compare also the definition Eq. (51)). The M ellin-inversion completes the solution in x-space: (1−n 1 ) 1 x (1−n 2 ) 2 M n 1 n 2 VuT 3 (Q 2 ) . (59) The procedure described for the example F VuT 3 , is valid for each F ij combination ((i, j) = V i , T 3 , T 8 ).

A.2 Examples of Flavor decomposition
A. Flavor decomposition at the generic scale Q 2 and the inverse F ΣΣ = F u + u + + F d + d + + F s + s + + + (F u + d + + F d + u + ) + (F u + s + + F s + u + ) + + (F d + s + + F s + d + ) ; (68) These relations are generally valid, not only at the specific Q 2 0 . B. Reduction at the µ 2 0 scale At the lowest scale one has: F s + s + = 0 ; (79) and the inverse These relations are valid when the contributions at µ 2 0 reduce to valence contributions only.