Double parton scattering and the proton transverse structure at the LHC

We consider double parton distribution functions (dPDFs), essential quantities in double parton scattering (DPS) studies, which encode novel non perturbative insight on the partonic proton structure. We develop the formalism to extract this information from dPDFs and present results by using constituent quark model dPDF calculations within the Light-Front approach, focusing on radiatively generated gluon distributions. Moreover, we generalize the relation between the mean transverse partonic distance between two active partons in a DPS process and the so called effective cross section to include correlations effects and the so called 2v1 mechanism contribution. Finally we investigate the impact of relativistic effects in digluon distributions providing an economic parametrization of their contributions which can be used to improve phenomenological dPDF ansatz.


Introduction
A proper description of the event structure in hadronic collisions requires the inclusion of the so called multiple parton interactions (MPI) which affect both the multiplicity and topology of the hadronic final state [1,2]. The Large Hadron Collider operation renewed the interest in MPI given the continuous demand for an increasingly detailed description of the hadronic final state which is crucial in many New Physics searches. In this rapidly evolving context, this type of studies have also received attention for their own sake since they might be sensitive to double partonic correlations in the colliding hadrons, see recent review in Ref. [3]. The simplest MPI process is the double parton scattering (DPS) [5,4]. In such a process, a large momentum transfer is involved in both scatterings which enables the use of perturbative techniques to calculate the corresponding cross section. The latter depends on two-body quantities, the so called double Parton Distribution Functions (dPDFs), which are interpreted as the number densities of parton pair at a given transverse distance, b ⊥ , and carrying longitudinal momentum fractions (x 1 , x 2 ) of the parent proton [1,6,7]. Double PDFs are not calculable from first principles, a feature shared with usual PDFs and other quantities in QCD. Moreover, due to their dependence upon the partonic interdistance [7], they contain information on the hadronic structure complementary to those obtained from one-body distributions such as generalized parton distribution functions (GPDs) and transverse momentum dependent PDFs (TMDs). Unfortunately the DPS cross section is obtained by integrating dPDFs over b ⊥ so that such a dependence is not directly measurable [1]. In this scenario, hadronic models have been used to obtain basic information on dPDFs and to gauge the phenomenological impact of longitudinal and transverse correlations, see Refs. [8,10,9,11,12,13], along with spin correlations [11,14,15,16,17]. We mention that quantities related to dPDFs, and encoding double parton correlations, have been recently calculated for pion by means of Lattice techniques [18]. Despite the wealth of information encoded in dPDFs, present experimental knowledge on DPS cross section is accumulated, up to now, into the so called effective cross section, σ ef f , for recent results see e.g.Refs. [19,20,21,22,23,24]. The latter is defined through the ratio of the product of two single parton scattering cross sections to the DPS cross section with the same final states. In the present paper, we continue the investigation of the relationship between σ ef f and the mean interpartonic distance pursued in Ref. [25]. We study its modification induced by including the so called splitting 2v1 term contribution in DPS processes [32,29,31,28,30,26,27,33,34,35] and, separately, the effects of longitudinal correlations in dPDFs. Numerical estimates will be shown and discussed in the kinematics of DPS processes initiated by digluon distributions, see e.g. Refs. [37,38,40,39,36] on this topic, which are perhaps the most interesting objects in the DPS context. Such distributions are radiatively generated by pQCD [33,34,6,42,32,43,27,28,44,26,45,41,35] evolution starting from dPDF model calculations at the initial hadronic scale. Thanks to this procedure, one can avoid additional modelization of the non perturbative see quarks and gluons contributions, such those proposed in Ref. [46]. In such a radiative scheme and given the structure of evolution equations, digluon interdistance will follow the pattern of valence quarks interdistance obtained from the underlying hadronic model used for the dPDF calculations. DPS measurements, sensitive to gluon initiated processes, will then provide test to our analyses. In the last part of the paper, we focus on relativistic effects in dPDF model calculations, in the relevant kinematic conditions of collider experiments, and already addressed in Ref. [12] for valence quarks at the hadronic scale. The aim of this analysis is to offer insight to unfactorized ansatz for dPDFs and to give an economical parametrization of these effects which could be implemented in phenomenological analyses of dPDFs.
The paper is organized as follows. In Sec. 2 we will discuss processes and corresponding kinematic conditions which we will focus upon in this analysis. In Sec. 3 we describe the formalism which allows to obtain physical information on the proton structure from dPDFs. In Sec. 4 we introduce the so called σ ef f , relevant quantity the DPS analyses and show how the latter is related to the geometrical properties of the proton. In Sec. 5 we discuss relativistic effects in dPDF calculations and in Sec. 6 we provide a parametrization for such effects. We present our conclusions in Sec. 7.

Analysis strategy and calculation details
In the present analysis we will focus on the digluon distributions and therefore we consider DPS prototype gluon-dominated-processes: where each final state particle is produced in a distinct parton-parton scattering. Double J/Ψ production has been already measured both at Tevatron and LHC [47,20,23,24]. Double Higgs production via DPS has been studied in the literature [30], but not yet measured given its rather low cross section. We mention here that it would be also interesting to consider the mixed process pp → HJ/ΨX with final state produced via DPS which, to the best of our knowledge, has never been considered in the literature. We also mention that, interesting information could be gained from the comparison of the double J/Ψ production with the double open charm one in the same kinematics. The combined measurements of these DPS processes, among many others with a less pure gluon initial state but larger cross section, give a wide coverage of digluon distribution both in hard scale and fractional momenta. We define the partonic subprocess in the two scatterings as where p's and k's are the relevant parton flavour and momenta, respectively. Since heavy particles appearing in Eq. (1) are produced by partonic annihilation in lowest order of perturbation theory, the fractional momenta of the incident gluons can be reconstructed from the mass m, transverse momentum k T and rapidity y of final state particles as In our calculations we set the centre-of-mass energy to its nominal value at the LHC, √ s=13 TeV and consider two rapidity region: the central one covered by ATLAS and CMS, |y| < 1.2 and the forward one covered by LHCb, 2 < y < 4.5. Neglecting transverse momentum, J/Ψ production gives access to fractional momenta in the range 10 −6 x 10 −2 while for Higgs production 10 −4 x 1. The factorization scale in each process is set equal to the mass of the particle, either the J/Ψ or Higgs boson, produced in the final state, µ F,A = m A and µ F,B = m B with m J/Ψ = 2m c . The differential DPS cross section, assuming that the two hard scatterings can be factorized [6,28,41,49,48], involves dPDFs through an integral over the transverse partonic distance b ⊥ and reads [1,6]: In Eq. (4) dσ are the differential partonic cross sections for processes with final state A or B respectively and the symmetry factor reads m = 1 if A = B and m = 2 otherwise. Double PDFs appearing in Eq. (4), are multidimensional distributions addressing non perturbative features of the proton structure and are therefore complicated to model. Some guidance in building appropriate initial conditions is offered by physical intuition at small x [29,6,50,44,28,41] and by sum rules [29,38,37,51]. Nevertheless, a large freedom is left in the gluon transverse spectrum, which is perhaps one of the most intriguing aspect for hadronic studies. In order to investigate some of these features, in the present paper we make use of dPDF calculations within constituent quark models (CQMs), e.g. Refs. [8,9,10,11]. Following the line of Ref. [12], we have adopted the hypercentral quark model (HP), in its relativistic version [52] and, in order to highlight model independent effects on dPDFs, the harmonic oscillator model (HO) [53]. In particular, for the latter, we considered the version described in Ref. [12], where the model parameter α, representing the width of the Gaussian, is set to be α 2 = 25 fm −2 in order to mimic a relativistic structure. These models differ in many dynamical aspects and offer a parametrization of the only non-vanishing valence-valance dPDF at the hadronic scale Q 0 . All other distributions are then radiatively obtained at higher scales by performing pQCD evolution in its homogeneous form, which is appropriate at fixed b ⊥ [44,6,54]. Since the initial scale of the model Q 0 , is located in the infrared region, dPDFs show a large sensitivity to its precise value. In order to overcome this hitch, Q 0 has been fixed according to the procedure outlined in Ref. [55] and resulting to be Q 2 0 = 0.26 GeV 2 . In order then to reduce the impact of this choice, we will often consider appropriate ratios which reduce, and in many cases almost cancel, this dependence. We mention here that single parton distributions have been obtained from the same HO and HP models and evolved from the valence initial condition to higher scales. In this case we choose the same value of Q 2 0 used for dPDFs. This feature is particularly relevant for the calculation of the effective cross section which we will be introduced in the next Section. It is worth to remark that the digluon distribution is likely to have a non perturbative contribution at Q 0 . In the present work, we make use of a pure radiative evolution scheme, and therefore it is our precise choice to neglect such an additional input in order to avoid any ambiguities arising from its modelization. Therefore all the presented results must then be interpreted accordingly.

Proton transverse structure from dPDFs
In this section we present the formalism necessary to extract physical information on the proton structure from dPDFs, i.e. the mean partonic distance between two partons in the transverse plane and its correlations with the longitudinal parton momenta. In fact, since dPDFs represent the number density of two parton with longitudinal momentum fractions x 1 and x 2 at a given transverse distance b ⊥ , see Refs. [1], they provide a new tool to access the 3D structure of the proton, complementary to that obtained from generalized parton distribution functions (GPDs). In particular, these two-body functions are sensitive to double parton correlations [11,46,12,13,44,8,9,3,10,45,18,14,35,36] that can not be accessed by means of one-body distributions such as GPDs. To this purpose we first introduce the effective form factor (eff) [56,25] as the first moment of the dPDFs: where i and j are flavor indexes and N ij factor is the dPDFs normalization, evaluated at k ⊥ = 0, according to sum rules discussed in Refs. [29,51,54]. Here F ij (x 1 , x 2 , k ⊥ ) is the Fourier transform (FT) of the dPDF F ij (x 1 , x 2 , b ⊥ ) in coordinate space and k ⊥ represents a transverse momentum imbalance between two partons in the amplitude and its conjugate [28]. Therefore, as for GPDs, F ij (x 1 , x 2 , k ⊥ ) does not admit a probabilistic interpretation in k ⊥ -space, which holds instead in b ⊥ -space. Thanks to the definition in Eq. (5), the eff represents the Fourier Transform of the probability of finding two partons at a given transverse distance. Moreover, as discussed in Ref. [25], the formal definition Eq. (5) allows to connect the eff in terms of the proton wave function. Such a condition suggests to deduce basic properties of the eff by the comparison with the standard electro-magnetic proton form factor. Thanks to this procedure, in Ref. [25], unknown information on the proton structure have been related to σ ef f . In the present section we extend the analysis of Ref. [25] by directly starting from dPDFs, being these distributions those entering the DPS cross section, see Eq. (4), rather than the eff. The present approach opens to the possibility of accessing new information on the proton structure such as correlations between the longitudinal momentum fractions of two partons and their mean partonic distance. The latter, for a pair of partons with flavour i and j and fractional momenta x 1 and x 2 , is defined as where Q 2 is a generic hard momentum scale at which dPDFs are evaluated and we have From Eq. (6), the mean partonic distance between partons is obtained as follows: where the factor 2 is due to the symmetry in the exchange of particles. As for the standard electro-magnetic nucleon form factor, such a relation can be equivalently obtained from dPDFs in momentum space. Since it follows that

Kinematics
HO model HP model From the above relation, Eq. (6) can be equivalently written in terms of dPDFs in momentum space, in analogy with the standard electro-magnetic form factor: Given the really limited knowledge on dPDFs, these quantities, together with the mean interpartonic distance, can all be evaluated via hadronic models. In Fig. 1 we present the digluon dPDFs, evaluated within the HO (left panel) and HP (right panel) models at Q 2 = m 2 H in b ⊥ -space. Since we consider unpolarized partons in an unpolarized proton, circular symmetry in b ⊥ is obtained, as apparent from the plot. Furthermore, the shape of the distributions are qualitatively similar to those shown in Ref. [12], where valence quark dPDFs have been evaluated within the same models but at the hadronic scale. By using these quantities, we have also evaluated the mean gluonic distance via Eqs. (6,10). The results, reported in Tab. 1, show that partonic correlations induce a dependence of the mean partonic distance upon the longitudinal momentum fractions carried by two partons. We recall that if correlations between x i and k ⊥ were absent, as in the non relativistic limit of dPDFs evaluated within the HO model (see Ref. [9]), the mean partonic distance does not depend on x's and reads b 2 = 0.283 fm. This discussion, however, is rather academic since the present accuracy of DPS measurements is far from being sensitive to this kind of effects. Nevertheless we have shown in Ref. [25] that physical information on the proton structure can still be obtained directly from σ ef f , a quantity which is often used in experimental analyses. In next sections we review the formalism that allows one to relate σ ef f to b 2 , and generalize it to more complicated cases.

Transverse Proton structure from effective cross section
Since dPDFs, the main non-perturbative ingredients appearing in Eq. (4), are basically unknown, a fully factorized ansatz is frequently assumed in phenomenological and experimental analyses: where q i (x) are ordinary single PDFs and f (k ⊥ ) is the effective form factor defined in Eq. (5). Usually, in such a simplified approach, f (k ⊥ ) does not depend on the parton flavors [5]. These assumptions allows to rewrite the DPS cross section as [7,57] being dσ SP S the single parton scattering cross sections with final state A(B). In this scenario σ ef f simply reads: where the last expression follows from rotational invariance. Eq. (12) shows that, in such an approximations, σ ef f enters the DPS cross section as an overall normalization factor. In Ref. [25], we have shown that, by using the formal definition of the eff in Eq. (5), appearing in Eq. (13), one can relate σ ef f to the mean partonic distance of two partons active in a DPS process. We will briefly review this procedure in the following. Since eff is the FT of the probability distribution of finding two parton at a given transverse distance in a confined quantomechanical system, two asymptotic conditions on eff can be imposed: By using the above relations and integration by parts, one obtains the following identities: with m ≥ 0, and with s ≥ 0. Since, as already mentioned, the eff is the FT of a probability distribution, it can be expanded as follows [25]: where the P J0 n are the expansion coefficients of the Bessel function and b 2n are weighted moments containing the dynamical information on the partonic proton structure. Combining these identities one finds see Ref. [25]. Then, by using variance properties, the second term on the right hand side turns out to be positive definite, so that one obtains A maximum to the mean partonic distance between the two partons can be found by solving the following integral inequality being N an unknown number and the function d 2 (k ⊥ ) defined by and normalised to In order to solve the inequality Eq. (20), one needs an additional assumption on the behaviour of the eff at large k ⊥ . Since f (k ⊥ ) can be considered a double form factor, it is reasonable to assume that it falls to zero at large k ⊥ at least as fast as a standard one-body form factor, i.e. faster then 1/k ⊥ , see Ref. [25], for additional details on this point. In such a case one finds a range of values for the arbitrary number N , 1/2 ≤ N ≤ 1, for which the following sufficient condition holds: Combining all these results one obtains an allowed range for the partonic interdistance: which is the main result of Ref. [25]. The above relation has been checked within all models of the eff in the literature. Let us remark that in order to make contact with experimental extraction of σ ef f , this result has been obtained under the approximation of Eq. (13). Thanks to this feature, data on σ ef f have been converted in the range of b 2 [25]. In the following Sections we will describe how σ ef f can be generalized to include partonic perturbative and non perturbative correlations, thus breaking the factorized ansatz in Eq. (13).

Generalization to 2v1 case
As discussed in Ref. [32,29,17,31,28,30,58,35], the DPS cross section might receive a contribution from the so called 2v1 mechanism. In this case, the parton pair active in the processes is perturbatively produced from the splitting of a single parton, see the right panel of Fig. 2. Given the large gluon flux at LHC energies, such a contribution can be non-negligible [35], e.g., for double quarkonia and/or Higgs production [30] with respect to the standard 2v2 mechanism shown in the left panel of Fig. 2. This contribution breaks the simple ansatz in Eq. (13) and it is of pure perturbative origin. Its presence in dPDF evolution equation and in DPS cross sections has been long debated, especially due to the double counting problem, see recent results on this issue in Refs. [31,54,49,58]. In this Section we consider the formalism developed in Refs. [30], where the σ ef f definition is generalized to include the contribution of the perturbative splitting of one parton of the ladder in to two, for example g → gg, which subsequently enter the DPS process. As discussed in Ref. [30], one can decompose the total DPS cross section in terms of the two leading 2v2 and 2v1 contributions as follows where here Ω 2v2 and Ω 2v1 represent the DPS cross sections calculated with longitudinal double PDFs for both mechanisms, and weighted by the corresponding σ ef f . In particular the Ω 2v1 term is calculated with dPDFs whose initial condition is given by the splitting term alone at the initial scale [30]. It is important to remark that in the 2v1 case, the parton pair generated by the splitting does appear only at b ⊥ = 0 [30,43,58,6]. Since in experimental analyses the 2v1 mechanism is neglected, it is usually assumed that σ DP S = Ω 2v2 /σ ef f . Nevertheless, it is possible to incorporate the 2v1 contribution in σ DP S by using the following generalization of σ ef f [30]: Under the assumption that the longitudinal dependence of dPDFs factorizes from the transverse one, the effective cross sections for the two mechanisms read [30]: 1 where it is worth noticing that both the above expression depends on the same effective form factor, f (k ⊥ ), so that they are not independent quantities. The first equation is the standard one, see Eq. (13). The second one reflects the perturbative production of the couple of partons, occurring at zero relative distance in transverse plane. In terms of the present notation, the main result of Ref. [25] reads: Figure 2: Diagrammatic representation of the two contributions to a DPS process: the so called 2v2 mechanisms is shown in the left panel and the 2v1 mechanism in the right panel. Small grey blobs represent the hard scattering processes.
where, by following Ref. [30], σ ef f,2v2 represents the usual definition of σ ef f if only the 2v2 mechanism is considered, see Eq. (27). In the case where also the 2v1 mechanism is included in the analysis, in order to relate b 2 to the experimentally extracted σ ef f Eq. (26), we need first to find a relation between the mean partonic distance and σ ef f,2v1 , defined in Eq. (28) and appearing in the full definition of σ ef f in Eq. (26). To this aim we use the identities in Eqs. (15,16) with m = 1 to obtain: Due to variance properties [25], the overall sign of the second term of the above equation is positive and consequently: 1 Furthermore, similarly to the case of Ref. [25], in order to estimate a reasonable maximum, one needs solve the following inequality: whereN is an arbitrary unknown number. Under the additional assumption that the eff falls to zero at large k ⊥ at least as fast as standard proton form factors, one finds the desired condition: Linking Eq. (31) and Eq. (33), the following result is found: Combining Eq. (29) and Eq. (34) in Eq. (26) one obtains the final inequality: where here we have defined the ratio r v = Ω 2v1 /Ω 2v2 , with r v ≥ 0. Let us remark that, in principle, the ratio r v could depend on the rapidities of particles produced in the final state and hence on parton fractional momenta in the initial state [30]. Such a dependence is not explored in the present analysis. The difference between the maximum and the minimum in Eq. (35) gives an estimate of the theoretical error on the transverse distance of the two active partons: The main effect of the inclusion of the 2v1 mechanism, is not only to shift the b 2 range towards higher values, but also to increase its theoretical error with respect to the case where r v = 0. In particular, the last comparison, between the r v = 0 and r v = 0 cases, makes sense only if σ ef f is equal in both scenarios. In principle, as observed in Refs. [28,30], in order to observe σ ef f ∼ 15 mb, one should expect σ ef f,2v2 ∼ 30 mb.
In general, if r v =1, from Eq. (26) one gets σ ef f ≤ σ ef f,2v2 . We find interesting to check the validity of Eq. (35) by using two phenomenological models for eff, such as those described in Refs. [28,30,59]. The first one is Gaussian eff of the type: In this case the mean partonic distance can be obtained in term of the width parameter a as: so that, according to Eqs. (27,28), By using the above expressions in Eq. (26), one gets the following result: which is included in the range Eq. (35). As a second example we consider an eff which is the square of the gluon form factor [59], i.e.: with the parameter m g has been fixed by fitting HERA data, i.e. m 2 g ∼ 1.1 GeV 2 [59]. In this case one obtains: and, according to Eqs. (27,28), By using the above expressions in Eq. (26), one gets the following result: which again lies in the range indicated in Eq. (35). Such a generalization of the inequality in Eq. (24) is however process dependent, in fact, as discussed in Refs. [30,28], r v is related to the kinematic conditions and to the type of the considered DPS process. Without a precise knowledge on r v , a determination of the range of the allowed mean partonic transverse distance is therefore prevented. Nonetheless, we note that, as already discussed in Ref. [30], 0 ≤ r v ≤ 1. In fact, r v is a ratio of cross sections which is positive definite. Furthermore r v ≤ 1 since the 2v1 mechanism is subdominant in the pQCD evolution of dPDFs w.r.t. the 2v2 one, being the former proportional to the gluon density and the latter proportional to its square. Thanks to these features, for r v = 0 one finds the absolute minimum in Eq. (35) while for r v = 1 one finds the absolute maximum: This result allows one to obtain information on the interpartonic distance of two active partons in a DPS process without knowing details on the relative size of the two mechanisms, 2v1 and 2v2, i.e. the exact knowledge of  (24). Vertical lines represent the full range of allowed partonic distances between two partons if the 2v2 and the 2v1 mechanism effects on the total σ ef f are not disentangled. The shadow between lines represents the additional theoretical error w.r.t. the case where only the 2v2 mechanism is considered, i.e.Eq. (24). The areas outside lines represent the exclusion region of the allowed transverse distance between two partons active in a DPS process. The red line stands for twice the transverse proton radius.

Ref.
Process  r v . This, of course, comes at the expense of an increased theoretical error. In order to quantify such an effect, we have plotted in Fig. 3  In addition, in Table 2, we report the interpartonic distances, calculated according to Eq. (35), for σ ef f values extracted from a selection of experimental analyses. It should be noted that, in all cases, b 2 < 2R ⊥ = 1.42 fm, where R ⊥ ∼ 0.71 fm is the transverse electro-magnetic proton radius. We close this Section by observing that Eq. (45) can be inverted to give: In such a form, given the value of b 2 associated to a particular f (k ⊥ ), the inequality (46) predicts the expected range in σ ef f associated to that specific model. Most importantly, Eq. (46) shows that, given an eff, characterized by b 2 , the σ ef f value does depend on the relative size of the 2v1 contribution. In particular, if r v is significantly larger than zero, the corresponding σ ef f will be lower than the one obtained if the 2v2 mechanism alone were considered (r v = 0).
Kinematics HO model HP model

Generalization to the unfactorized ansatz
As shown in several constituent quark models calculations of dPDFs, double parton correlations may survive at high momentum scales, thus breaking the factorized ansatz in Eq. (11), see Refs. [8,9,11,46,55,12]. Therefore in this Section we investigate how Eq. (24) is generalized to the case in which the factorized ansatz is not assumed thus allowing the presence of longitudinal and mixed longitudinal-transverse partonic correlations. We therefore consider an unapproximated scenario in which σ ef f depends on the longitudinal momentum fractions of the active partons, as suggested in Refs. [56,60]. Within this improved framework, the relationship between σ ef f and the mean partonic distance will be sensitive to x 1 − x 2 correlations. For this purpose we consider the next-to-simple generalization of the results presented in Section 3, namely we consider non-factorizable dPDFs.
In the zero rapidity case, i.e. x i = x i , σ ef f can be written as follows [56]: where, for clarity, we suppress flavor indeces, F (x) is the standard PDF. As discussed in Ref. [25], similarly to the case of GPDs, whose transverse momentum dependence is dominated by the associated form factors, in this analysis we consider that the k ⊥ dependence of dPDFs is dominated by that of the eff. Such a feature is also discussed in Sec. 3. Thanks to the latter condition, identities in Eqs. (15,16) can be generalized as follows: In addition, the expansion in Eq. (17) can be also used in the dPDF case and it reads: By using the generalized identities in Eqs. (48,49,50) and general arguments on the k ⊥ dependence, previously discussed and addressed in Ref. [25], one can relate σ ef f to the mean partonic distance. The inequality (24) is then generalized to: and additionally depends on the ratio r ij defined by Such a ratio encode longitudinal correlations in the proton structure, and therefore so does b 2 x1,x2 . Such a generalized inequality can be tested against dPDFs obtained via model calculations. To this aim we have evaluated the terms appearing in Eq. (51), i.e. σ ef f (x 1 , x 2 ), r ij (x 1 , x 2 ) and b 2 x1,x2 within CQMs. Since we are interested to kinematic regions close those experimentally accessed, we have calculated the above quantities In addition, in order to assess the hadronic model dependence of the results, Eq. (51) has been calculated with digluon distribution obtained within two different CQMs. The results are reported Table 3 and, as one may notice, the inequality Eq. (51) is verified in all kinematic conditions. One should also notice that, at variance with the case where the factorization ansatz in Eq. (11) is assumed, in this new scenario the effects of correlations in dPDFs, embodied in the r gg (x 1 , x 2 ) factor, play a crucial role in verifying the identity. This generalized inequality effectively allows one to estimate the impact of double parton correlations on the range of allowed parton transverse distances. Before closing this Section, we find interesting to discuss the role of r ij in the case in which double and single PDFs are evaluated within CQM calculations. When dPDFs are built upon the factorized ansatz in Eq. (11) such a ratio reduces, by construction, to unity at small x. However, when both single and double PDFs are calculated via CQMs, a spurious effect related to their respective normalization arises. Let us consider single and double PDFs constructed from a given hadronic model at Q 2 0 . At that scale, the model predicts only the valence contribution: for single PDF it gives F uv (x, Q 2 0 ) which is normalized to give the number of u v -quark in the proton: Similarly, model calculations return the double distributions for valence quarks F uv,uv (x 1 , x 2 , , k ⊥ = 0, Q 2 0 ), which is normalized to [29] 1 0 dx 1 Both single and double PDFs are then evolved to higher scales according to their respective evolution equations.
If, at such low scale a double parton distribution is constructed as the product of valence quark single PDFs, which indeed is a very bad approximation [9,11,46,29], one gets the following normalization The above normalization is determined by the dPDF support condition, i.e. x 1 + x 2 ≤ 1, and it weakly depends on the specific hadronic model. When single and double PDFs are then used to calculated, for istance, σ ef f , such a normalization mismatch is propagated via pQCD evolution to sea quarks and gluons giving a ratio in Eq. (52) substantially lower than one, especially in the small x region, where r ij is expected to reflect the factorization of the dPDFs into single PDFs, i.e r ij ∼ 1. We show in Fig. 4 the ratio r gg computed within the HO and HP models in different kinematics conditions at two different scales. From the plot is clear that at the smallest accessed x, the ratio approaches the value 0.6, which is indeed the ratio of the two above normalizations, i.e. 2/3.4 ∼ 0.6. In this case, r ij = 1, in the small x region, is therefore an artefact of the normalization and does not imply any type of correlation. This discussion is particularly relevant for the present investigation since both single and double distributions are obtained from hadronic models in which only valence quarks are present at the hadronic scale. On more general grounds such an effect is simply traced back to the assumption that double PDFs, constructed at any scale as a product of single PDFs, reflects the inherent normalization of the latter at Q 0 . The ratio in Eq. (52) displayed in Fig. 4 shows some interesting additional features. Its scale dependence, indicated by the shadow area between the various curves, increases as moving toward higher fractional momenta, i.e. approaching the boundary. In the very same limit r gg is a decreasing function of x, meaning that dPDFs fall on the boundary much faster than a product of single PDFs. Moreover, the dependence of the ratio on the adopted hadronic model is rather weak. A final remark on the interpretation of σ ef f is in order. On the theoretical side, σ ef f has a direct connection to the mixed transverse-longitudinal structure of the proton through its dependence on dPDFs. On the experimental side, this quantity merely represents a measure of the suppression of the DPS cross section with respect to the product of the SPS ones. Therefore, as suggested by the results outlined at the end of this Section, the numerical comparison between values of this quantity extracted from theoretical and experimental analyses requires care. In experimental analyses which extract σ ef f under the factorization assumptions for dPDFs, r ij is unity by construction, and σ ef f is almost completely associated to the transverse structure of the proton. In CQM calculations such a values does depend both on the transverse distance distribution of the two partons, the difference between dPDFs and PDFs product and spurious effects due to normalization. In such a case σ ef f predicted within constituent quark models should be rescaled by a factor r 2 ij in order to compare the physical information on the transverse structure of the proton to the one extracted in the experimental analyses.

Relativistic effects in dPDFs
In this Section we consider relativistic effects on dPDFs, already addressed in Ref. [12], and study their relevance when propagated at high momentum transfer in typical LHC kinematics, with a special emphasis on the gluon sector. For this purpose we follow the procedure described in Ref. [12] where dPDFs have been evaluated within the Light-Front (LF) formulation. Among the three forms of relativistic dynamics [61], the LF one has the maximum number of kinematical generators, such as LF boosts [61]. This feature makes the LF approach suitable to implement special relativity for strongly interacting systems [62,63,64] and therefore it has been extensively used to evaluate other kind of parton distributions [65,66,67,68]. We consider the dPDFs expression presented in Ref. [11], i.e.: where k i is the intrinsic three-momentum of the i parton whose flavor is determined by τ i , k ⊥ is the relative transverse momentum unbalance in the parton pair, Ψ is the proton canonical (instant form) wave function in momentum space and |S ⊗ F is a generic spin-flavor state. M 0 is the proton mass with constituent quarks treated as free particles and whose dependence on x i and k i⊥ is given by: being m i and x i the constituent quark mass and longitudinal momentum fraction carried by the i quark, respectively. Here, as in Ref. [11], we consider for simplicity a factorized dependence between the spin-flavor and the spatial part of the proton wave function. Thanks to the LF approach, momentum conservation is preserved, i.e. dPDFs vanish in the unphysical region x 1 + x 2 > 1. The canonical proton wave function appearing in Eq. (56) can be calculated within constituent quark models, see e.g. Refs. [9,11]. Nevertheless, the price for the use of the canonical proton wave function is the inclusion of boosts from the Light-Front centre of mass frame to the instant form one, i.e. the so called Melosh operators [69], which appear in the second line of Eq. (56) and are defined asD where σ x and σ y are Pauli sigma matrices. In particular, the Melosh operators allow to rotate Light-Front spin into the canonical one. We emphasize that for unpolarized PDFs, for which the initial proton state is equal to the final one in the light-cone correlator, the product of Meloshs reduce to the unity,D †D = 1. However, as shown in Ref. [11], in the case of dPDFs, for which in general k ⊥ = 0, Melosh operators contribute also in the case of unpolarized partons. In the present analysis, we are interested in (x i − k ⊥ ) correlations induced by Melosh operators on dPDFs. However, given the complicated structure of Eq. (56), it is non trivial to single out their effects, since they mix with the proton wave function. In order to determine to which extent their effects on dPDFs are independent of the chosen hadronic model, we compare dPDF calculations performed within the HO and the HP models and build appropriate ratios in order to highlight relativistic effects alone.
In Fig. 5 we present the double gluon distribution in coordinate space,F gg (x 1 , x 2 , b ⊥ , Q 2 = m 2 H ), evaluated in different configurations of x 1 and x 2 , including (black full lines) and neglecting (orange dashed lines) Melosh operators within different hadronic models, the HP in the upper panels and the HO in the lower ones. Results are consistent with those of Ref. [12], where only valence quark dPDFs was evaluated at the low hadronic scale HO model HP model  Table 4: Calculations of σ ef f in the relevant experimental rapidity range of the process pp → J/ΨJ/ΨX. Results are presented for digluon distribution evaluated at Q 2 = 4m 2 c and obtained within the HO and HP models, including and neglecting Melosh operators.
of the models. We observe in the plots that there exists a value of b 0 ∼ 1.5 GeV −1 , which slightly depends upon the kinematics and the hadronic model used in the calculations, such that the inclusion of Melosh operators strongly decrease dPDFs for b ⊥ < b 0 and slightly increase them for b ⊥ > b 0 . It is worth noticing that Melosh operators reduce to the identity for k ⊥ = 0, so that dPDFs with and without Melosh coincide in this limit. Since the latter condition corresponds to an integral of dPDFs over d 2 b ⊥ , it follows that dPDFs with and without Melosh are normalized to the same number. The digluon distributions in Fig. 5 show a marked dependence on the specific proton wave function built-in the CQMs. It is therefore instructive to present the ratio: where we indicate withF N M gg (x 1 , x 2 , b ⊥ , Q 2 ) the dPDFs in Eq. (56) evaluated neglecting Melosh operators. The ratio in Eq. (59) is shown in Fig. 6 with calculations performed within the HP model at the final scales Q 2 = m 2 H (full lines) and Q 2 = 4m 2 c (dot-dashed lines), and the HO model at the final scales Q 2 = m 2 H (dashed lines) and Q 2 = 4m 2 c (dotted lines) in three configurations of x 1 and x 2 . As one can see, up to b ⊥ < b 0 , Melosh operators induce a sizeable reduction of dPDFs which is almost a kinematical and scale independent effect. These conclusions hold for both the considered CQMs, which give rather close results. It is also interesting to study the impact of Melosh rotations directly on experimental related observables, such as σ ef f . We first consider the production of double J/Ψ via DPS at the LHC. Calculations are performed in the rapidity range |y| < 1.2 for ATLAS and CMS kinematics and 2 < y < 4.5 for the LHCb one. The calculation of σ ef f is performed via digluon distribution evaluated at Q 2 = 4m 2 c . In both these rapidity ranges, the involved parton momenta are quite small and we found that σ ef f is nearly constant. For this reason we just quote the averaged results in Tab. 4. The inclusion of Melosh operators determines an increase in σ ef f by almost 60%, whereas there is only a slight dependence on the chosen hadronic model. Then we consider double Higgs production via DPS in the same kinematic range. In this case, the digluon distribution is evaluated at Q 2 = m 2 H . The results for σ ef f , as a function of final state particle rapidities, are shown in Figs. 7 and 8. We note that σ ef f is almost constant in the central rapidity region, as already observed at Q 2 = 4m 2 c . However, for Q 2 = m 2 H in LHCb kinematics, the involved x i are substantially higher with respect to those addressed in the Q 2 = 4m 2 c case and σ ef f starts to show a non trivial x dependence. From this plots it is clear that the production of heavy objects in the forward rapidity region represents a way to access the kinematic region where longitudinal correlation are the strongest. For both the considered final scales, the inclusion of Melosh operators increase the value of σ ef f , as they act to reduce the size of dPDFs at small b ⊥ , as shown in Fig 5. The above results are similar, in quality, to that discussed in Ref. [56]. In order to further explore the role of Melosh operators in σ ef f , we consider the following ratio [12]: where in the denominator the effective cross section has been evaluated by means of gluon dPDFs calculated without Melosh rotations. Results of numerical calculations are presented in Fig. 9 for three fixed typical values of x 1 . Such a ratio shows a very weak dependence on x and the chosen model, and a weak dependence on the hard scale. Moreover its numerical value is found to be quite close to that obtained with valence quarks dPDFs evaluated at the hadronic scale of models described in Ref. [12]. It is interesting to note that Melosh's effects on σ ef f by far exceed the dependence induced by using different hadronic models.

Parametrization of Melosh effects
As shown in previous sections, Melosh effects on dPDF calculations are rather independent with respect to the adopted CQM, see as shown in Fig. 9, Melosh effects mildly depend on typical hard scales involved in the hard scatterings, either the J/Ψ or the Higgs mass in the present analysis. These features suggest that one may attempt to find a phenomenological parametrization of these effects which can be eventually implemented in dPDF analyses. For this purpose we define the ratio R between dPDFs calculated within CQM in a fully LF calculation and its approximation obtained neglecting relativistic effects expressed by Melosh operators: The above ratio is properly built in order to suppress dynamical effects encoded in the chosen hadronic wave function. Given a generic ansatz for dPDFs, indicated by F pheno , as those discussed in Refs. [28,30] and consistent with Eq. (13), under the assumption of a mild hard scale dependence of R, such an input could be improved to include relativistic effects by using the following prescription: Since we are interested in (x i − k ⊥ ) correlations induced only by Melosh operators, we have evaluated the ratio (61) at Q 2 = m 2 H within the only model which does not include any additional (x i −k ⊥ ) correlation generated by its wave function, i.e. the HO model [9,12]. We display in Fig. 10 the ratio R(x 1 , x 2 , k ⊥ ) for two representative values of k ⊥ as a function of x 1 and x 2 . We found that a suitable parametrization in (x 1 , x 2 ) space, which is able to describe the ratio R at fixed k ⊥ , is the following one: The coefficient w controls the overall normalization of the parametrization of R(x 1 , x 2 , k ⊥ ), t its small x slope and e and h its behaviour on the boundary. We note here that such a form goes beyond the standard factorized ansatz often used for dPDFs. By using the functional form in Eq. (63), we perform a series of fit of R(x 1 , x 2 , k ⊥ ) at fixed values of k ⊥ . This procedure gives us access to the k ⊥ dependence of the parameters which is displayed in Fig. 11. Then, we have fitted the k ⊥ dependence of the four parameters by using a fourth order polynomial of the type: where i = {w, e, t, h}. Since for k ⊥ = 0 Melosh operators reduce to unity, R(x 1 , x 2 , k ⊥ = 0) = 1 and therefore e(k ⊥ = 0) = h(k ⊥ = 0) = t(k ⊥ = 0) = 0. The latter condition is fulfilled by d e = d h = d t = 0 and d w = 1 which are held fixed at those values during the fit. The corresponding results are displayed as solid lines in   Table 5. The obtained parametrization reproduces with good accuracy (at the percent level) the ratio R(x 1 , x 2 , k ⊥ ) calculated within the HO model, as a function of the different variables and kinematics ranges, as shown by the R-sections in Fig. 12. We close this section discussing the scale dependence of the ratio R(x 1 , x 2 , k ⊥ ). For this purpose we present in the bottom sector of Fig. 12 the ratio R obtained with the full calculation of dPDFs (dotted lines), at different hard scales, compared to the R parametrizations (full lines) discussed above. In the bottom-left panel, the ratio R is presented at Q 2 = m 2 H as a function of k ⊥ . In the bottom-central panel, the full calculation of R has been obtained via dPDFs evaluated at the initial scale Q 2 = Q 2 0 in the valence region. Finally, in the bottom-right panel, the ratio R has been computed through dPDFs evaluated at the final scale Q 2 = 4m 2 c . As one can see, the parametrization of R fit well with its full calculation. Thanks to this feature, we conclude that, to high degree of accuracy, the R(x 1 , x 2 , k ⊥ ) parametrization is almost scale independent. This feature is consistent with the calculation of the ratio Eq. (59) by means of the full evaluation of dPDFs, see Fig. 6. Given this result, Melosh operators can be effectively accounted for by in any phenomenological dPDF parametrization a posteriori, via Eq. (62). Alternatively it can be directly supplemented to dPDFs input distributions at Q 0 prior to evolution, therefore completely removing the problem of the residual scale dependence of Eq. (63). We also notice that the functional form used to describe their effects is far beyond the standard factorization ansatz and therefore provide guidance in building non-factorized dPDF input.

An application of the parametrization of Melosh operators
In this section we describe a possible application of the parametrization of the Melosh effects obtained in the previous Section. We consider a standard, factorized, ansatz for the digluon distribution of the type: F pheno gg (x 1 , x 2 , k ⊥ , Q 2 ) = g(x 1 , Q 2 )g(x 2 , Q 2 )f (k ⊥ ) , where g(x, Q 2 ) stands for ordinary single gluon PDF eventually taken from parton distribution sets available in the literature. We assume that such a factorized ansatz works at the typical scale Q involved in a measured DPS processes. We indicate withf (k ⊥ ) a generic effective form factor. In the "Model I" scenario we use the gaussian form in Eq. (37) where the parameter a is chosen to give a mean value of σ N M ef f in the range 16 − 32 mb, in line with values extracted in experimental analyses. In the "Model II" scenario we consider the effective form factor which is proposed in Ref. [59] and related to the gluon form factor. Such a model is parameter free and it gives σ N M ef f = 32 mb. Combining Eq. (62) with the factorized ansatz in Eq. (65), the full dPDF reads Within this almost factorized scenario, the ratio (60) simplifies to: We show the ratio r σ calculated within the two described models, in Fig. (13) as a function of x 1 for two selected values of x 2 . The effect of including Melosh operators induces an enhancement on σ ef f close to the 30% for σ N M ef f ∼ 16 mb and 12% for σ N M ef f ∼ 32 mb. Furthermore, the ratio does depend on the numerical value of σ N M ef f , but is flat against kinematics variables and almost independent on the chosen functional form for the effective form factor, given a fixed value of σ N M ef f .

Conclusions
In the present analysis we have investigated to which extent information on the partonic proton structure, complementary to that obtained via other parton distribution, can be accessed via dPDFs. In particular we have focused our attention on the connection between the mean transverse partonic distance between two partons and σ ef f . We have discussed how this relation is modified when correlation of perturbative and nonperturbative origin are included in the calculation. In the former we have considered perturbative correlations induced by the so called splitting term in dPDF evolution. In the latter we have considered non perturbative correlation beyond the factorized ansatz for dPDFs. We proved that also in these two cases, the mean value of σ ef f provides new indications on the structure of the proton in the non perturbative regime of QCD, again indicating dPDFs as a valuable tool to investigate partonic longitudinal and transverse correlations. In the last part of this work we took advantage of CQM calculations of dPDFs within the Light-Front relativistic approach, to study model independent correlations between x 1 , x 2 and k ⊥ induced by the so called Melosh operators. We have investigated their effects on the digluon dPDF, perturbatively obtained at high momentum scales relevant for DPS studies at the LHC. We have shown that Melosh operators produce a non-negligible reduction of dPDFs and generate kinematical, model independent, (x i − k ⊥ ) correlations. Thanks to these results, we have also provided an economic parametrization of these effects which allows to to incorporate such relativistic corrections in dPDF ansatz.