Perturbative QCD correlations in multi-parton collisions

We examine the role played in double parton interactions (DPI) by the parton--parton correlations originating from perturbative QCD parton splittings. Also presented are the results of the numerical analysis of the integrated DPI cross sections at Tevatron and LHC energies. To obtain the numerical results the knowledge of the single-parton GPDs gained by the HERA experiments was used to construct the non-perturbative input for generalized double parton distributions. The perturbative two-parton correlations induced by three-parton interactions contribute significantly to resolution of the longstanding puzzle of an excess of multi-jet production events in the back-to-back kinematics observed at the Tevatron.

Multi-parton interactions can serve as a probe for nonperturbative correlations between partons in the nucleon wave function and are crucial for determining the structure of the underlying event at LHC energies. They constitute an important background for new physics searches at the LHC. A number of experimental studies were performed at the Tevatron [18][19][20]. New measurements are underway at the LHC [21,22].
Double hard parton scattering in hadron-hadron collisions (DPI) can contribute to production of four hadron jets with large transverse momenta p ⊥ Λ QCD , of two electroweak bosons, or "mixed" ensembles comprising three jets and γ, two jets and W , etc. In this paper we present the numerical results for a variety of final states. However, for the sake of definiteness, in what follows we refer to production of four final state jets in the collision of partons 1 and 2 from one incident hadron with partons 3 and 4 from the second hadron: 1 + 3 → J 1 + J 3 , 2 + 4 → J 2 + J 4 .
Double hard interaction is a process difficult to approach theoretically. Formally speaking, it calls for analysis of four-parton operators that emerge in the squared matrix element describing a two-parton state in a hadron. Relevant objects -quasi-parton operators -were introduced and classified, and evolution of their matrix elements studied by Bukhvostov et. al in the eighties in a series of papers [23][24][25].
An approach to MPI based on the operator product expansion and on the notion of transverse momentum dependent parton distributions (TMD) is being developed in [5][6][7].
MPI is a multi-scale problem. Not only because the separate parton-parton interactions may differ in hardness. More importantly, each single hard interaction possesses two very different hardness scales. The distinctive feature of DPI is that it produces two pairs of nearly back-to-back jets, so that Hence, in the collision of partons 1 and 3 the first (larger) scale is given by the invariant mass of the jet pair, Q 2 = 4J 2 1⊥ 4J 2 3⊥ , while the second scale is the magnitude of the total transverse momentum of the pair: δ 2 = δ 2 13 . It is important to stress that in the MPI physics there is no factorization in the usual sense of the word. The cross sections do not factorize into the product of the hard parton interaction cross sections and the multi-parton distributions depending on momentum fractions x i and the hard scale(s).
In [13] we have introduced necessary theoretical tools for approaching the problem by introducing the generalized double parton distributions ( 2 GPDs) in the momentum space. The double-parton GPD depends on one extra variable as compared to the double-parton distribution -the transverse momentum mismatch ∆ between the partons in the wave function and the wave function conjugated. In the mixed space of longitudinal momenta and transverse coordinates, an object equivalent to 2 GPD has been introduced by Treleani and Paver in the early eighties [26], and was long present in the literature since, see, in particular, [5,27,28].
In order to construct a viable model for the two-parton distribution stemming from the non-perturbative par-arXiv:1306.3763v3 [hep-ph] 29 May 2014 2 ton wave function of the proton and, in particular, for the ∆-dependence of the corresponding part of 2 GPD, we have exploited the existing experimental information obtained by the HERA experiments. This model disregards (rather arbitrarily, for lack of anything better) any long-range correlations between the partons in the proton wave function, which approximation seems reasonable in a limited range of parton momenta 10 −1 > x i > 10 −3 (see Section VI for details). This allows one to expresses the non-perturbative 2 GPD via the standard generalized parton distributions (GPDs) studied at HERA.
In [14] we have studied how do perturbative QCD (pQCD) phenomena affect two-parton correlations. One obvious effect is "scaling violation" due to parton multiplication processes in higher orders. Another important effect is short-range correlations induced by perturbative parton splitting, giving rise to three-parton collisions. Such contributions we will label as 1 ⊗ 2 ("3 → 4" of [13,14], "1v2" of [29]). This contribution to double hard interactions emerges as a result of collision of two partons from one hadron with the two offspring of the perturbative splitting of a single parton taken from another hadron. Perturbative splitting of one parton into two as a contributor to double-parton distributions was being discussed in the literature for a long time, see, e.g., [30,31]. However, for a relevant object -double-parton GPD -it was embedded into the parton evolution picture only recently.
The pQCD expressions for fully differential distributions were derived in the leading collinear approximation in [14]. In the present paper we examine how various contributions to the differential distribution, being integrated over transverse momentum imbalances, give rise to the expression for the integrated cross section in terms of the product of 2 GPDs of colliding hadrons.
The results of numerical studies of the integrated DPI cross sections based on the theory and the model for nonperturbative two-parton distributions in the proton developed earlier are reported. We demonstrate that the 1 ⊗ 2 processes contribute significantly to the DPI cross section and, in particular, are capable of explaining the longstanding Tevatron puzzle [18][19][20]. Namely, that the DPI cross section extracted by the the CDF and D0 at x i ∼ 0.01 turned out to be a factor of two larger than the expectation based on the approximation of independent partons in the proton with the transverse distance spread extracted from the HERA data [12].
This observation does not exclude the presence of genuinely non-perturbative correlations between partons in the proton. However, we find it interesting that the pQCD induced correlations alone could explain the scale of the enhancement found by the Tevatron experiments.
We also discuss the pattern of the x-and Q 2dependence of the effective interaction area (aka "effective cross section") induced by pQCD correlations at Tevatron and LHC energies.
The paper is organized as follows.
In Section II we reflect upon intrinsic complexity of theoretical MPI analysis. In Section III we remind the reader the basics of our approach to generalized twoparton distributions, 2 GPD, and discuss the approximations adopted for constructing the model for 2 GPD in terms of the standard GPDs.
In Section IV we discuss importance of perturbative parton correlations and give a semi-quantitative estimate of their magnitude. Section V is devoted to three-parton interactions (1 ⊗ 2 subprocesses). Here we reexamine the role of the so-called "short split" term in 1 ⊗ 2 subprocesses which is concentrated in the kinematical region of nearly equal transverse momentum imbalances, see [14]. We show here that the "short split" contribution to the integrated cross section is actually contained in σ 3 that one obtains integrating the simple DDT-like formula derived in the complimentary kinematical region δ 2 13 δ 2 24 (δ 2 13 δ 2 24 ). Thus we correct (and simplify) the expression for the integrated DPI cross section (original Eq. (28) of [14] was plagued by double counting).
In Section VI we discuss the independent parton model for the non-perturbative part of two-parton correlations. The numerical results for Tevatron and LHC energies are presented in Section VII. We conclude and discuss the results and perspectives of DPI studies in Section VIII.

II. HIDDEN REEFS OF DPI PHYSICS
An important question is, whether MPI admit an intuitive probabilistic parton interpretation as do the classical single hard interaction processes.
When one considers inclusive one-parton distribution in a hadron (pdf), the quantum state of the parton in the light-cone hadron wave function (w.f.) coincides with that in the wave function conjugated (w.f.c.). This lays down the foundation for the probabilistic QCD improved parton picture.
In the case of the two-parton correlation the situation is different. Here only the overall quantum characteristics of the parton pair as a whole (its total energymomentum, its spin and color states) should be identical in the w.f. and w.f.c. As a result, the momentum, spin and color state of a single parton in the pair do not necessarily match in the wave function and in the wave function conjugated, thus endangering the very probabilistic interpretation of the process under consideration.
A general approach to double hard interactions has been developed in [13]. It turned out that the transverse momentum of the parton in the w.f. and that of its counterpart in the w.f.c. are indeed necessarily different, with their difference ∆ being conjugate to the relative transverse distance between the two partons in the hadron. This has led to introduction of the notion of the generalized double parton distribution, 2 GPD, which depends on 3 a new momentum parameter ∆. It is important to stress that the cross section of the double hard process does not factorize into the product of the hard cross section and parton distributions. Instead, it contains a convolution of the product of two 2 GPDs over d 2 ∆.
Non-diagonality in the longitudinal momentum fractions is there too. Representing incident partons as plain waves with definite momenta, one ignores the fact that the two partons originate from one and the same finite size hadron. Therefore, when one picks two partons with momenta x 1 , x 2 from the hadron wave function, one has to integrate over x 1 − x 2 at the level of the amplitude, which integration ensures that the longitudinal separation between the partons does not exceed the size of the parent hadron. When one considers independent hard interactions of two pairs of partons, (x 1 , x 3 ) and (x 2 , x 4 ), taken from colliding hadrons, integrals over r = x 1 − x 2 and r = x 3 − x 4 do not manifest themselves, since all four longitudinal momentum fractions x i are uniquely determined by the kinematics of the two produced hard systems.
Not so for 1 ⊗ 2 processes. In this case a flow of large momenta between the two hard vertices is possible, and the integration over r is instrumental in getting rid of a fake singularity of the scattering matrix element (a simple explanation of the origin of this unphysical singularity can be found in [32], see also [14]).
As a result, a small offset between the parton longitudinal momenta in the w.f. and w.f.c. emerges, |x 3 − x 3 | ∝ δ 2 /Q 2 . However, this mismatch turns out to be negligible in the dominant kinematical region where the squared total transverse momentum δ 2 of the jet pair is much smaller than the overall hardness Q 2 of the process: δ 2 13 ≡ ( p ⊥1 + p ⊥3 ) 2 Q 2 1 4p 2 ⊥1 . As noticed by Gaunt [33], the non-diagonality in the longitudinal momentum space is likely to get induced in 1 ⊗ 2 processes by higher order QCD effects. This happens when incoming partons, in the w.f. and w.f.c., exchange (real or virtual) gluons with transverse momenta in the interval Λ 2 QCD k 2 ⊥ δ 2 . In [33,34] arguments were raised in favor of smallness of the crosstalk effects in the DPI cross section. Otherwise, this would have led to unwelcome complication of the problem: one would be forced to deal with an unknown function of four longitudinal momentum fractions (three independent variables) in place of the two of 2 GPD(x 1 , x 2 ). Hence for the time being we disregard this complication.
We do not dwell into potential non-diagonality in color and in spin variables either. One may argue that such non-diagonal configurations are likely to be suppressed as they can be related with form factors for proton transition between two states with different quantum numbers of the proton constituents. For example, swapping inside the proton the colors of two partons that sit at distances of the order of the nucleon radius, corresponds to excitation that can be visualized as adding an extra piece of color string whose energy, O (1 GeV), would excite and destroy the proton.

III. GENERALIZED TWO-PARTON DISTRIBUTION
In [13,14] we have developed a formalism to address the problem of multi-parton interactions. QFT description of double hard parton collisions calls for introduction of a new objectthe generalized two-parton distribution, Here the index h refers to the hadron, x 1 and x 2 are the light-cone fractions of the parton momenta, and Q 2 1 , Q 2 2 the corresponding hard scales. As has been mention above, the two-dimensional vector ∆ is Fourier conjugate to the relative distance between the partons 1 and 2 in the impact parameter plane. The distribution obviously depends on the parton species; we suppress the corresponding indices for brevity. The double hard interaction cross section (and, in particular, that of production of two dijets) can be expressed through the generalized two-parton distributions 2 GPD.
2 GPDs enter the expressions for the differential distributions in the jet transverse momentum imbalances δ ik in the kinematical region (1), as well as for the "total" DPI cross section (integrated over δ ik ). In the latter case the hardness parameters of the 2 GPDs are given by the jet transverse momenta Q 2 i , while in the differential distributions -by the imbalances δ 2 ik themselves. The corresponding formulae derived in the leading collinear approximation of pQCD can be found in Ref. [14].
It is important to bear in mind that the DPI cross section does not factorize into the product of the hard parton interaction cross sections and the two two-parton distributions depending on momentum fractions x i and the hard scales, Q 2 1 , Q 2 2 . Instead, .

(2b)
The effective interaction area σ eff is given by the convolution of the 2 GPDs of incident hadrons over the transverse momentum parameter ∆ normalized by the product of single-parton inclusive pdfs.
For expression (2) for the DPI cross section to make sense, the integral over ∆ has to be convergent. This is well the case when the two partons are taken from the non-perturbative (NP) proton wave function. Indeed, a typical inter-parton distance in the proton is large, of the order of the hadron size R. Accordingly, one expects the corresponding correlator in the momentum space to be concentrated at a finite NP scale ∆ 2 ∼ R −2 and to fall fast at large ∆ 2 (exponentially or as a sufficiently high power of ∆).
However, there is another source of two-parton correlations. This is purely perturbative (PT) mechanism when the two partons emerge from perturbative splitting of one parton taken from the hadron wave function. In this scenario the production of the parton pair is concentrated at much smaller distances. As a result, the corresponding contribution to 2 GPD turns out to be practically independent on ∆ 2 in a broad range, up to the hard scale(s) characterizing the hard process under consideration (∆ 2 only affects the lower limit of transverse momentum integrals in the parton cascades, causing but a mild logarithmic dependence).
Given essentially different dependence on ∆, one has to treat the two contributions separately by casting the 2 GPD as a sum of two terms: Here subscripts [2] D and [1] D mark the first and the second mechanisms, correspondingly: two partons from the wave function versus one parton that perturbatively splits into two.

IV. PERTURBATIVE TWO-PARTON CORRELATIONS
In this Section we discuss the role of the PT parton correlations and show that, given a sufficiently large scale of hard interactions, they turn out to be as important as NP ones.
A. Estimate of the PT correlation.
Let us chose a scale Q 2 sep that separates NP and PT physics to be sufficiently low, so that parton cascades due to evolution between Q 2 sep and Q 2 i are well developed. To get a feeling of relative importance of the PT correlation, as well as to understand its dependence on x and the ratio of scales, Q 2 vs. Q 2 sep , the following lowest order PT estimate can be used.
Imagine that at the scale Q 2 sep the nucleon consisted of n q quarks and n g gluons ("valence partons") with relatively large longitudinal momenta, so that triggered partons with x 1 , x 2 1 resulted necessarily from PT evolution. In the first logarithmic order, α s log(Q 2 /Q 2 sep ) ≡ ξ, the inclusive spectrum can be represented as where we suppressed x-dependence as irrelevant. If both gluons originate from the same "valence" parton, then while independent sources give [2] D: Recall that the ∆-dependence is different in (4a) and (4b). However, at ∆ = 0 the second terms cancel in the sum and we get for the correlator .
The correlation is driven by the gluon cascade --the first term in (4a) -and is not small (being of the order of unity). It gets diluted when the number of independent "valence sources" at the scale Q 2 sep increases. This happens, obviously, when x i are taken smaller. On the other hand, for large x i ∼ 0.1 and increasing, the effective number of more energetic partons in the nucleon is about 2 and decreasing, so that the relative importance of the 1 ⊗ 2 processes grows.
We conclude that the relative size of PT correlations is of order one, provided ξ = O (1).
Moreover, the PT parton correlations cannot be disregarded without running a risk of violating general principles. This can be illustrated by looking at the momentum sum rule for double parton distributions.

B. Momentum sum rule
An obvious momentum sum rule should be satisfied. Namely, that the integral over dx 2 with the weight x 2 (summing over all parton species) should produce in the end (1 − x 1 ) times the inclusive one-parton distribution D(x 1 ), that is, the total longitudinal momentum carried by all the partons but the triggered one: (here we have explicitly restored the parton species indices i 1 , i 2 ). This sum rule, together with other ones concerning valence quantum numbers, has been discussed by Gaunt and Stirling in [11] as means for restricting the form of double parton distributions. Setting the ∆ argument of 2 GPD to ∆ = 0 corresponds to taking integral over the relative transverse distance between partons in the proton. Derivation of (6) is carried out in Appendix A. It explicitly demonstrates that the PT parton splitting enters on equal grounds with the contribution due to two partons taken both from the initial NP hadron wave function.

V.
1 ⊗ 2 DPI PROCESS Actually, the NP and PT contributions do not enter the physical DPI cross section in arithmetic sum (3), driving one even farther from the familiar factorization picture based on universal (process independent) parton distributions. As explained in [14], a double hard interaction of two pairs of partons that both originate from PT splitting of a single parton from each of the colliding hadrons, does not produce back-to-back dijets. In fact, such an eventuality corresponds to a one-loop correction to the usual 2 → 4 jet production process and should not be looked upon as multi-parton interaction. The term [1] D h1 × [1] D h2 has to be excluded from the product D h1 × D h2 , the conclusion we share with Gaunt and Stirling [29].
So, we are left with two sources of genuine twoparton interactions: four-parton collisions described by the product of (PT-evolved) 2 GPDs of NP origin (2 ⊗ 2), and three-parton collisions due to an interplay between the NP two-parton correlation in one hadron and the two partons emerging from a PT parton splitting in another hadron (1 ⊗ 2 ), described by the combination Given that [2] D falls fast at large ∆, a mild logarithmic ∆-dependence of [1] D can be neglected in the product in (7b).
A. On separation of NP and PT parts of 2GPD. Parameter Q 2 0 .
Separation of PT and NP contributions is a delicate issue. By definition of the perturbative correlation function, [1] D vanishes when Q 2 1 , Q 2 2 are taken equal to the separation scale Q 2 sep that one chooses to set the lower limit for applicability of pQCD calculations. Strictly speaking, Q 2 sep can be chosen arbitrarily: both the NP input function [2] D and the PT-calculable correlation [1] D contain Q 2 sep -dependence, but their sum does not depend on this formal parameter. At the same time, the character of the ∆-dependence of [2] D depends, obviously, on the choice of the Q 2 sep scale. Indeed, by increasing the value of Q 2 sep one will shuffle a part of the perturbative splitting contribution from [1] D into [2] D. As a result, the "NP correlator" [2] D, contaminated by a short range PT correlation, would acquire a "tail" at large ∆ 2 which would spoil convergence of the ∆ integration in (2).
Thus, in order to preserve the logic of the NP-PT separation, one is led to introduce a specific resolution scale, 0 is no longer an arbitrary "factorization scale" but a phenomenological parameter whose value (which one expects to be of order of 1GeV) should be established from the data.

B. Composition of the 1 ⊗ 2 DPI cross section
In order to derive the DPI cross section, one has to start with examination of the double differential transverse momentum distribution and then integrate it over jet imbalances δ ik . Why this step is necessary?
The parton distribution D(x, Q 2 ) -the core object of the QCD-modified parton model -arises upon logarithmic integration over the transverse momentum up to the hard scale, k 2 ⊥ < Q 2 . Analogously, the double parton distribution D(x 1 , x 2 , Q 2 1 , Q 2 2 ; ∆) embeds independent integrations over parton transverse momenta k 2 1⊥ , k 2 2⊥ up to Q 2 1 and Q 2 2 , respectively. However, the 1 ⊗ 2 DPI cross section contains a specific contribution ("short split", see below) in which the transverse momenta of the partons 1 and 2 are strongly correlated (nearly opposite). This pattern does not fit into the structure of the pQCD evolution equation for 2 GPD where k 1⊥ and k 2⊥ change independently. Given this subtlety, a legitimate question arises whether the expression for the integrated 1 ⊗ 2 cross section (7b) based on the notion of the two-parton distribution [1] D takes the short split into account. Below we demonstrate that in fact it does.
The differential distribution over jet imbalances was 6 derived in [14] in the leading collinear approximation of pQCD. It resembles the "DDT formula" for the Drell-Yan spectrum [43] and contains two derivatives of the product of 2 GPDs (7) that depend on the corresponding δ ik as hardness scales, and the proper Sudakov form factors depending on (the ratio of) the Q 2 i and δ 2 ik . In particular, in the region of strongly ordered imbal-ances, ; δ 2 13 δ 2 24 , δ 2 13 δ 2 24 , (8) the differential 1 ⊗ 2 cross section reads The differential distribution for the 2 ⊗ 2 DPI mechanism has a similar structure, see Eq. (25) of [14].
In addition to (9), there is another source of double collinear enhancement in the differential 1 ⊗ 2 cross section. It is due to the kinematical region where the two imbalances nearly compensate one another, and the dominant integration region is complementary to that of (8): This enhancement characterizes the set of 1 ⊗ 2 graphs in which there is no accompanying radiation with transverse momenta exceeding |δ |. In this situation, the parton that compensates the overall imbalance, k ⊥ = − δ is radiated off the incoming, quasireal, parton legs as shown in Fig. 1. At the same time, the virtual partons after the core splitting "0"→ "1"+"2" enter their respective hard collisions without radiating any offspring on the way. The 1 → 2 splitting neighbors the hard vertices, therefore the name "short split" (aka "endpoint contribution", [14]).
A closed expression for the differential distribution in jet imbalances due to short split, derived in the leading collinear approximation, is given by Eq. (27) of [14]: The short split becomes less important when the scales of the two hard collisions separate. Indeed, the logarithmic integration over δ 2 is kinematically restricted from above, δ 2 < δ 2 max min{Q 2 1 , Q 2 2 }. As a result, when transverse momenta of jets in one pair much exceed those of the second pair, e.g., Q 2 1 Q 2 2 (see (1)), the contribution of the short split becomes suppressed as Here S 1 and S 3 are the double logarithmic Sudakov form factors of the partons "1" and "3" that enter the hard interaction with the larger hardness scale.
Short split induces a strong correlation between jet imbalances which is worth trying to look for experimentally.
The relative weight of the short split depends on the process under consideration. For most DPI processes in the kinematical region we have studied, it typically provides 10-15% of the R value. However, it happens to be more important when the nature of the process favors parton splitting. This is markedly the case of the double Drell-Yan pair production where the short split contribution reaches 30-35%. On the contrary, it turns out to be practically negligible for the same-sign double W -meson production (see discussion in the Conclusions section below).

C. Short split in the integrated cross section
In the preceding publication we have treated the contribution of the short split to the total DPI cross section as an addition, 1/σ short , to the 1 ⊗ 2 effective interaction area 1/σ 3 , see Eqs. (28), (32c) of [14]. However, this was not a right thing to do. As it turns out, the short split contribution to the integrated cross section is already contained in (17b).
To see how this happens, one has to examine the structure of the DDT formula for the 1 ⊗ 2 differential cross section (9) more closely. The distribution function [1] D in (9) in the leading collinear approximation is given by Eq. (18) of [14]: Here a, b mark the registered partons, and a , b , c are the indices of the partons involved in the splitting c → a (z) + b (1 − z) with the DGLAP probability P (a ,b ) c (z). The distribution D c h (y, k 2 ) describes the standard probability of finding a parton c inside the incident hadron h at the transverse momentum scale k 2 , and the functions D i i (x, q 2 ; k 2 ) in the second line stand for the distribution of parton i, probed at scale q 2 , in the initial parton i at lower virtuality scale k 2 .
Let us examine the structure of the terms one gets applying two derivatives to (13) substituted into (9): Taking derivative over q 2 i of the product D(x, q 2 i ) × S(Q 2 , q 2 i ) corresponds to picking up from the parton chain an accompanying parton (with the largest transverse momentum) which compensates the imbalance, = −q i .
Apart from logarithmic dependences of D and S, the derivative in (9) may act upon the upper limit of the k 2 integration in (13). α s q 2 1 α s q 2 2 min (q 2 1 ,q 2 2 ) dk 2 k 2 α s (k 2 )F (k 2 ; q 2 1 , q 2 2 ). (15a) Differentiating the upper limit of the k 2 -integral de-scribes another legitimate situation when one of the partons "1" and "2" does not radiate before entering the hard interaction. In this case the smaller of the jet imbalances is determined by the splitting momentum k: Finally, applying both derivatives to the integration limit (13) gives rise to Contrary to the first two contributions (15a) and (15b), the last term (15c) obviously violates the condition of applicability (8) of the DDT formula. Instead, in the region δ 2 13 ∼ δ 2 24 it is the short split that contributes in the leading collinear approximation, so that Eq. (12) has to be used to describe the differential spectrum in place of (15c).
However, as far as the total cross section is concerned, the integrals over imbalances of the short split and of the fake singular term (15c) turn out, as by miracle, to be the same. Indeed, integrating the short split (12) over δ 2 up to δ 2 , we get . On the other side, taking the integrand of [1] D (13) in the point q 2 1 = q 2 2 = k 2 , and using D a a x 1 zy , q 2 1 ; k 2 we evaluate the momentum integrals, to arrive at the very same expression (16).
It is worth noticing that this correspondence does not depend on the precise form of the upper integration limit in (13). The result does not change, within the leading logarithmic accuracy, if one replaces a sharp ϑ-function cut by a smooth damping factor that cuts the logarithmic k 2 integration at k 2 ∼ min{q 2 1 , q 2 2 }.
Thus, for the integrated DPI cross section we obtain two contributions to the effective interaction area: Let us stress in conclusion that a compact and intuitively clear expression containing the product of the 2 GPDs [2] D and [1] D in (17b) applies only to the integrated 1 ⊗ 2 cross section. When addressing the differential distributions, one has to employ the "DDT-like formula" (9) in the region of strongly ordered transverse momenta (8), and a quite different expression (12) in the kinematical region of nearly opposite jet pair imbalances (11).

VI. MODELING [2] D
To proceed with quantitative estimates, one needs a model for the non-perturbative two-parton distributions in a proton.
A priori, we know next to nothing about them. The first natural step to take is an approximation of independent partons. It allows one to relate 2 GPD with known objects, namely [13] [2] D(x 1 , x 2 , Q 2 1 , Q 2 2 ; ∆) G(x 1 , Q 2 1 ; ∆ 2 )G(x 2 , Q 2 2 ; ∆ 2 ). (18a) Here G is the non-forward parton correlator (known as generalized parton distribution, GPD) that determines, e.g., hard vector meson production at HERA and which enter in our case in the diagonal kinematics in x (x 1 = x 1 , see Fig. 2).

FIG. 2: GPD in the vector meson electroproduction amplitude
The modeling by (18a) is not perfect. First of all, it does not respect an obvious restriction D(x 1 + x 2 > 1) = 0. So, x i have to be taken not too large (say, x i < 0.5).
Actually, x i must be taken even smaller. The GPD in Fig. 2 is an elastic amplitude, while the corresponding block in the DPI represents the inclusive cross section (the cut-through amplitude). For the analogy to hold, the interaction amplitude has to be close to imaginary. This condition calls for x i < 0.1.
On the other hand, x i should not be too small to stay away from the region of the Regge-Gribov phenomena where there are serious reasons for parton correlations to be present at the non-perturbative level (see discussion in [35]).
Thus, we fix the domain of applicability of the model (18a) for 2 GPD as 10 −1 > x i > 10 −3 .
The GPD, in its turn, can be modeled as with D the usual one-parton distribution determining DIS structure functions and F the so-called twogluon form factor of the hadron. The latter is a nonperturbative object; it falls fast with the "momentum transfer" ∆ 2 . This form factor can be parametrized differently. For example, by a dipole formula where an effective parameter m 2 g extracted from the FNAL and HERA J/ψ exclusive photoproduction data lies in the ballpark of m 2 g (x ∼ 0.03, Q 2 ∼ 3 GeV 2 ) 1.1 GeV 2 and decreases with further decrease of x [36].
Turning to the 1⊗2 term, we neglect a mild logarithmic ∆-dependence of [1] D in (17b) and use the model (18) for where we substituted the value of the integral (cf. (19)) We will parametrize the result in terms of the ratio For the effective interaction area, we then have (the phenomenological value m 2 g = 1.1 GeV 2 was used).
Within the framework of the NP two-parton correlations model (18a), there is but one free parameter Q 2 0 . The DPI theory applies to various processes and holds in a range of energies and different kinematical regions. Therefore, having fixed the Q 2 0 value, say, from the Tevatron data, one can consider all other applications (in particular, to LHC processes) as parameter-free theoretical predictions.

A. Calculation framework
In numerical calculations we used the GRV92 parametrisation of gluon and quark parton distributions in the proton [37]. We have checked that using more advanced GRV98 and CTEQ6L parametrisations does not change the numerical results. The explicit GRV92 parametrisation is speed efficient and allows one to start the PT evolution from rather small virtuality scales. The combination (Q 2 0 + ∆ 2 ) was used as the lower cutoff for logarithmic transverse momentum integrals in the parton evolution, which induced a mild (logarithmic) ∆dependence on top of the relevant power of the two-gluon form factor F 2g (∆ 2 ).

CDF experiment
In Fig. 3 we show the profile of the 1 ⊗ 2 to 2 ⊗ 2 ratio R for the γ + 3jets process in the kinematical domain of the CDF experiment [18]. The calculation was performed for the dominant "Compton scattering" channel of the photon production: g(x 2 ) + u(ū)(x 4 ) → γ + u(ū). The longitudinal momentum fractions of two gluons producing second pair of jets are x 1 and x 3 . The typical transverse momenta were taken to be p ⊥1,3 5 GeV for the jet pair, and p ⊥2,4 20 GeV for the photon-jet system. In Fig. 3 R is displayed as a function of rapidities of the photon-jet, η 2 = 1 2 ln(x 2 /x 4 ), and the 2-jet system,  (21) in the CDF kinematics for the process pp → γ + 3 jets + X.
In our previous report [35] we have included a plot of the x-dependence of R for central production at Tevatron and LHC energies. That plot turned out to be confusing: a rather sharp x-dependence it has demonstrated seemed to contradict the CDF findings of approximate constancy of σ eff . In fact, this variance with x was in the major part resulting from the kinematical link between x and Q 2 for a given collision energy (Q 2 = x 2 s/4).
The results of numerical calculation for a fixed hardness Q 2 shown in Fig. 3 for the CDF kinematics exhibit a very mild x-dependence of the R factor and thus of the σ eff .

D0 experiment
The ratio R is practically constant in the kinematical domain of the D0 experiment on photon+3 jets production [19,20] and is very similar to that of the CDF experiment shown above in Fig. 3. So, for the D0 kinematics we instead display in Figs. 4, 5 the enhancement factor 1+R in dependence on p ⊥ of the secondary jet pair for photon transverse momenta 10, 20, 30, 50, 70, and 90 GeV(from bottom to top).

B. Models of parton spatial density
In this section we study the limits that can be obtained on the parameters of three phenomenological models of parton spatial density using the measured effective cross section (16). In the discussion below we follow a simple classical approach. For a given parton spatial density inside the proton or antiproton ρ(r), one can define a (time-integrated) overlap O(β) between the parton distributions of the colliding nucleons as a function of the impact parameter β [10]. The larger the overlap (i.e. smaller β), the more probable it is to have at least one parton interaction in the colliding nucleons. The single hard scattering cross sections (for example, γ + jets or dijet production) should be proportional to O(β) and the cross section for the double parton scattering is proportional to the squared overlap, both integrated over all impact parameters β [28,36]: First, we consider the "solid sphere" model with a constant density inside the proton radius r p . In this model, the total hard scattering cross section can be written ton density ρ(r) = constant for r ≤ r for r > r p and found to be f = 2.1 the enhancement factor can be seen bet Eq. (1) as σ DP = f σ A σ B /σ hard . The ha parton interaction is the more it is bia central hadron-hadron collision with a s rameter, where we have a larger overlap sities and, consequently, higher proba ond parton interaction [5]. Using the for the solid sphere model we extract dius r p = 0.53 ± 0.06 fm and proton rm 0.41 ± 0.05 fm. The latter is obtained r 2 as R 2 rms ≡ ∞ 0 r 2 4πr 2 ρ(r)dr = 4π The results are summarized in the line of Table VI. The Gaussian model with and exponential model with ρ(r) ∝ e −r/ tested. The relationships between the (r p , a or b) and rms-radius for all the mo Table VI. The relationships between th section σ eff and parameters of the Gau nential models are taken from [38], negl that represent correlations in the transv scale parameters and rms-radii for both given in Table VI. In spite of difference the proton rms-radii are in good agree other, with average values varied as 0.41 about 12% uncertainty. On the other h tained rms-radius from other sources (fo and using the measured σ eff , the size o correlations [38] can be estimated.

IX. SUMMARY
We have analyzed a sample of γ + 3 lected by the D0 experiment with an nosity of about 1 fb −1 and determined events with hard double parton scatter a single pp collision at √ s = 1.96 TeV. are measured in three intervals of the in p T ) jet transverse momentum p jet2  The corresponding prediction for σ eff is shown in Fig. 6 in comparison with the D0 findings. Both the absolute value and (a hint at) the p ⊥ -trend look satisfactory.

C. LHC energies
In Fig. 7 we show the 1⊗2 to 2⊗2 ratio for production of two pairs of back-to-back jets with transverse momenta 50 GeV produced in collision of gluons at the LHC energy √ s = 7 TeV.  Dependence on the hardness parameters of the DPI process of double gluon-gluon collisions is illustrated in Fig. 8. For the sake of illustration, we have chosen the value of the PT cutoff parameter Q 2 0 = 0.5 GeV 2 , and calculated the enhancement factor 1 + R for five values of the transverse momenta of the jets in one pair, p ⊥1 = 20, 40, 60, 80, 100 GeV. Fig. 8 demonstrates its dependence on the transverse momenta of jets in the second pair, p ⊥ ≡ p ⊥2 .
We observe that within the chosen range, R increases by about 15-25% with increase of the hardness of one of the jet pairs. This corresponds to approximately 10% drop in σ eff .

12
Finally, in Fig. 9 we show the rapidity profile of the R ratio for the process of production of the vector boson, ud → W + , accompanied by an additional pair of (nearly back-to-back) jets with transverse momenta p ⊥ = 30 GeV produces in a gluon-gluon collision. It is interesting to notice that the effect of perturbatively induced parton-parton correlations is maximal for equal rapidities of the W and the jet pair, and diminishes when they separate. This feature is more pronounced when the cutoff parameter Q 2 0 is taken larger, so that the PT correlation becomes smaller and, at the same time, exhibits a stronger rapidity dependence.
The recent ATLAS study [39] reported for this process the value σ eff = 15 ± 3 +5 −3 mb which is consistent with the expected enhancement due to contribution of the 1 ⊗ 2 DPI channel, see (21).

D. Q0-dependence
The dependence of the enhancement factor 1 + R on the Q 0 parameter is shown for the typical kinematics of the CDF photon+3 jets experiment in Fig. 11, and for central production of two pairs of p ⊥ 50 GeV jets at the LHC (η 1 = η 2 ) in Fig. 12.

VIII. CONCLUSIONS AND DISCUSSION
In the previous paper [14] we have analyzed the perturbative correlation that arises due to 1 ⊗ 2 splitting of the parton in one of the colliding hadrons and derived the corresponding expressions in the leading collinear approximation of the pQCD. Here we presented results of numerical evaluation of this contribution to the DPI cross section measured at the Tevatron and found the theoretical results to be consistent with the data for the value of the model parameter Q 2 0 0.5 GeV 2 . With Q 2 0 fixed, theoretical expectations for certain exemplary DPI processes at LHC energies become parameter-free predictions.
Theoretical derivation of the effective interaction area σ eff ("effective cross section") relied on certain assumptions and approximations. Our approach to perturbative QCD effects in DPI developed in [14] was essentially probabilistic. In particular, we did not discuss the issue of possible interference between 1 ⊗ 2 and 2 ⊗ 2 twoparton amplitudes. One can argue that such eventuality should be strongly suppressed. Indeed, spatial properties of accompanying radiation produced by so different configurations make them unlikely to interfere, since in the 2⊗2 mechanism a typical transverse distance between two partons from the hadron w.f. is of order of the hadron size, while in the 1 ⊗ 2 case it is much smaller and is determined by a hard scale. Moreover, we disregarded potential contributions from non-diagonal interference diagrams that are due to crosstalk between partons in the amplitude and the amplitude conjugated. Such contributions are absent in 2 ⊗ 2 collisions but emerge specifically in the 1⊗2 DPI process in which the two partons from the hadron wave function are relatively close to one another in the impact parameter plane [33].
Our prediction for σ eff was obtained as the ratio of the DPI and SPI cross sections derived in the leading collinear approximation of pQCD (LLA). This means that evolution of perturbative parton cascades was treated at the one-loop level, and the matrix ele-ments of hard parton interactions treated in the Born approximation. In the DPI problem, the subleading nonlogarithmic corrections to the LLA are bound to be sizable. Indeed, when deriving the total DPI cross section within the LLA accuracy, one integrates the differential distribution over jet imbalances, δ 2 ik Q 2 i , up to the scale given by the transverse momenta of the jets, Q 2 i . In reality, due to experimental cuts that are imposed in order to extract jets in the back-to-back kinematics the true hard scale of the DPI cross section is lower. Being formally a subleading O (α s ) correction, it will affect both the 1 ⊗ 2 and 2 ⊗ 2 cross sections. Which way the subleading pQCD effects will change the ratio is so far unknown. To establish the true hard scales of the parton distributions entering the DPI cross section formula, one has to carry out the NLLA analysis which would include taking into consideration concrete details of the jet finding algorithms employed in the experimental setup.
Finally, our prediction for the DPI cross sections was based on a model assumption of the absence of NP twoparton correlations in the proton. This assumption is arbitrary. One routinely makes it for the lack of any firsthand information about such correlations. In [35] we have pointed out a source of genuine non-perturbative two-parton correlations that should come onto the stage for very small x values, x 10 −3 , and estimated its magnitude via inelastic diffraction in the framework of the Regge-Gribov picture of high energy hadron interactions. Also it was argued in [40] that strong quarkantiquark correlations may arise from dynamical chiral symmetry breaking.
In order to be able to reliably extract the DPI physics, one has to learn how to theoretically predict 1 ⊗ 1 parton collision processes with production of two hard systems (four jets in particular). This is the dominant channel, and it is only in the back-to-back kinematics that the 2⊗2 and 1⊗2 DPI processes become competitive with it. Among first subleading pQCD corrections to the 1⊗1 amplitude, there is a loop graph that looks like two-by-two parton collision. This resemblance is deceiving though. Unlike the 2 ⊗ 2 and 1 ⊗ 2 contributions, this specific correction does not depend on the spatial distribution of partons in the proton (information encoded in σ eff ), is not power enhanced in the region of small transverse momenta of hard systems, and therefore does not belong to the DPI mechanism [14,29]. Treating the 1 ⊗ 1 amplitude at the one-loop level corresponds to the twoloop accuracy for the cross section. Until this accuracy is achieved, the values of σ eff extracted by experiments should be considered as tentative.
Our first conclusion is that in the kinematical region explored by the Tevatron experiments, the x-dependence of σ eff turns out to be rather mild. This by no means implies, however, that σ eff can be looked upon as any sort of a universal number. On the contrary, we see that the presence of the perturbative correlation due to the 1 ⊗ 2 DPI mechanism results in dependence of σ eff not 14 only on the parton momentum fractions x i and on the hardness parameters, but also on the type of the DPI process.
For example, in the case of golden DPI channel of production of two same sign W bosons [41] the discussed mechanism leads to expectation of significantly larger σ eff than for, say, W plus two jets process. Indeed, the comparison of the values of R for central production of two gluon jet pairs, W jj and W + W + (with jet transverse momenta p ⊥ M W /2), gives ( √ s = 7 TeV, η 1 = η 2 = 0) As a result of the varying magnitude of the perturbative correlation, the effective interaction areas σ eff come out to be significantly different for the three processes: The smaller value for each effective interaction area corresponds to more developed perturbative parton cascades than the larger one (Q 2 0 = 0.5 GeV 2 versus Q 2 0 = 1.0 GeV 2 ).
Contrary to the W + W + channel, the double Drell-Yan process favors the 1 ⊗ 2 mechanism, g → uū. As a result, the effective interaction area here turns out to be even smaller.
For central production of two Z bosons at √ s = 7 TeV we get R(ZZ) = 1.03 (0.73), ZZ : σ eff = 15.9 ÷ 18.5 mb. (26) An important feature of the 1 ⊗ 2 mechanism is its dependence on the hardness of the process. With increase of Q 2 i , the 1 ⊗ 2 to 2 ⊗ 2 ratio R should increase rather fast thus pushing σ eff to smaller values. At the same time, with decrease of the p ⊥ of the jets this contribution decreases. As we have seen above, such a trend is consistent with the D0 data for x ∼ 10 −2 .
By pushing the hardness scales down to p ⊥ ∼ 3÷4 GeV, one enters the domain of the physics of minijets. Here one should have σ eff ∼ 25 mb for x i ∼ 10 −2 , the much larger value than the Q-independent σ eff which was assumed in the Monte Carlo models like PYTHIA for a long time.
It would be interesting to implement in the MC models a more realistic account of MPI in which σ eff would decrease with increase of p ⊥ .
Due to the momentum sum rule for the parton wave function of the proton, summing over parton species B produces the total energy carried by all initial partons but the parent of the second registered parton i 1 : Bȳ B = 1 − y 1 .
The parent parton energy y 1 is not observable. What is fixed by the measurement is x 1 , while y 1 is being inte-15 grated over: where we have introduced an average conditional parent parton energy y 1 which depends on x 1 and Q 2 . It is important to notice that this quantity depends also on the unphysical scale Q 2 sep which separates the domains of PT and NP description (unlike the physical inclusive parton distribution in the proton D on the r.h.s. of (A1)).

One parton splitting into two
For the perturbative correlation we have We split the factor (1 − z) into two pieces, (1) + (−z). The first one gives Here we have used the evolution equation for the second D function differentiated over the smaller scale k 2 : where S A is the Sudakov form factor of the parton A depending on the two scales, the overall Q 2 and the floating splitting scale k 2 [42,43].
An alternative evolution equation where the derivative is applied to the upper scale of the parton distribution in the proton reads This equation allows us to analogously represent the second piece (−z) as the derivative of the first D-function over the upper scale: with y ≡ zy. Combining (A2) and (A5), we get a full logarithmic scale derivative of the product of the Dfunctions: Once the two-and one-parton contributions (A1) and (A6) are taken together, the unphysical quantity y 1 cancels out, and we arrive at the desired sum rule (6). We conclude that for consistency of the DPI picture, the perturbative parton correlation (and thus the 1 ⊗ 2 subprocesses) should be taken into full consideration at that very moment when one allows distributions of partons picked from the hadron wave function to evolve with the hard scale(s).