Correlations in double parton distributions at small x

We present a dynamical study of the double parton distribution in impact parameter space, which enters into the double scattering cross section in hadronic collisions. This distribution is analogous to the generalized parton densities in momentum space. We use the Lund Dipole Cascade model, presented in earlier articles, which is based on BFKL evolution including essential higher order corrections and saturation effects. As result we find large correlation effects, which break the factorization of the double scattering process. At small transverse separation we see the development of"hot spots", which become stronger with increasing Q^2. At smaller x-values the distribution widens, consistent with the shrinking of the diffractive peak in elastic scattering. The dependence on Q^2 is, however, significantly stronger than the dependence on x, which has implications for extrapolations to LHC, e.g. for results for underlying events associated with the production of new heavy particles.


Introduction
Analyses of high energy pp collisions show that multiple hard parton collisions are quite common. Four-jet events, in which the jets balance each other pairwise, were observed at √ s = 63 GeV by the AFS collaboration at the CERN ISR [1], and at the Tevatron events with four jets, or with three jets + γ, have been observed by the CDF [2,3] and D0 [4,5] collaborations. These results also show that the hard subcollisions are correlated, and double interactions occur with a larger probability than expected for uncorrelated hard interactions.
Knowing the cross section for double subcollisions is important for an estimate of the underlying event. A good understanding of the correlations is therefore also essential for the interpretation of signals for new physics at the LHC. The aim of this paper is to study what kind of dynamical correlations in momentum space and impact parameter space are to be expected at higher energies, as a result of parton evolution to small x and of saturation effects.
Correlations which follow from energy and flavour conservation have been discussed e.g. in refs. [6][7][8], and from an impact parameter picture in ref. [9]. Nontrivial correlations which are consequences of DGLAP evolution have been studied in refs. [7,10,11]. In ref. [12] it is shown that this implies an increasing correlation for higher Q 2 . However, in these references the connection between correlations in momentum and b-space have not been studied. Correlation effects also follow from fluctuations in the number of particles in the parton cascades, which is discussed in refs. [6,7]. For related problems in connection with production of two heavy gauge bosons, see e.g. refs. [13,14].
In the model by Sjöstrand and van Zijl [8], later modified in ref. [15] and implemented in the PYTHIA event generator [16][17][18], it is assumed that the dependence of the double parton density on the kinematic variables (x and Q 2 ) and the separation in impact parameter space (b) factorizes. The correlation is described by a distribution in the separation b, where a denser central region implies that central collisions have on average more hard interactions, and peripheral collisions fewer. The relative shape of this distribution is kept constant, but the width is scaled with energy, to accommodate the growth of the total (non-diffractive) and the multiple interaction cross sections with increasing energy. For the dependence on x and Q 2 effects of energy and flavour conservation are taken into account. In a recent modification of the model, the shape is also varying with x [19], in a way which gives smaller correlations for smaller x. At a fixed energy, higher Q 2 is related to larger x, which implies that the production of heavy mass objects is more common in high multiplicity events. Similar approaches are also implemented in the HERWIG++ event generator [20][21][22].
While it is straight forward to compare the models in PYTHIA and HERWIG with experimental data on double parton scattering, it is not trivial to disentangle exactly the sources of possible correlations in the underlying models. This makes it difficult to extrapolate to higher energy and/or higher Q 2 . Here we will instead use a detailed dynamical model for parton evolution, for a study of the correlations between gluons inside hadrons.
In a series of papers we have presented a model based on BFKL [23,24] evolution and saturation, which well reproduces data on total, elastic and diffractive cross sections in DIS and pp collisions [25][26][27][28][29]. It is based on Mueller's dipole cascade model [30][31][32], but includes also non-leading effects from energy conservation and running coupling, as well as confinement effects and saturation within the evolution. The model is implemented in a Monte Carlo (MC) program, which makes it easy to study in detail the types of correlations and fluctuations in x and Q 2 as well as in impact parameter space, which are consequences of the parton evolution.
The correlations come from two different sources. First the distribution in transverse separation will vary with energy and virtuality, as the impact parameter profile widens at higher energy, but at the same time the development of "hot spots" makes the effective double scattering profile more narrow. Secondly, BFKL dynamics implies large fluctuations in the cascade evolution. This implies that double scattering is more probable in events with a higher than average number of partons. This effect is neglected in most phenomenological analyses, which are based on average parton densities.
The outline of this paper is as follows. First we will go through experimental aspects and the theoretical framework for double parton scattering in section 2. In section 3 we will then briefly describe the dipole model in the DIPSY Monte Carlo, followed by a description in section 4 of how we will use this model to investigate correlations in the double parton densities. The results are presented in section 5 followed by conclusions and outlook in section 6.

Experimental results
A measure of the correlation is given by the quantity σ eff defined by the relation Here σ D (A,B) is the cross section for the two hard processes A and B, σ S A and σ S B the corresponding single inclusive cross sections, and (1 + δ AB ) −1 is a symmetry factor equal to 1/2 if A = B. If the hard interactions were uncorrelated, σ eff would be equal to the total non-diffractive cross section. The CDF and D0 measurements give instead σ eff ∼ 15 mb, which thus is significantly smaller.
As mentioned in the introduction, experimental results on double parton interactions have been published from the ISR and the Tevatron. (The UA2 collaboration at the CERN SppS collider has presented a lower limit for σ eff .) The ISR experiment [1] studied events with four jets with a summed transverse energy above 29 GeV in pp collisions at √ s = 63 GeV. The observed rate corresponds to σ eff = 5 mb, indicating a very strong correlation between the partons. As the sum of E ⊥ for the four jets is about 30 GeV, the typical x-values for the partons is about 0.25. These partons must be dominated by valence quarks, and therefore not representative for the partons involved at higher energies and smaller x.
The CDF collaboration studied 4-jet events with p ⊥jet > 140 GeV at √ s = 1.8 TeV [2], which implies x-values ∼ 0.04. The result obtained was σ eff = 12.1 +10.7 −5.4 mb. CDF also studied events with γ+3jets [3]. This gave a more clear signal than the 4-jet events, and the result is σ eff = 14.5 ± 1.7 +1.7 −2.3 mb. The photons had E ⊥γ > 16 GeV, and the jets E ⊥jet > 5 GeV, and thus the x-values are smaller. No dependence on the Feynman scaling variable for the two pairs was observed within the ranges 0.01 < x F (γ + jet) < 0.4 and 0.002 < x F (dijet) < 0.2. We note, however, that the 4-jet and γ+3jet processes are not equivalent, as the production of a photon are particularly sensitive to quark distributions.
While the D0 4-jet signal is rather weak [4], this collaboration also measures a clear signal for γ+3jets [5]. The transverse momenta, and thus the parton x-values, are larger than those in the CDF analysis; p ⊥γ > 60 GeV and p ⊥jet2 > 15 GeV. The result obtained is σ eff = 16.4 ± 0.3 ± 2.3 mb. D0 also studied the dependence on p ⊥ of the second jet. Here σ eff dropped by 24% when p ⊥ increased from 17.5 to 27.5 GeV, but this variation was fully within the errors of the measurement.

Formalism
Following ref. [11] we define the "double parton distribution" Γ ij (x 1 , x 2 , b; Q 2 1 , Q 2 2 ), which describes the inclusive density distribution for finding a parton of type i with energy fraction x 1 at scale Q 2 1 , together with a parton of type j with energy fraction x 2 at scale Q 2 2 , and with the two partons separated by a transverse distance b. The distributions Γ ij are via a Fourier transformation related to the "two-parton generalized parton distributions" in transverse momentum space, D(x 1 , x 2 , Q 2 1 , Q 2 2 , ∆), studied by Blok et al. [33]. Assuming factorisation of two hard subprocesses A and B, the cross section for double scattering is given by 1 Hereσ is the cross section for a parton-level subprocess. The approximation used in the event generators PYTHIA and HERWIG means that, apart from effects of energy and flavour conservation (which are quite important), Γ factorizes in the form Here D i denotes the single parton distribution for parton type i, and F (b; s) describes the distribution in the transverse separation between the two partons, and is assumed to be independent of x i and Q 2 i , and to be the same for quarks and gluons. (As mentioned in the introduction, in a recent option in PYTHIA, F depends also on x [19].) The width is adjusted by tuning the model to reproduce the total cross section and multiple interactions, and thus depends on the collision energy √ s, as indicated in eq. (2.3).
As mentioned in the introduction, the relation in eq. (2.3) is not consistent with DGLAP evolution. Assume that the distribution of two gluons satisfies the factorization relation in (2.3) for some specific values of x 1 , Q 2 1 , x 2 , Q 2 2 . Evolving to higher virtualities, gluon 1 and gluon 2 can form separate cascades by emitting new gluons. A pair of gluons can then be formed, either with one from each of these cascades, or by two gluons from the same cascade. The latter contribution will then necessarily break the factorization relation [11]. Gaunt and Stirling assume instead the weaker relation The relations in eq. (2.2) and eq. (2.3) or (2.4) imply that the effective cross section in (2.1) is determined by the relation (2.5) In the following we want to use the Lund Dipole Cascade model to study the correlations and fluctuations, which follow from the parton evolution in a proton, and see how well the approximations in eqs. (2.3) and (2.4) are satisfied. We note, however, that as our model is based on the BFKL evolution, it contains only gluons and should only be trusted at small x-values where the gluons dominate.
We define the distribution F (b; x 1 , x 2 , Q 2 1 , Q 2 2 ) by the relation Thus F is a density in transverse coordinate space b, which may depend on all four variables x 1 , x 2 , Q 2 1 , and Q 2 2 , and which contains all information about correlations between the two partons. In case e.g. some kind of "hot spots" develop for small x and/or large Q 2 , this will show up as an increase in F for small b-values. With this definition F is related to the double scattering cross section via eq. (2.1) and the relation with the constraint For hard subcollisions at midrapidity we have x 1 ≈ x ′ 1 and x 2 ≈ x ′ 2 , and thus recover eq. (2.5), with the difference that σ eff may now depend on the variables x i and Q 2 i (with Q 2 1 /x 2 1 = Q 2 2 /x 2 2 = s). We note, however, that for subcollisions away from midrapidity we get different x-values in the two F -distributions in eq. (2.7). This feature will be further discussed in sec. 5.

Correlations
There are two different sources for correlations between the partons, which both give contributions to σ eff .

Distribution in impact parameter space
More central collisions are expected to have more hard subcollisions than peripheral collisions. This feature causes correlations, which depend on the matter distribution within the colliding protons. To understand the effect of the matter distribution, we here study a simplified version of the model implemented in the PYTHIA event generator. A discussion of the results in our model, and a comparison with more refined versions of PYTHIA, are presented in sec. 5.
Assume that the parton density inside a proton is given by a Gaussian distribution ρ ∝ exp(−r 2 /a 2 ). Integrating over the longitudinal coordinate z gives a factor √ π a independent of the impact parameter b = x 2 + y 2 . Thus the density in impact parameter space is also given by The two parton density is then given by The F -distribution is normalized to 1, and has therefore a normalization constant equal to (π2a 2 ) −1 . This implies that In PYTHIA it is assumed that for two colliding protons, the average number of hard subcollisions is proportional to the overlap of the two distributions. For a collision at impact parameter b, this overlap is given by (2.12) We note that the integral in eq. (2.12) is exactly the same as the one determining F in eq. (2.10). The distribution F is normalized to 1, and we therefore find Since in PYTHIA the average number of hard subcollisions at a fixed impact parameter, n(b), is proportional to the overlap, O(b), it can be written as where the cross section for parton-parton collisions is approximated byσ δ (2) (b). The number of subcollisions is assumed to be given by a Poisson distribution with averagen(b), and a non-diffractive event is obtained if n = 0. This implies that the total non-diffractive cross section is given by (2.15) and the average number of subcollisions is In PYTHIAσ depends on a parameter p ⊥0 , which acts as a smooth cutoff for small p ⊥ in parton-parton scattering. We have therefore two adjustable parameters, a and p ⊥0 , which can be used to fit σ ND and n . A result from such a fit is presented in ref. [19], with the result a ≈ 0.48 fm, giving σ eff ≈ 28 mb, at the Tevatron and σ eff ≈ 34 mb at LHC.
These values could be regarded as a kind of benchmarks. If the evolution gives smaller regions with higher parton density, the width of the F -distribution will be wider, and σ eff will be smaller. We here note that, although the single Gaussian fit gives the proper number of subcollisions, the correlations appear to be too week. More successful tunes to data have a distribution ρ (and thus also O), which is given by a sum of two Gaussian. This enhances small and large b-values, and gives a somewhat smaller σ eff . In a recent modification to the model [19], the width of the density ρ is assumed to vary with x, which also has the effect of enhancing small and large b-values increasing the correlations and reducing σ eff .

Fluctuations in the parton cascade
Another source for correlations is coming from fluctuations within the cascades, and affects σ eff even if the double parton distribution factorizes as in eq. (2.3). This effect is discussed in refs. [6] and [7], but is normally neglected in phenomenological analyses. Double hard scattering is more likely in collisions with protons that have more than the average number of partons. A measure of this effect is given by the integral of F as follows: The cascade evolution is a random process, which can lead to different partonic states. We label the states by the parameter n, and the parton distribution in state n is denoted D n (x, Q 2 ). The probability to obtain this state is denoted P n , with P n = 1. For a fixed state n the density of partons in a small interval δx i around x i at resolution Q 2 i , is given by D n (x i , Q 2 i )δx i . As we have a single state, the existence of a parton at x 1 is fixed, and it is not possible to define a correlation with the existence of another parton at x 2 . Therefore . Averaging over all states n we thus get . (2.17) Here . . . denotes an average over the different cascades n. For the special case of Q 2 1 = Q 2 2 and x 1 = x 2 , we get Thus the integral of F equals 1 + V / D(x, Q 2 ) 2 , with the variance V equal to the square of the width of the distribution. We see that the fluctuations imply that the F -distribution is enhanced, which also implies a larger correlation and a smaller σ eff . When presenting our results in section 5, we will also include values for the integral of F , in order to facilitate the interpretation of the results. Note that in the Gaunt-Stirling approximation in eq. (2.4), F is normalized to 1, and all correlation effects are included in the function

Mueller's dipole cascade
Mueller's dipole cascade model [30][31][32] is a formulation of the leading logarithmic (LL) BFKL evolution in transverse coordinate space. Gluon radiation from the colour charge in a parent quark or gluon is screened by the accompanying anti-charge in the colour dipole. This suppresses emissions at large transverse separation, which corresponds to the suppression of small k ⊥ in BFKL. For a dipole (x, y) the probability per unit rapidity (Y ) for emission of a gluon at transverse position z is given by This emission implies that the dipole is split into two dipoles, which (in the large N c limit) emit new gluons independently. The result is a cascade, where the number of dipoles grows exponentially with Y .
In a high energy collision, the dipole cascades in the projectile and the target are evolved from their rest frames to the rapidities they will have in the specific Lorentz frame chosen for the analysis. The scattering probability between two elementary colour dipoles with coordinates (x i , y i ) and (x j , y j ) in the projectile and the target respectively, is given by 2f ij , where in the Born approximation The optical theorem then implies that the elastic amplitude for dipole i scattering off dipole j is given by f ij . Summing over i and j gives the one-pomeron elastic amplitude The growth in the number of dipoles also implies a strong growth for the interaction probability, but the total scattering probability is kept below 1 by the possibility to have multiple dipole-dipole subcollisions in a single event. In the eikonal approximation the unitarized elastic amplitude is given by the exponentiated expression 4) and the total, elastic, and diffractive cross sections are given by

The Lund dipole cascade model
In refs. [25,26,28] we describe a modification of Mueller's cascade model with the following features: • It includes essential NLL BFKL effects.
• It includes non-linear effects in the evolution.
• It includes effects of confinement.
The model also includes a simple model for the proton wavefunction, and is implemented in a Monte Carlo simulation program called DIPSY. As discussed in the cited references, the model is able to describe a wide range of observables in DIS and pp scattering, with very few parameters.

NLL effects
The NLL corrections to BFKL evolution have three major sources [34]: The running coupling: This is relatively easily included in a MC simulation process. The scale in the running coupling is chosen as the largest transverse momentum in the vertex [35].

Non-singular terms in the splitting function:
These terms suppress large z-values in the individual parton branchings, and prevent the daughter from being faster than her recoiling parent. Most of this effect is taken care of by including energy-momentum conservation in the evolution. This is effectively taken into account by associating a dipole of transverse size r with a transverse momentum k ⊥ = 1/r, and demanding conservation of the light-cone momentum p + in every step in the evolution. This gives an effective cutoff for small dipoles, which also eliminates the numerical problems encountered in the MC implementation by Mueller and Salam [36].

Projectile-target symmetry:
This is also called energy scale terms, and is essentially equivalent to the so called consistency constraint [37]. This effect is taken into account by conservation of both positive and negative light-cone momentum components, p + and p − . The treatment of these effects includes also effects beyond NLL, in a way similar to the treatment by Salam in ref. [34]. Therefore the power λ eff , determining the growth for small x, does not turn negative for large values of α s .

Non-linear effects and saturation
As mentioned above, dipole loops (or equivalently pomeron loops) are not included in Mueller's cascade model, if they occur within the evolution. They are only included if they are cut in the Lorentz frame used in the calculations, as a result of multiple scattering in this frame. The result is therefore not frame independent. (The situation is similar in the Color Glass Condensate [38][39][40] or the JIMWLK [41,42] equations.) As for dipole scattering the probability for such loops is given by α s , and therefore formally colour suppressed compared to dipole splitting, which is proportional toᾱ = N c α s /π. These loops are therefore related to the probability that two dipoles have the same colour. Two dipoles with the same colour form a quadrupole field. Such a field may be better approximated by two dipoles formed by the closest colour-anticolour charges. This corresponds to a recoupling of the colour dipole chains, favouring the formation of small dipoles. We call this process a dipole "swing". The swing gives rise to loops within the cascades, and makes the cross section frame independent to a good approximation. We note that a similar effect would also be obtained from gluon exchange between the two dipoles.
In the MC implementation each dipole is assigned one of N 2 C colour indices, and dipoles with the same colour index are allowed to recouple [26]. The weight for the recoupling is assumed to be proportional to (r 2 1 r 2 2 )/(r 2 3 r 2 4 ), where r 1 and r 2 are the sizes of the original dipoles and r 3 and r 4 are the sizes of the recoupled dipoles. Dipoles with the same colour are allowed to swing back and forth, which results in an equilibrium, where the smaller dipoles have a larger weight. We note that in this formulation the number of dipoles is not reduced, and the saturation effect is obtained because the smaller dipoles have smaller cross sections. Thus in an evolution in momentum space the swing would not correspond to an absorption of gluons below the saturation line k 2 ⊥ = Q 2 s (x); it would rather correspond to lifting the gluons to higher k ⊥ above this line. Although this mechanism does not give an explicitly frame independent result, MC simulations show that it is a very good approximation.

Confinement effects
Confinement effects are included via an effective gluon mass, which gives an exponential suppression for very large dipoles [27]. This prevents the proton from growing too fast in transverse size, and is also essential to satisfy Froisart's bound at high energies [43].

Initial dipole configurations
In DIS an initial photon is split into a qq pair, and for larger Q 2 the wavefunction for a virtual photon can be determined perturbatively. The internal structure of the proton is, however, governed by soft QCD, and is not possible to calculate perturbatively. In our model it is represented by an equilateral triangle formed by three dipoles, and with a radius of 3 GeV −1 ≈ 0.6 fm. The model should be used at low x, and when the system is evolved over a large rapidity range the observable results depend only weakly on the exact configuration of the initial dipoles, or whether the charges are treated as (anti)quarks or gluons.

Application to double parton distributions
In principle we could estimate the gluon density in the proton from the cross section for γ * p collisions. The photon would then be treated as a superposition of dipoles of varying sizes, with weights determined by QED. However, in order to more easily isolate the correlations and fluctuations in the proton, without the complications from the additional fluctuations in the photon wave functions, we prefer instead to calculate the cross section for a dipole with a fixed transverse size. The photon coupling to a qq pair contains a factor (z(1 − z)) 2−i K 2 i ( z(1 − z) Qr). Here K i are generalized Bessel functions, r is the transverse separation between the quark and the antiquark, z is the fractional energy taken by the quark, and the index i is 1 for transverse and 0 for longitudinal photons. Important contributions are obtained for z ∼ 0.5, and since the Bessel functions fall off exponentially when the argument is larger than 1, characteristic r-values are given by the relation r ≈ 2/Q.
To determine the parton distributions D(x, Q 2 ) we thus calculate the cross section for scattering of a single dipole of size r = 2/Q colliding with a proton. The projectile dipole has a large cross section when colliding with equally large or larger dipoles in the proton, while the interaction with smaller dipoles is suppressed. The dipoles in the proton are also connected to gluons with k ⊥ of the order 1/r, and therefore the D-distributions obtained correspond to the integrated gluon distributions.
The double parton density Γ(x 1 , x 2 , b; Q 2 1 , Q 2 2 ) is in the same way proportional to the cross section, σ (1,2) , for simultaneous scattering of two dipoles with size r i = 2/Q i , separated by a distance b, and colliding with a proton. Thus also Γ is defined as an integrated density. From these relations we conclude that the F -distribution is directly related to the ratio between the corresponding cross sections, and in an obvious notation we have (2) . In principle the double scattering cross section, and thus also F , depends in the dipole model on the three vectors, the dipole sizes r 1 , r 2 , and the distance between them b. This implies that the estimate of σ eff should be given by a 6-dimensional integral over these variables. For the results presented in the next section we have calculated F keeping the separation between two gluons constant equal to b, but averaging over the directions of the dipoles, r 1 and r 2 . To check this approximation we have also calculated the result obtained when keeping the distance between the centers of the dipole fixed. The results presented in sec. 5 show that the correlations obtained are insensitive to how the averages are taken.

Results
In this section we show some results relevant for pp collisions at √ s ≈ 1.5 and 15 TeV, qualitatively representing Tevatron and LHC energies. We calculate the cross sections for one or two dipoles against a proton, where the rapidity separation between the projectile and the target is given by y = ln(1/x i ) + ln(Q i /2). Thus for subcollisions at midrapidity, with y = ln √ s, we have x = Q/(2 √ s).

Subcollisions at midrapidity
In figs. 1 and 2 we show the Q 2 -dependence for x-values corresponding to central production at a fixed energy. Fig. 1 shows the result for √ s ≈ 1.5 TeV and three combinations 2 of Q 2 1 and Q 2 2 . The corresponding x-values are given by x = Q/(2 √ s). Fig. 2 shows similar results for √ s ≈ 15 TeV.
We see that the distribution is not well described by a Gaussian. Instead there is an almost exponential dependence in the range 0. distribution drops faster , and this effect is stronger for higher Q 2 . The position of the break is related to the size of the initial proton configuration in the model. For small b-values a spike is developing, growing stronger with increasing Q 2 . This is associated with a faster drop for larger b for high Q 2 . Results with one softer and one harder subcollision lie (as expected) in between the results for two soft or two hard collisions, but closer to the result for two softer subcollisions (not shown for √ s ≈ 15 TeV).
In fig. 3 we show the energy dependence by comparing the results for Q 2 = 10 3 GeV 2 at the two different energies. We see here that the peaks at small b are almost identical, but at the higher energy the distribution becomes a little wider, with a slightly larger tail out to large b-values.
The corresponding values for σ eff = ( d 2 b F (b) 2 ) −1 and d 2 b F (b) were calculated numerically, and are presented in table 1. The rather slow fall off for F at large b gives a large σ eff for small Q 2 , but the more narrow distributions for higher Q 2 imply a very strong Q 2 -dependence. At 15 TeV σ eff drops by a factor 2 when Q 2 is changed from 10 to 10 2 GeV 2 . For fixed Q 2 the wider distribution at higher energy, implies a slightly larger σ eff . The variation with √ s is, however, much weaker than the variation with Q 2 . We also note that the effect of fluctuations in the cascades, represented by the difference from 1 of d 2 b F , is of the order 5-10% (contributing 10-20% to σ eff ), and is smaller for large Q 2 .

Sub-collisions off midrapidity
As mentioned in sec. 2.2, for subcollisions away from midrapidity, the effective cross section is given by eq. (2.7), which contains a product of two F -distributions with different x-values. In fig. 4 we show results relevant for two subcollisions with Q 2 = 10 GeV 2 at 1.5 TeV. One collision is at rapidity zero, with   and x ′ 2 = 0.0001, which corresponds to a rapidity for the pair of jets (or a produced massive particle) given by y pair = ln(x 2 /x ′ 2 )/2 ≈ 2.3. The corresponding F -distributions are called F S (for small x 2 = 0.0001) and F L (for large x 2 = 0.01). We see that in the tail F S and F L lie on opposite sides of the distribution F central , which is the relevant one for two subcollisions at midrapidity (both with x i = x ′ i = 0.001). As σ eff is given by ( d 2 bF S F L ) −1 , this implies that although the F -distributions vary with x, σ eff is approximately independent of y pair . This result is further illustrated in fig. 5, which shows the ratio √ F S F L /F central . We see that this ratio is close to 1 except for small b-values, which have a low weight in the integral over d 2 b. As a consequence also the result for σ eff varies very little with rapidity, in this  (see table 1). Actually, for small shifts δy pair , the difference between the product F L F S and F 2 central must be only second order in δy pair . This must then also be the case for the corresponding values for σ eff .

Comparison with experiment
The result for σ eff is larger than the effective cross section observed for γ +3jet events at the Tevatron. We note, however, that these data depend on quark distributions, and not only on the gluon distributions. It would be interesting to try to estimate the difference, but this goes beyond the scope of this paper. For Q 2 = 10 3 GeV 2 at 1.5 TeV we find σ eff = 23 mb. Within the errors this is in agreement with the 4-jet results from CDF, which however, have large uncertainties. It would also be interesting to have, for comparison, the values for σ eff in the tuned versions of PYTHIA, not only for the less successful single Gaussian approximation used in section 2.3.1.
We also note that there are two observed features, which are well reproduced by our model: i. The D0 collaboration [5] has observed a reduction in σ eff for growing p ⊥ of the jet pair in γ + 3 jet-events. Although the errors are large, the central values for σ eff drop from 18.2 to 13.9 mb, when p jet2 T increases from 17.6 to 27.3 GeV. This agrees qualitatively with the variation with Q 2 found in our model.
ii. The CDF collaboration [3] observes in the same reaction, that σ eff is very insensitive to a change in the rapidity y pair (or equivalently x pair F ) for one of the subcollisions. As discussed in sec. 5.2 this is also consistent with our model. We note, however, that this result is mainly a consequence of the fact that the product F L F S is close to F 2 central , and does not neccesarily show that F does not vary with x.
These qualitative agreements with experimental data may indicate that, also if our model possibly underestimates the correlations, we may believe that the qualitative features, and variations with Q 2 and √ s, are correct. The indicated strong dependence on Q 2 should then be very important for the interpretation of multiple hard interactions at LHC. We here note a possible qualitative difference between our model and the model with x-dependent matter densities in ref. [19]. Our result shows a larger variation with Q 2 for fixed energy, than with x (or energy) for fixed Q 2 . At fixed energy a variation with Q 2 is equivalent to a variation with x = Q/ √ s. However, the two models may give different results when data at different energies are compared. We would like to study this question more in the future.

Comment on the definition of b
As mentioned in sec. 4, in the dipole model F depends in principle on three vectors: the size of the two dipoles (r 1 , r 2 ) and the separation between them. For the results presented above we have defined the separation as the distance between two gluons, and averaged over r 1 and r 2 . When calculating d 2 b F (b) this gives exactly the correct result. This is, however, not the case for the integral d 2 b F 2 (b), which determines σ eff . When the separation is of the same order of magnitude as the size of the dipoles, the averaging does not give the correct result. To estimate this error we have changed the definition of b, to the separation between the centers of the dipoles, keeping this fixed when averaging over angles for r 1 and r 2 . The ratio between the two different F -distributions is shown in fig. 6. In this example we have Q 2 1 = Q 2 2 = 10 GeV 2 and √ s ≈ 15 TeV. We can here see a difference when b is of the same order of magnitude as the dipoles, around 0.1 fm.
The contribution from small b-values is, however, suppressed in the integral over d 2 b, and as mentioned above difference is even smaller, as the region where the two b-definitions give different results moves towards smaller b-values. This effect can therefore be safely neglected.

Conclusions and outlook
The Lund Dipole Cascade model offers unique possibilities to study the evolution of gluons inside hadrons at small x. The formalism is based on BFKL evolution including essential higher-order corrections and saturation effects. By following the evolution, emission-byemission, in rapidity and in transverse position, we can investigate the correlations and fluctuations of the gluon distribution in great detail. We have here concentrated on the double parton distribution, which enters into the multiple parton scattering cross section in proton-proton collisions. The double parton distribution in transverse coordinate space is the analogy of the generalized parton density in momentum space. The correlations can be described by a distribution in impact parameter space, F (b; x 1 , x 2 , Q 2 1 , Q 2 2 ). In many applications this distribution is assumed to be independent of the energy fractions and scales of the partons, and the same for quarks and gluons. In the PYTHIA event generator the width of the distribution does vary with energy, and in a recent study it is assumed to vary with x. Earlier analyses have demonstrated that DGLAP evolution implies non-trivial correlations depending on x and Q 2 , but this analysis does not give information about the b-dependence. In our analysis we find that the two-parton correlation function F (b) depends in a non-trivial way on all the kinematic variables x 1 , x 2 , Q 2 1 , Q 2 2 . For subcollisions at midrapidity, F is directly connected to the "effective cross section" used in experimental analyses, via the relation d 2 b F 2 = σ −1 eff . The result that F depends on x i and Q 2 i implies, however, that for subcollisions away from y = 0, σ eff is determined by an integral of a product of two F -functions with different x-values, connected by the relation x i x ′ i = Q 2 /s. Therefore σ eff is very insensitive to the rapidity of a hard subcollision; when one of the two F -functions in the product is above, the other is below the F -distribution relevant for midrapidity subcollisions. An experimental result showing a weak dependence on rapidity is therefore not a proof that F does not depend on the scaling parameters x i . In our formalism we can also study event-by-event fluctuations in the density of partons, which are usually neglected in other analyses. Neglecting the fluctuations F satisfies the normalization condition d 2 b F = 1, but when taking them into account we obtain a larger value for this integral, and therefore larger two-parton correlations. In our analysis this effect contributes 10-20% to the value of σ eff .
By studying the correlation function F , we can see explicitly how the emission of high-p ⊥ gluons result in so-called hot spots developing for small b-values at large Q 2 . For fixed Q 2 the emission of larger dipoles, related to low-p ⊥ gluons, implies that the distribution widens as we go down in x. (This effect is related to the shrinking of the diffractive peak at higher energies.) As a result σ eff is reduced for higher Q 2 at fixed energy, but increases at higher energy for fixed Q 2 . We see, however, that the result appears to depend much stronger on Q 2 than on x or collision energy.
It is not straight-forward to compare our results for the effective cross section with data from the Tevatron and ISR, which correspond to lower values for σ eff and thus stronger correlations. The ISR data are obtained for very high x-values, and the more exact Tevatron results are for γ + 3 jet events, which necessarily involves quarks. Our calculations include only gluons, and should therefore only be applied at small x. It is also difficult to estimate a possible difference between quark and gluon distributions. Our results are, however, compatible within errors with the four-jet measurement at CDF, where gluons should be dominating, and they show qualitative agreements with the observed rapidity independence at CDF and the hints of strong Q 2 -dependence at D0. We also note that our results appear to be in qualitative agreement with results from PYTHIA tuned to multiple interactions and underlying events.
Even if it is possible that we overestimate the result for σ eff to some degree, we are more confident about the variation with x and Q 2 , which should have relevance for extrapolations to the LHC. Our model here predicts a very strong decrease of the effective cross section with increasing Q 2 (especially when both Q 2 1 and Q 2 2 increases), while the dependence on the total energy (∼ 1/x) for fixed Q 2 is rather weak. Also the dependence on the rapidity of the two subcollisions is quite weak, which follows as σ eff here is given by a product of two F -distribution, one which is smaller and one which is larger than the F -distribution relevant for midrapidity subcollisions.
When comparing to the recent option in PYTHIA, with an x-dependent impact parameter profile [19], our result may be quite similar if we study the dependence on Q 2 at a fixed cms energy √ s. In our model the correlations increase as a direct consequence of higher Q 2 , and in the PYTHIA model they increase because x is higher for large Q 2 at fixed energy. The difference between the two models might therefore show up in the extrapolation of Tevatron data to LHC energies. It should be noted, however, that the extraction of the effective cross sections from data is very hard, and it will still be difficult to compare future LHC results to the numbers we have presented here. For this reason it is important to compare to event generator predictions for the experimental observables. Fortunately we have now developed the DIPSY MC to also generate exclusive final states [44], which means that in the future we can in detail simulate exclusive final state observables related to double parton scattering.