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é 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).

JHEP10(2016)063
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", σ eff , 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 σ eff is constant w.r.t. the center-of-mass energy of the collision. In ref. [23] we have presented a predictive study of σ eff , 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 σ eff 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 JHEP10(2016)063 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][32][33], 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 [34,35], read (see, e.g., [12,13]) 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. [36]). 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, Instant-Form The two descriptions are related by Melosh rotations [37]. Following our previous developments (e.g. refs. [24,[31][32][33]) 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 [38]. 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 JHEP10(2016)063 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 The internal wave function is an eigenstate of the baryon mass operator and where the interaction term V must be independent on the total momentumP 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. (2.1) to the first (valence) Fock components and specializing the result to the u quarks as an example, one has (λ 1 , The spin projector values are determined by the Melosh rotationsD i :

JHEP10(2016)063
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 figure 1 the numerical results of the eqs. (2.8)-(2.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 (cf. eq. (2.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 [39]) as indicated by analogous Leading − Order calculations (see, e.g., refs. [24,27,28]). That scale is clearly indicated in the resulting expressions (2.11) and (2.12) and in both panels of figure 1.

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. [40], 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. [34,35]. The authors of ref. [40] 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 [41]), while eq. (2.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. [40] 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. (2.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. (2.17) does not factorize into separate contributions from each of the two partons, a and b (the values of the parameters in eq. (2.17) can be found in ref. [40]). 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., [42] 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. (2.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 (cf. eq. (2.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 (2.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][32][33] some years ago, one can check directly the accuracy of eq. (2.18). 1 Because of the natural normalization of the expression eq. (2.11): 19) and the normalization of the H u V GPDs the comparison holds for x 1 x 1 can have some validity in the restricted regions x 1 + x 2 < 1, the range where the dPDFs do not vanish. For x 1 + x 2 > 1 the dPDFs must vanish while the single parton responses H and E do not.
A detailed comparison is shown in figure 3 and figure 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 figure 3), but it extends to the kinematical region up to k 2 ⊥ = 0.5 GeV 2 (see figure 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.

Scale impact on correlations
3.1 Analysis of the approach of ref. [40] Let us first analyze the scale dependence of the correlations introduced at Q 2 0 within the assumptions eqs. (2.15) and (2.16), as proposed by Diehl et al.. in ref. [40]. 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. [40]. As it is done also in ref. [40], 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 [35,43]; its analysis is beyond the scope of the present paper.
In figure 5 we show the ratios eqs. (3.1) for different quark and gluon combinations: 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) 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 5. Effects of evolution on correlations according to the scheme of ref. [40]. Upper panel:   Figure 6. Effects of evolution on correlations for the LF-Hypercentral approach. Upper panel with the same notations. The gluon distribution F gg vanishes identically at Q 2 = µ 2 0 (cf. eqs. (3.3)).

JHEP10(2016)063
In particular the upper panel of figure 5 shows the valence components of . For those distributions, only the non-singlet evolution is relevant. In the case of , the singlet components are contributing in a substantial way; F gg (lower panel) is purely singlet. The distributions are defined as functions of the distance | y|, the Fourier transform would give the distributions as functions of k ⊥ without adding more information. The choice to show the y 2 -dependence makes easier the comparison of the results shown in figure 5 with the calculation of Diehl at al. as illustrated in their figures 1(a), 1(c), and 1(e).

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 In figure 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 figure 6, following the same notations and criteria of figure 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.

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 (non-perturbative) meson degrees of JHEP10(2016)063 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][32][33]44]).
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. (2.8) and (2.11), while the non-perturbative sea and gluons components are evaluated by means of a factorized approximation of the kind discussed in section 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: 2) The pure valence (and dominant) term, the expression (4.1) 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. (4.2) and (4.2), one can assume factorized forms (see e.g. ref. [3]). The complete (approximate) expression for F uu becomes: 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 ( The exponent n has to be fixed phenomenologically, as seen in section 2.3.1 in the case of the model of ref. [40] 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.

JHEP10(2016)063
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 give three examples: i) the identification of the exponent n to fix the correlating function eq.

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 LFapproach, 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. (4.5), 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 section 2.3, therefore one cannot expect eq. (4.7) 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. [45] and references therein). More recent arguments (see, e.g., refs. [40], 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. (4.7). 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. [41]
The advantages of the approximation eq. x 1 x 1 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. [41]). 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. (4.6) 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 x 1 x 1 which contains the additional sea contributions fromū (u = u V +ū, see text) (continuous lines). x 1 complete Mellin procedure we are proposing for both Singlet and Non-Singlet sectors is illustrated in some detail. The result is shown in figure 7 for four selected values of x 2 .
In figure 8 and figure 9 we show the complete set of dPDFs involving the u-quark at the starting scale Q 2 0 . In particular, the combination figure 8 giving explicit evidence to the contribution due tō u quarks. The largest effects of the sea component are clearly evident for smallest values of the momentum fraction. In figure 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 (cf. figure 7).

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 (cf. for example, figures 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. (4.5) and (4.7), valid at k ⊥ = 0. For instance, eq. (4.7) becomes and simple choices are (cf. eq. (2.16)), Within scenario A of eq. (4.11), 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. (4.12) the exponent depends on x 1 , x 2 , similarly to eq. (2.16). The knowledge of u V u V (x 1 , x 2 , k ⊥ , µ 2 0 ) LF from eq. (2.11) and of the sPDF u V (x 1 , µ 2 0 ) LF with the additional information n = 0.2 from the analysis of section 4.1.1, can be used to optimize the fit eq. (4.10). The results of the optimization procedure are shown in figure 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. (4.12)). 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. (4.10) and scenario B.

Sea and gluon contribution at k ⊥ > 0
The valence-valence dPDFs do not allow for a simple factorization. The approximation proposed in eq. (4.10) is therefore quite poor in that case, as figure 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 (cf. figure 1 and eqs. (2.8), (2.11)) and, using the best factorization scheme (n = 0.2 and b A = b B = 0.6 GeV −1 in eqs. (4.11), (4.12)) 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 section 4 and the procedure described in section 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 figure 11 the distribution F uu (x 1 , x 2 , k ⊥ , Q 2 0 ) at  the scale Q 2 0 which generalizes eq. (4.5) and figure 8, and whose resulting expression reads (from now on, we will discuss the scenario B of eq. (4.12) only): where u V u V (x 1 , x 2 , k ⊥ , Q 2 0 ) is obtained evolving at the scale of the MSTW parametrization, Q 2 0 = 1.0 GeV 2 , the LF result (2.8), (2.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 [41].

JHEP10(2016)063
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, (figures 12 and 13), we compare the results obtained evolving directly from the lowest scale µ 2 0 where only valence-valence dPDFs are present (cfr. figure 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 (cf. figure 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 figures. The Singlet components parametrized by means of the factorization procedure appear to play a relevant role at low x, where the dPDFs can be more easily studies by means of proton-proton 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.

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, nonsinglet and singlet, LO QCD evolution to the very high Q 2 scales typical of pp scattering at the LHC. As in ref. [40], 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. For the seek of clarity let us stress that all calculations performed with the initial scale µ 0 are associated to the LF model where valence quarks are the only non perturbative degrees of freedom while those performed starting from Q 0 are related to the model where non perturbative gluons and sea quarks have been taken into account in the analysis, see eq. (4.14).

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 within our LF scheme for both sPDFs and dPDFs, starting from the hadronic scale µ 2 0 , where only valence degrees of freedom are present. Results for the evolution which includes gluons and sea degrees of freedom within the factorization scheme of section 4 will be discussed in the next section 5.2.
It is clear that the ratio (5.1) 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 − Singletindex, 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. (5.2) for the flavor combinations gg, u V u V , u Vū ,ūū, u V g are shown in figures 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. figure 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. figure 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 JHEP10(2016)063 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 ratios 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    figure 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 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. Numerator and denominator are evolved by means of dPDF evolution and single parton evolution, respectively. The starting scales differ for the two curves: the low momentum one µ 2 0 (continuous) and the larger Q 2 0 (dashed), where gluons and sea quark, generated non perturbatively according to eq. (4.14), are taken into account. The differences are artifacts due to numerical uncertainties. Right panel: the same of the left panel but for the ratio gg . 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.

JHEP10(2016)063
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.

A Perturbative evolution of dPDFs in Mellin space
Following Diehl and Kasemets [46] 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. [47,48]) and have been recently retaken [49][50][51] 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 [52]).
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. [53] 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. [54]) The convolution integrals appearing in eq. (A.1) have the same structure of the integrals appearing in the evolution of the single parton distributions, namely the renormalization JHEP10(2016)063 group equation (RGE). In order to solve evolution equations, one can perform a M ellintransformation of eqs. (A.1), in particular for the first two terms ij in the context of dPDFs can be found (e.g.) in appendix A of ref. [46]).

A.1 dPDF (flavor) decomposition and evolution
In order to solve eqs. (A.2) 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 and similar combination if one includes heavier quarks (e.g. ref. [55]  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 The procedure described for the example F VuT 3 , is valid for each F ij combination ((i, j) = V i , T 3 , T 8 ).
On the other hand, the double-distributions containing gluons and Σ evolve mixing the two and each index must be evolved in the appropriate way. For example:
B. Reduction at the µ 2 0 scale At the lowest scale one has: These relations are valid when the contributions at µ 2 0 reduce to valence contributions only.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.