Correlations in double parton distributions: effects of evolution

We numerically investigate the impact of scale evolution on double parton distributions, which are needed to compute multiple hard scattering processes. Assuming correlations between longitudinal and transverse variables or between the parton spins to be present at a low scale, we study how they are affected by evolution to higher scales, i.e. by repeated parton emission. We find that generically evolution tends to wash out correlations, but with a speed that may be slow or fast depending on kinematics and on the type of correlation. Nontrivial parton correlations may hence persist in double parton distributions at the high scales relevant for hard scattering processes.


Introduction
An intriguing aspect of proton-proton collisions at high energies is double parton scattering (DPS), where two partons from each proton interact in two separate hard subprocesses. While the description of single hard scattering has become an area of precision calculations, our understanding of double hard scattering (and of its extension to three or more hard subprocesses) remains fragmentary, both at the conceptual and at the quantitative level. The initial state of double parton scattering is quantified by double parton distributions (DPDs), which quantify the joint distribution of two partons in a proton, depending on their quantum numbers, their longitudinal momentum fractions and their relative transverse distance from each other. A better knowledge of these distributions is important because DPS processes can contribute to many final states of interest at the LHC [1][2][3] and, furthermore, because they quantify characteristic aspects of proton structure beyond the information contained in the familiar parton distribution functions (PDFs) for a single parton.
DPDs depend on a scale, which in a physical process is given by the typical scale of the hard scattering, just as for PDFs. The scale dependence of DPDs is described by a generalization of the familiar DGLAP evolution equations. Two versions of this have been discussed in the literature: a homogeneous equation describing the separate evolution of each of the two partons and an inhomogeneous equation, which has an additional term describing the perturbative splitting of one parent parton into the two partons that will undergo a hard scattering [4][5][6][7][8]. Which version is adequate for the description of double hard scattering processes remains controversial in the literature [9][10][11][12][13][14][15][16]. In this work we use the homogeneous equation, having in mind that a systematic theory of double hard scattering will treat the physics associated with the inhomogeneous splitting term separately.
The joint distribution of two partons in the proton can be subject to various correlations. A number of arguments suggest a nontrivial interplay between the dependence of DPDs on the longitudinal momentum fractions x 1 , x 2 of the partons, as well as between their momentum fractions and their relative transverse distance y [17]. Moreover, the polarizations of two partons can be correlated even in an unpolarized proton [12,18,19]. A recent study in the MIT bag model [20] finds indeed large quark spin correlations in its range of validity, i.e. at x i above, say 0.1, and at a low scale.
Given the large range of energy scales relevant in LHC processes, it is important to understand how correlations in the distribution of two partons evolve with the scale. One might for instance expect that spin correlations that exist at a low scale become diluted by the subsequent parton radiation that is described by DGLAP evolution to higher scales. The purpose of the present paper is to study the evolution behavior of different correlations in DPDs in a quantitative manner. First results of our study have been presented in [21]. This paper is organized as follows. In section 2 we recall some basics of DPDs and in section 3 some details of their scale evolution. The behavior of correlations between x 1 , x 2 and y under evolution is examined in section 4. A large part of our investigation, presented in section 5, is the study of how spin correlations evolve. In section 6 we investigate how correlations between x 1 and x 2 are generated by the phase space available for parton radiation, which is described at the technical level by the integration limits in the evolution equations. Our conclusions are given in section 7. In appendix A we motivate the choice of PDFs used in our studies, and in appendix B we collect some analytical expressions needed in section 5.

Double parton distributions
If one assumes factorization, then the cross section for a double parton scattering process can be written as dσ dx 1 dx 2 dx 3 dx 4 = 1 C p 1 ,p 2 ,p 3 ,p 4 σ p 1 p 3 (x 1 x 3 )σ p 2 p 4 (x 2 x 4 ) d 2 y F p 1 p 2 (x 1 , x 2 , y) F p 3 p 4 (x 3 , x 4 , y) + {color, flavor and fermion number interference terms} , where C is a combinatorial factor,σ p i p j is the tree-level cross section for the hard scattering initiated by partons p i and p j and F p i p j is a DPD for partons p i and p j in the proton. The formula can be extended to include radiative corrections forσ p i p j and then involves convolution integrals over parton momentum fractions, just as for single hard scattering.
There is no complete proof that factorization as in (2.1) actually holds, but many important ingredients to such a proof can be found in [12,22]. We shall not be concerned with the interference terms alluded to in (2.1), but note that the DPDs describing color or fermion number interference are accompanied by Sudakov double logarithms and have a different scale evolution than the one we are studying in this work. In the parlance of [12], the distributions F p i p j in (2.1) are color singlet DPDs. It is understood that the DPS cross section in (2.1) needs to be added to the one for single hard scattering. (One also needs to add the interference between single and double hard scattering, which has received only little attention in the literature so far and will not be discussed here). Double parton scattering can compete with the single scattering mechanism in parts of phase space and even in inclusive cross sections when the single parton scattering contribution is suppressed by multiple small coupling constants. In particular, DPS is enhanced for small momentum fractions x i , because the joint distribution of two small-x partons increases faster with 1/x than the distribution of a single one. As an immediate consequence, DPS tends to be more important at the LHC than at previous hadron colliders.
We can include the effects of parton spin correlations in (2.1) by denoting with p i the type of a parton and its polarization at the same time. F p i p j (σ p i p j ) are then sums or differences of DPDs (subprocess cross sections) for different polarization states of the partons. Following [12], we write q,q, g for unpolarized partons, ∆q, ∆q, ∆g for longitudinally polarized ones, δq or δq for transversely polarized quarks or antiquarks, and δg for linearly polarized gluons. For each label δq or δq the corresponding DPDs and hard-scattering cross sections carry one Lorentz index in the transverse plane (corresponding to the transverse polarization vector), whereas for each index δg we need two transverse indices. In DPDs for two quarks, the polarization combinations allowed by parity and time reversal invariance are [12] F qq (x 1 , x 2 , y) = f qq (x 1 , x 2 , y) , where we write y = y 2 and introduce the proton mass M so that all functions f have the same mass dimension. Furthermore, we useỹ j = ǫ jj ′ y j ′ with the antisymmetric symbol ǫ jj ′ in two dimensions, and As shown in [23], one can write for DPDs of one quark and one gluon, and for two gluons. Expressions analogous to (2.2) and (2.4) hold if one or two quarks are replaced by antiquarks.

Polarization effects in double parton scattering
DPDs for polarized partons contribute to the cross section (2.1) if the cross section differencesσ p i p j for the relevant hard subprocesses are nonzero. A systematic discussion would go beyond the scope of this work, but let us mention a few important examples. A detailed discussion of the impact of parton spin correlations on the production of two electroweak gauge bosons (γ * , Z or W ) has been given in [22,24]. One finds nonzero cross section differencesσ ∆q∆q , and for Z and W production alsoσ ∆qq andσ q∆q thanks to their parity violating couplings. If the corresponding DPDs for longitudinal quark and antiquark polarization are nonzero, this influences both the overall rate of DPS and the distribution in transverse momentum and rapidity of the leptons into which the gauge bosons decay. For γ * and Z production there is a nonzero cross section differenceσ δq δq , which leads to an azimuthal correlation between the decay planes of the two bosons, provided that the transverse polarizations of two quarks or antiquarks in the proton are correlated as well. The cross section differencesσ ∆a∆b for longitudinal polarization in jet production are nonzero for most combinations a, b of quarks, antiquarks and gluons, and the same holds for prompt photon production (see e.g. table 4.1 in [25]). This will impact the overall rate of DPS as well as the transverse-momentum and rapidity distributions of the jets if there are longitudinal spin correlations between two partons in the proton. For transverse (anti)quark polarization, there are nonzero cross section differencesσ δq δq andσ δq δq for jet production andσ δq δq for the prompt photon channel qq → gγ [26], which together with transverse polarization correlations in the proton induce azimuthal correlations between the relevant jet planes. Azimuthal correlations in the final state can also be induced by linearly polarized gluons, with a nonzero cross section differenceσ δg δg for jet production [27]. The cross section differenceσ gδg is zero in that case, but it is nonzero for the production gg → QQ of heavy quarks [28] and the production gg → γγ of a photon pair [29]. We thus see that a number of important DPD channels will be impacted by spin correlations of partons in the proton.

Evolution of double parton distributions
As discussed in the introduction, we consider the homogeneous evolution equation of DPDs. For two unpolarized quarks we then have is a convolution in the first or second argument of the DPDs with the splitting functions P ab known from the DGLAP evolution of single parton distributions. We use the leading-order (LO) approximation of the splitting functions throughout this work. The explicit evolution equations for all unpolarized and polarized DPDs, as well as a list of the associated splitting functions, are collected in appendix A of [23]. We note that the splitting functions for antiquarks are identical with those for quarks. Let us briefly recapitulate the pattern of evolution for small x i , starting with gluon distributions. For small argument x, we have where N c is the number of colors and C F = 2N c /(N 2 c − 1). In each case the correction to the asymptotic behavior is one power higher in x. For small x 1 (and x 2 not too large) the first convolution integral in (3.2) can receive a substantial contribution from the region x 1 ≪ z ≪ 1 − x 2 where the splitting functions take their asymptotic forms (3.3) while the DPDs are far away from the kinematic limit z = 1 − x 2 , where they become small. An analogous statement holds of course for the second integral in (3.2) at small x 2 (and not too large x 1 ). This explains the steep rise of the unpolarized gluon distribution with Q 2 . For longitudinal gluon polarization, this rise is weaker since P ∆g∆g and P ∆g∆q lack the 1/x singularity of their unpolarized counterparts.
For linearly polarized gluons (which do not mix with quarks under evolution) the situation is special. The small-x limit of the splitting function reads where n F is the number of active quark flavors. Here we have included the small-x limit of the NLO contribution computed in [30] because it has a 1/x enhancement while the LO term vanishes like x for x → 0. Unless this NLO effect is very large (i.e. unless one considers scales where α s is large) one hence expects that distributions for linearly polarized gluons have at most a moderate growth with Q 2 at small momentum fractions. For quark distributions the splitting functions P qq (x), P qg (x), P ∆q∆q (x) and P ∆q∆g (x) all tend towards constant values at small x, whereas P δqδq (x) → 2C F x. Compared with unpolarized gluons, one thus expects a much milder growth with Q 2 for quarks at small x i , irrespective of their polarization. The strongest increase is to be expected for unpolarized quarks since they mix with the large unpolarized gluon distribution.

Numerical implementation
To solve the evolution equations numerically, we use a modified version of the code described in [7,31]. The original code was written to solve the inhomogeneous evolution equations of Refs. [4,5] for unpolarized DPDs. We have modified the code by removing the inhomogeneous term and by adding the splitting functions for polarized partons.
The code solves the double DGLAP equations in a variable flavor number scheme. It works in x-space, on a grid in x 1 , x 2 and t = ln Q 2 , performing the evolution stepwise in t by a fourth-order Runge-Kutta method. The x i grid points are evenly spaced in log x i /(1 − x i ) , with an equal number of points in both directions. They are bounded from above by the kinematic limit x 1 +x 2 ≤ 1 and from below by the choice of x min = 10 −6 . The grid points in t are evenly spaced, ranging from t 0 to t max for which we chose different values in different parts of our investigation. We used 240 grid points in each of the x i directions and 60 grid points in t. The code is supplemented with a routine that interpolates between different grid points. We made some small changes to this routine, making in particular sure that it can handle polarized distributions, which may have zero crossings and negative values.
The accuracy of the original code was investigated in [7], with estimated errors below 1% for x i ≤ 0.3 and evolution from Q 2 = 1 GeV 2 to Q 2 = 10 4 GeV 2 . We have updated these estimates after our modification of the code and with our grid settings. We find again an accuracy better than 1% for the evolution of unpolarized DPDs, whereas for polarized distributions the accuracy increases up to 4% in some regions at moderate x i . In the vicinity of zero crossings, the relative accuracy diverges and is of course no longer a useful measure for numerical uncertainties.
For the solution of the evolution equations, the code transforms the DPDs to a straightforward generalization of the single parton "evolution basis". This basis is defined by [32] and similar combinations including heavier quarks, where q ± i = q i ±q i . Analogous linear combinations are formed for polarized partons. In this basis single parton evolution is particularly simple, since mixing only occurs between the singlet (Σ) and the gluon while the other combinations evolve separately. For the up and down quarks, V i corresponds to the valence contributions u v and d v . The evolution code makes use of the basis (3.5) for both partons. 4 Correlations between x 1 , x 2 and y Various studies of generalized parton distributions suggest a nontrivial interplay between the distribution of partons in longitudinal momentum and in transverse space [17]. Specifically, the impact parameter dependent single parton distribution f a (x, b), i.e. the probability density to find parton a with momentum fraction x at a transverse distance b from the proton center, is not simply the product between a function of x and a function of b. It is therefore natural to assume that there is also a correlation between the longitudinal variables x 1 , x 2 and the transverse distance y in DPDs. In this section we investigate how such a correlation behaves under scale evolution. We consider only unpolarized partons and focus on the region of small momentum fractions x i , which is relevant for a large range of DPS processes at the LHC.

Initial conditions
As a model for the DPD at the starting scale of evolution, we take the simple ansatz that follows if one assumes the two partons to be independent. The DPD can then be written as a convolution of impact parameter dependent single parton distributions f a (x, b). For these distributions, we assume a Gaussian b dependence with an x dependent width, namely Here f a (x) denotes the usual parton densities, for which we take the LO set of the MSTW 2008 analysis [33]. We return to the choice of this PDF set below. For the starting scale where the ansatz (4.2) is assumed we take Q 2 0 = 2 GeV 2 . We should note that the form (4.3) is tailored for the region of x up to 10 −1 and not meant to be realistic for larger x, see e.g. the discussion in section 7.3 of [34]. We take different parameters in (4.3) for gluons and for the sum q + = q +q and difference q − = q −q of quark and antiquark distributions, The parameter values for q − were obtained in a model dependent determination of generalized parton distributions from electromagnetic form factor data [34]. For the remaining parameters we use input from hard exclusive scattering processes. At leading order in α s , the scattering amplitudes for exclusive J/Ψ photoproduction and for deeply virtual Compton scattering (DVCS) are described by generalized parton distributions for gluons and for the sum q + of quarks and antiquarks, respectively. Up to an uncertainty from the so-called skewness effect, one can thus connect the measured t dependence in those processes with the Fourier transform of the distribution (4.2) to transverse-momentum space. The values of α ′ g and B g given in (4.4) have been determined in [35] to match the measurement of elastic J/Ψ photoproduction in [36]. Experimental uncertainties do not allow us to extract a value for α ′ q + from DVCS, and we take the same value as for gluons in this case. The value of B q + in (4.4) has been chosen to correspond to a DVCS cross section dσ/dt ∝ e bt with b ≈ 7 GeV −2 at x ≈ 10 −3 and Q 2 ≈ 2 GeV 2 , guided by a fit to the t dependence in [37].
Let us emphasize that the functional form and numerical values in (4.2) to (4.4) are not meant to be a precision extraction of impact parameter dependent parton densities, but as a simple ansatz in rough agreement with phenomenology. The focus or our interest is how correlations of this type are affected by scale evolution.
Inserting (4.2) into (4.1) we obtain our ansatz for the unpolarized DPDs, Note that the y dependence in (4.5) does not factorize into separate contributions from each of the two partons. In the following we will consider the combinations u − and u + as representatives of the quark sector for definiteness. To specify the mixing between gluon and quark singlet distributions, we take the parameters α ′ q + and B q + in (4.4) for all light quark flavors, u, d, s at the starting scale Q 0 . The charm distribution is negligibly small there, since for the MSTW 2008 distribution we have m c ≈ Q 0 .

Change under evolution
According to (3.1) DPDs evolve independently at each value of the interparton distance y. However, the nontrivial interplay between y and the momentum fractions x 1 and x 2 in the starting conditions (4.5) has consequences for the scale evolution at different values of y. The exponential factor in (4.5) leads to a suppression of the large x i region, which becomes more important as y increases. Furthermore, the relative size of gluon and q + distributions, which mix under evolution, changes with y because our ansatz implies h g (x) < h q + (x) and thus has a broader y profile for quarks than for gluons. Let us first see to which extent the Gaussian y dependence of our starting condition (4.5) is changed by evolution. To this end, figure 1(a) shows ln F u − u − as a function of y 2 for different values of the scale. A Gaussian y dependence of the DPD translates into a straight line in this plot. We see that the shape remains approximately Gaussian even up to the high scale of Q 2 = 10 4 GeV 2 , and that the slope in y 2 becomes steeper with Q 2 , corresponding to a narrowing of the Gaussian profile. This is seen more quantitatively in figure 1(b), where we show the slope of the curves as a function of y 2 , multiplied with an overall minus sign. The departure from a Gaussian y dependence after evolution is reflected in a slow decrease of the slope with y, but overall the effect is rather mild.
The corresponding plots for the DPDs for two u + or two gluons are shown in figure 1(c) to (f). For two u + we observe a similar trend as for two u − , with a slight departure from a Gaussian behavior and an overall narrowing of the y profile at higher scales. For two gluons, the y dependence also remains approximately Gaussian after evolution, but the local y slope shown in figure 1(f) shows a different behavior than for quarks, with a tiny increase at low y and a weak decrease at higher y. The overall size of the effect is, however, quite small.
Let us now take a closer look at the evolution of the width of the y dependence. We quantify this by taking the difference quotient between y = 0 and y = 0.4 fm, a region where according to figure 1 the approximation of a linear y 2 dependence of ln F aa (x, x, y) works very well. The function h eff aa (x, x) thus defined may be regarded as an effective Gaussian width. Figure 2(a) shows the evolution of h eff aa (x, x) at x = 0.01 for a = u − , u + and g. The width for the double u − distribution starts at a larger value than for the other partons, while the starting value of h eff gg is the smallest and h eff u + u + is found almost half way in between. As we already saw in figure 1, the effective Gaussian width decreases under evolution for both u − and u + , whereas it barely changes for the gluon. We also observe that h eff u + u + approaches h eff gg with increasing scale, although it does so rather slowly. The difference between the transverse distribution of gluons and quarks, which we have assumed at Q 0 , thus persists over a wide range of scales.
The dependence of h eff aa (x, x) on x is shown in figure 3(a), (b) and (c) for the different parton types. We see that evolution is faster at small momentum fractions x than at large ones. At low x, there is a rapid decrease of h eff aa (x, x) with Q 2 for all parton types. For u + this results in a region of intermediate increases with x at high Q 2 . For u − and g the curves for h eff aa (x, x) are approximately linear in ln(x) as long as we stay away from the large-x region. This allows us to extract an effective shrinkage parameter α ′ eff a by fitting the effective Gaussian width to in an appropriate region of x, which we choose as 0.004 ≤ x ≤ 0.04. At the starting scale of evolution, we recover of course the value of α ′ a in our original ansatz (4.6) for h aa (x 1 , x 2 ) at a is shown in figure 2(b). We find that α ′ eff a decreases quite rapidly for a = g and more gently for a = u − .
Turning our attention to the double u + distribution, we see in figure 3 that at large Q 2 a local fit of h eff u + u + (x, x) to the form (4.8) would result in a negative α ′ eff u − whose value would strongly depend on the chosen range of x. To see whether this is a particular feature of DPDs, we have investigated the evolution of the impact-parameter dependent single parton distributions f a (x, b), which proceeds according to the usual DGLAP equations at each value of b. Using the program QCDNUM [38] to perform the evolution, we find that the initial Gaussian b dependence in our ansatz (4.2) approximately persists at higher scales, so that we can extract an effective Gaussian width h eff a (x) from the difference quotient of ln f a (x, b) between b = 0 and b = 0.4 fm in full analogy to (4.7). The scale dependence of h eff u + (x) determined in this way is shown in figure 3(d) and shows the same qualitative behavior as h eff u + u + (x, x). We conclude that an increase with x of the effective Gaussian width is not special to the evolution of DPDs. It is likely connected to the strong mixing of the u + distributions with the large unpolarized gluon at small x.
Our studies described so far have been done with the MSTW 2008 parton distribution in the initial conditions, and it is natural to ask how much our findings depend on this choice. In appendix A we show a selection of recent LO PDF sets at the scale Q 2 0 = 1 GeV 2 and find that among the sets suitable for our purposes (namely those that are positive at that scale) the parameterizations of MSTW 2008 and GJR 08 [39] represent two extreme choices, with a very slow or a very fast increase of the gluon at small x, respectively. We have therefore repeated the studies reported in this section by replacing MSTW 2008 with GJR 08 in our ansatz (4.2) at the scale Q 2 0 = 2 GeV 2 . We obtain similar results, regarding both the qualitative effects of evolution (including the behavior in figure 3(c) and (d)) and the rate of change of the parameters describing the y dependence of the DPDs.

Evolution of polarized double parton distributions
We now investigate the evolution of spin correlations between two partons inside a proton. For simplicity we assume in this section a multiplicative y dependence of the DPDs, which is stable under scale evolution. Since our focus is on the degree of parton polarization rather than on the absolute size of the DPDs, we set the y dependent factor G(y) = 1 in all plots and omit the tilde inf p 1 p 2 (x 1 , x 2 ; Q) from now on.
For the unpolarized DPDs we content ourselves with a simple factorizing ansatz at the starting scale, which we take as Q 2 0 = 1 GeV 2 unless specified otherwise. For the single parton densities in (5.2) we consider the two LO sets MSTW 2008 and GJR 08, hereafter referred to as MSTW and GJR for brevity. This ansatz is clearly unsatisfactory close to the kinematic limit x 1 + x 2 = 1, where the DPDs are expected to vanish, but since we are not particularly interested in that region we have refrained from taking a more sophisticated form. As discussed in section 6, evolution to higher scales approximately conserves the factorized form (5.2) for sufficiently low x 1 and x 2 .
In figures 4 and 5 we show our unpolarized model DPDs for a uū pair and for two gluons, either as functions of x 1 at x 2 = x 1 or as functions of ln(x 1 /x 2 ) at fixed x 1 x 2 = 10 4 . Notice that ln(x 1 /x 2 ) is related to the rapidity difference between the systems of particles produced in the two hard-scattering subprocesses. Denoting the total four-momenta of the final states in the two subprocesses by q 1 and q 2 , we have For equal c.m. energies of the two subprocesses we simply have ∆Y = ln(x 1 /x 2 ). In that case, rapidity differences in the range −4 ≤ ∆Y ≤ 4 correspond to longitudinal momentum fractions 0.0014 ≤ x i ≤ 0.074 for x 1 x 2 = 10 −4 . In figure 4 we see that at low scales the double gluon distributions constructed from MSTW and GJR PDFs strongly differ in size and shape. Only for larger scales do the two sets approach each other, which is plausible since at high scales the single parton densities are more directly constrained by data than at low scales.
To model the polarized DPDs is much more difficult. There is no reason to believe that the single parton distributions for a polarized parton in a polarized proton should be suitable even as a starting point to describe the DPDs for two polarized partons in an unpolarized proton. In other words, it is far from obvious how to connect the spin correlations between one parton and the proton with the spin correlations between two partons, and we will not try to do so.
Instead, we pursue two scenarios. In the first one, which we call the "max scenario", we make use of the positivity bounds for DPDs derived in [23]. At the starting scale Q 0 of evolution, we maximize each polarized DPD individually with respect to its unpolarized counterpart. For the combinations we will investigate, this gives and (yM ) 2 |f aδg | ≤ f ag (5.5) for a, b = q,q, g at equal values of x 1 , x 2 and y on the left-and right-hand sides. As follows from equation (4.6) in [23], the bounds in (5.4) can be satisfied simultaneously, as well as  the bounds in (5.5), but not the two sets together. The distribution f t δaδb is subject to the same bound as f δaδb in (5.4). Since it also follows the same evolution equation and does not mix with any other distribution, we will not discuss it further.
As shown in [23], leading-order evolution to higher scales preserves the above bounds. By contrast, if the bounds are saturated at some scale, they will in general be violated at lower scales. For this reason, we take a rather low value of Q 0 in this study. Evolved to high scales, results in the max scenario show how large polarization effects can possibly be if one assumes that the density interpretation of DPDs and thus their positivity holds down to the scale Q 0 .
The sign of the distributions on the l.h.s. of (5.4) and (5.5) can be either positive or negative. For DPDs involving only transverse or linear polarization, such as f δqδq or f δgδg , this is of no consequence for the evolution behavior and we take the positive sign for definiteness. For polarized DPDs that mix with others under evolution, relative signs are important. In the following, we will always assume the positive sign for all polarized distributions in the max scenario. With other choices, one typically obtains lower polarization after evolution.
Our second scenario, called the "splitting scenario", contains more detailed dynamical input. At small distances y, DPDs can be calculated perturbatively in terms of single parton distributions as discussed in [12]. We then have and analogous relations for the other DPDs, which are collected in appendix B. At large y these relations will no longer hold. In the splitting scenario, we assume that at scale Q 0 the ratio of polarized and unpolarized DPDs computed in the perturbative regime is valid up to large values of y, even if the form (5.6) is not. We thus continue to use the factorized form (5.1) and (5.2) for the unpolarized DPDs, while the polarized ones at scale Q 0 are given by and (yM ) 2 f q δg = 2z 1 + z 2 f qg , (yM ) 2 f g δg = z 2 z 2 +z 2 + z 2z2 f gg , (5.8) where z = x 1 /(x 1 + x 2 ) andz = 1 − z. Further non-zero distributions are obtained by interchanging parton labels and momentum fractions or by interchanging quarks with antiquarks; other combinations such as f qq or fq δq vanish at Q 0 in this scenario. We note that the appearance of the factors (yM ) 2 in (5.5) and (5.8) is a consequence of the factors y k y k ′ M 2 in the DPD definitions (2.4) and (2.5). In cross sections, the distributions f q δg and f g δg always appear multiplied with (yM ) 2 . For convenience we will set this factor equal to 1 in our plots.

Quark and antiquark distributions
We start our examination of spin correlations with the DPDs for longitudinally or transversely polarized quarks and antiquarks.
We will show a series of figures which are all using the MSTW parton distributions in the initial conditions, with curves for the three scales Q 2 = 1, 16 and 10 4 GeV 2 . In each figure, the upper row shows the polarized DPDs and the lower row shows the ratio between polarized and unpolarized DPDs for the same parton type. The latter ratio is restricted to the range from −1 to 1 according to (5.4) and will be called the "degree of polarization". It is the size of this ratio that indicates how important spin correlations are in the cross sections of DPS processes. As we did for unpolarized DPDs in figures 4 and 5, we show the polarized distributions both as functions of x 1 at x 1 = x 2 and as functions of ln(x 1 /x 2 ) for x 1 x 2 = 10 −4 . Results will be given both for the max scenario and for the splitting scenario in the initial conditions.
The distribution for longitudinally polarized up quarks and antiquarks in the max scenario is shown in figure 6. The polarized distribution f ∆u∆ū evolves very slowly and hardly changes with Q 2 . However, the degree of polarization decreases with the evolution scale. This is due to the increase of the unpolarized DPDs, which can be seen in figure 4. At low x i the degree of longitudinal polarization decreases rapidly with Q 2 , whereas at intermediate and larger x i values it does so rather slowly. We find a degree of polarization around 50% at Q 2 = 16 GeV 2 and above 20% at Q 2 = 10 4 GeV 2 for x 1 x 2 = 10 −4 and a wide range of ln(x 1 /x 2 ).
In the splitting scenario we have f ∆q∆q = −f qq for all x 1 and x 2 at the starting scale. In order to facilitate comparison with the max scenario, we multiply the polarized distribution with −1 in figure 7. The mixing of f ∆q∆q with distributions involving gluons and the zero starting value of all distributions where the quark and antiquark have different flavors induce some differences in evolution compared to the max scenario. We observe a small decrease of f ∆u∆ū with the evolution scale and a somewhat lower degree of polarization. At large momentum fractions, the distribution does not quite reach 100% polarization in the x i range shown in the figure. The dependence of the degree of polarization on ln(x 1 /x 2 ) is slightly tilted towards larger polarization when the quark has a bigger momentum fraction than the antiquark.
We now turn to transverse quark and antiquark polarization, which leads to characteristic azimuthal correlations in the final state of DPD processes [24]. Transversely polarized quarks or antiquarks do not mix with gluons under evolution, nor with quarks or antiquarks of different flavors. Figure 8 shows the DPD for transversely polarized up quarks and antiquarks in the max scenario. There is a slight decrease of the DPD with Q 2 over the entire x i range, but the suppression of the degree of polarization is mainly due to the increase in the unpolarized distributions. The evolution of the degree of polarization is similar to the case of longitudinal polarization in the max scenario, with a somewhat faster decrease. At intermediate and large x i values, the degree of polarization decreases slowly. For x 1 x 2 = 10 −4 it amounts to 40% at Q 2 = 16 GeV 2 and to 10% at Q 2 = 10 4 GeV 2 over a wide rapidity range.
In the splitting scenario we have maximal negative polarization f δqδq = −f qq at the starting scale for x 1 = x 2 , but the degree of polarization decreases when the two partons have different momentum fractions and tends to zero for both x 1 ≪ x 2 and x 1 ≫ x 2 . The dependence of the polarized DPD on ln(x 1 /x 2 ) slightly flattens under evolution, and its overall size at x 1 = x 2 evolves in a similar way as in the max scenario. The same holds for the degree of polarization.  The polarization for other combinations of light quarks and antiquarks is of similar size and shows a similar evolution behavior as for the case of a uū pair just presented. Generically, the evolved distributions have a slightly larger degree of polarization for quarks compared with antiquarks, and for up quarks compared with down quarks. Obvious exceptions in the splitting scenario are distributions that do not have a quark and an antiquark of equal flavor. These distributions start at zero. For transverse polarization they hence  remain zero at all scales, while for longitudinal polarization they become nonzero due to the mixing with other distributions. The ratio f ∆u∆u /f uu for example can reach a few percent in the intermediate x i region for equal momentum fractions after evolution. In summary, we find that both longitudinal and transverse polarization remains sizeable up to large scales for intermediate and large x i values, provided it is large at low scales. In particular the distributions for longitudinally polarized quarks, which enter linearly in electroweak cross sections, can thus play a significant role. In the small x i region, however, both longitudinal and transverse polarization are strongly suppressed at high scales.
We have repeated our study with the MSTW distributions replaced by those of GJR in the initial conditions. Naturally, this has some effect on the quark distributions, but the differences on the degree of polarization are comparably small and do not change the conclusions we have just drawn.

Gluon distributions
Gluons can be polarized longitudinally or linearly. As discussed in section 3 the unpolarized (single or double) gluon density increases rapidly at small momentum fractions due to the 1/x behavior of the gluon splitting kernel. The absence of this low-x enhancement in the polarized gluon splitting kernels lead us to expect that the degree of gluon polarization will vanish rapidly in the small x region. As can be seen in figure 10 for longitudinally polarized gluons, this is indeed the case. The distribution f ∆g∆g does increase with evolution scale, but at a much lower rate than f gg . Evolution thus quickly suppresses the degree of longitudinal gluon polarization in the small x i region. In figure 10 we see that in the max scenario with MSTW starting distributions this suppression stays rather constant between Q 2 = 4 GeV 2 and Q 2 = 10 4 GeV 2 . For this range of scales, the degree of longitudinal polarization is between 10% and 15% at x 1 x 2 = 10 −4 , with a very weak dependence on ln(x 1 /x 2 ). Our knowledge of the single gluon distribution at the low scale remains, however, quite poor as is documented in appendix A. As an alternative to the MSTW distributions used in figure 10, we show in figure 11 the corresponding results obtained with the GJR parton densities. The degree of polarization at Q 2 = 10 4 GeV 2 is nearly twice as large as for MSTW and amounts to almost 20% at x 1 x 2 = 10 −4 . At Q 2 = 16 GeV 2 the difference is even more striking, with a degree of polarization equal to 30%, a factor of three larger than for MSTW. To understand this difference, we recall from figures 4 and 5 that at high scales the unpolarized gluon DPDs obtained with the two PDF sets are relatively similar. The  difference in the degree of longitudinal polarization between the two cases is hence mainly due to the polarized DPDs. At the starting scale, these are much larger if we use the GJR set instead of MSTW, and this large difference persists after evolution.
In the splitting scenario the differences between the results obtained with the two PDF sets are similar to the differences we just described for the max scenario. We show in figure 12 the results obtained with the GJR set and note that the degree of polarization obtained with MSTW distributions is significantly smaller. At the starting scale, the degree of longitudinal polarization has a maximum of 78% for x 1 = x 2 in the splitting scenario and quickly decreases when the two gluons have different momentum fractions. Evolution decreases the degree of polarization in a similar manner as in the max scenario, leaving us with polarization around 10% for central rapidities and Q 2 between 16 and 10 4 GeV 2 .
Linearly polarized gluons give rise to azimuthal asymmetries in DPS cross sections, in a similar way as transversely polarized quarks and antiquarks. The effect of evolution on the distribution of two linearly polarized gluons in the max scenario is shown in figure 13. We see that even the polarized distribution f δgδg itself decreases with the scale. Together with the rapid increase of the unpolarized two-gluon DPD this results in a rapid decrease of the degree of linear polarization, especially at small x i . As in the case of longitudinal gluon polarization, using the MSTW distributions at the starting scale (not shown here) results in an even faster suppression. In that case the degree of polarization is tiny already at Q 2 = 16 GeV 2 . In the splitting scenario, the ratio f δgδg /f gg is at most 11% at the starting scale, which leads of course to even lower polarization after evolution. We must hence conclude that even in the most optimistic scenario shown in figure 13, the correlation between two linearly polarized gluons is quickly washed out by evolution and can only be appreciable at rather large x i or rather low scales.  Figure 13: Distribution for two linearly polarized gluons in the max scenario, using the GJR PDFs in the initial conditions. Color (line style) coding as in figure 10. Here and in the following, the LO expression of the evolution kernel P δgδg is used, except in figure 16.
So far we have only considered the case when both partons have the same type of polarization. There is however the possibility to have an unpolarized gluon and a linearly polarized one, whose polarization direction is correlated with the interparton distance y. In the max scenario, the corresponding DPD at x = x 1 = x 2 is well approximated by for not too large x. At the starting scale, this is trivial, and at higher scales it reflects the fact that double DGLAP evolution proceeds approximately independently for the two partons, as long as x 1 + x 2 is not close to 1. In the splitting scenario, (5.9) does not hold even at the starting scale. Figure 14 shows f g δg in the max scenario. The presence of one unpolarized gluon increases the distribution for small x i values and results in a significantly larger degree of polarization. However, for very large Q 2 linear polarization is still strongly suppressed at small x i . For unequal x i we observe that the degree of polarization is enhanced when the unpolarized gluon has the smaller momentum fraction, reaching 30% for ln(x 1 /x 2 ) = −4 even at the high scale Q 2 = 10 4 GeV 2 . A significant degree of polarization for one linearly polarized gluon is also found in the splitting scenario, as shown in figure 15. The main difference to the max scenario appears for unequal x i . According to (5.8) the splitting g → gg is such that for very asymmetric kinematics x 1 ≫ x 2 the slow gluon carries maximal linear polarization. As is seen in the figure, evolution weakens this trend, and at very high scales we find again a higher degree of polarization if x 1 < x 2 rather than Linear polarization: effect of NLO corrections As we noted in section 3, the evolution kernel for linearly polarized gluons has a qualitatively different small-x behavior at leading and next-to-leading order in α s . To examine how this impacts the fraction of linearly polarized gluons in DPDs at low x i , we have incorporated the leading low-x term in the NLO kernel as given in (3.4). The results in the splitting scenario with either MSTW or GJR input distributions are shown in figure 16. Although the NLO corrections increase the polarization as one may expect, they do not significantly change the overall picture. The NLO enhancement is largest for the case of MSTW distributions around x 1 = x 2 = 10 −3 , where the polarized DPD has a dip and hence the increased migration of partons from larger to smaller x i is most important.

Quark-gluon distributions
In the two previous subsections we have seen that in general gluon polarization is washed out under evolution at a faster pace than the polarization of quarks. In this section we will see that the corresponding decrease of polarization for quark-gluon distributions is in between the pure quark and gluon cases. For the unpolarized DPDs, our factorized ansatz (5.2) results in the approximate relation as long as x is not too large. The mixed quark-gluon distribution for longitudinal polarization is shown in figure  17 for the max scenario using GJR distributions as input. There degree of polarization remains large up to high scales except for very small x i , and in figure 17(b) we find more than 20% polarization for Q 2 = 10 4 GeV 2 and x 1 x 2 = 10 −4 . The corresponding results   for the splitting scenario are shown in figure 18. For x 1 = x 2 the degree of polarization is 60% at the starting scale and decreases moderately fast as long as x i is not too small. In asymmetric kinematics, the polarization remains sizeable up to high scales if x 1 ≫ x 2 , while the polarization decays rather quickly for x 1 ≪ x 2 . Let us finally discuss the distributions for an unpolarized quark and a linearly polarized gluon. Figure 19 shows that in the max scenario we have a moderate degree of polarization after evolution, in particular when x 2 > x 1 . In the splitting scenario, shown in figure 20, the starting conditions provide maximal linear polarization in the limit where the gluon momentum is soft (x 2 ≪ x 1 ), in analogy to what we observed earlier for f g δg . This trend is preserved by evolution up to moderately high scales.

Very low starting scale
The studies presented so far in this section have taken an initial scale of Q 2 0 = 1 GeV 2 for evolution. One may ask how our findings change if we assume the starting conditions of the max or the splitting scenario to hold at a much lower scale. This question can be addressed if we construct the initial conditions from the GJR PDFs, which are available down to Q 2 = 0.3 GeV 2 . A very low scale is also typically associated with quark models, which may be used to calculate unpolarized and polarized DPDs for two quarks [20,40,41]. We note that very large spin correlations were found in the bag model study [20].
In figure 21 we compare the DPDs for two longitudinally polarized up quarks obtained in the max scenario with the two different starting scales just mentioned. Worth noting is that the peak in the valence-like distribution f ∆u∆u (x, x) at Q 2 = 0.3 GeV 2 migrates only slowly to smaller x under evolution. In the case of two transversely polarized quarks (not shown here) the position of that peak stays at x ≥ 0.1 for all scales, and the degree of transverse polarization is smaller than the degree of longitudinal polarization after evolution. At Q 2 = 10 4 GeV 2 and x 1 = x 2 = 10 −2 the degree of polarization for transversely polarized up quarks is 15%, compared with 35% for longitudinal quark polarization. Comparing figures 21(a) and (b) we note that with the higher starting scale, f ∆u∆u (x, x) shows a moderate growth with Q 2 at small x, whereas with the low starting scale it barely evolves at all for x below 10 −2 . The degree of polarization shows little difference between the two cases for Q 2 ≥ 16 GeV 2 and is largely controlled by the rise of the unpolarized DPD with Q 2 . Figure 22 shows the distribution for two longitudinally polarized gluons obtained by  starting evolution in the max scenario with the lower or the higher starting scale. With Q 2 0 = 0.3 GeV 2 the degree of polarization becomes suppressed already at Q 2 = 1 GeV 2 and then remains rather stable. At low x i it becomes altogether negligible, in contrast to the case where we start evolution at Q 2 0 = 1 GeV 2 . While f gg increases dramatically at low x i when evolved from Q 2 0 = 0.3 GeV 2 , its polarized analog f ∆g∆g rises only slowly and never even reaches the size it has in the max scenario at Q 2 0 = 1 GeV 2 . Only for x i well above 10 −2 do we find a significant degree of gluon polarization in the scenario with a low starting scale. Comparing the evolution of f δgδg from Q 2 0 = 0.3 GeV 2 (not shown here) with the one of f ∆g∆g in figure 22(b), we find that f δgδg does not increase with Q 2 at small x i . The resulting degree of linear polarization is hence even smaller than the degree of longitudinal polarization.
Taking the initial conditions of the splitting scenario at Q 2 0 = 0.3 GeV 2 we find the same general pattern as in the max scenario: little change in the degree of polarization for quarks at Q 2 ≥ 16 GeV 2 and a very small degree of polarization for gluons at low x i .

Effect of the kinematic limit
Perhaps the most immediate difference between single and double DGLAP evolution is the maximal momentum fraction that can be carried by a parton. While the evolution equation for a single PDF involves an integral of f b (z; Q) over momentum fractions z all the way up to 1, the corresponding integration in the double DGLAP equation (3.1) is limited by momentum conservation to for the evolution of the parton with momentum fraction x 1 . It is obvious that the reduced integration limit has an impact for very large x i values, but through evolution the effect  Figure 23: The ratio R defined in (6.2), which quantifies the effect of the kinematic limit on parton radiation in the double DGLAP equation. The starting distribution at Q 2 = 1 GeV 2 is the product of two PDFs from the MSTW set.
can propagate down towards smaller x i . We investigate this effect by evolving a product of MSTW distributions, f ab (x 1 , , with the DPD evolution equations and comparing the result with the product of evolved single parton distributions. The ratio would be equal to 1 at all scales if the phase space effect just described was absent. Figure 23(a) shows R uū and R gg at x 2 = x 1 for three different values of Q 2 . We can see how the effect of the integration limit propagates down towards lower x i , especially for the gluon distribution. In figure 23(b) we show R as a function of x 1 for three different x 2 values at high Q 2 . Overall, the effect of the integration limit is not large, except at large x i . Indeed, at the kinematic limit x 1 + x 2 = 1 the DPD should vanish, which is ensured by evolution even though it is not satisfied for our oversimplified starting conditions. Our exercise illustrates that even if one assumes that a DPD factorizes into the product of single parton distributions at some scale, this factorization cannot strictly hold at higher scales.

Conclusions
Correlations between partons can have a large impact in double parton scattering processes, both on the overall cross section and on the distribution of particles in the final state. We have shown that the effect of double DGLAP evolution, where each of the two partons develops its own parton cascade, generally suppresses such correlations at higher scales. The strength of this suppression varies widely, with a rapid decrease of correlations in some cases and a slow decrease in others.
At a certain degree of accuracy, the dependence of DPDs on the transverse distance y between the two partons is expected to depend on the type of the partons and on their momentum fractions. We have studied the evolution of a y dependence motivated by the phenomenology of generalized parton distributions at the initial scale. We find that a Gaussian y dependence at the initial scale is approximately preserved under evolution, with a noticeable but relatively slow change of the effective Gaussian width. Despite the mixing between gluons and quarks in the singlet sector, the differences between their distributions persist up to high scales.
Spin correlations between two partons in the proton are described by polarized DPDs. Positivity constrains these to be at most as large as the unpolarized DPDs for the same parton types, a property that is preserved under evolution to higher scales. For the initial conditions of evolution, we have either assumed maximum polarization or a degree of polarization as it is obtained when the two partons originate from the perturbative splitting of a single, unpolarized one. In the latter scenario, the degree of polarization strongly depends on the ratio x 2 /x 1 of momentum fractions for certain parton combinations. We find that the DPDs for two longitudinally or two transversely polarized quarks decrease slowly with the evolution scale Q 2 or even remain approximately constant. Given the rise of the corresponding unpolarized DPD, the degree of longitudinal or transverse quark polarization shows a rather pronounced decrease with Q 2 . The DPD for two longitudinally polarized gluons rises slowly with the scale, as does the DPD for a longitudinally polarized gluon and a longitudinally polarized quark. However, due to the very rapid increase of unpolarized gluon distributions with Q 2 , the degree of longitudinal polarization decreases with the scale both for gg and qg distributions, in particular at small x. The distribution for two linearly polarized gluons decreases under evolution, and the associated degree of polarization quickly becomes negligible for x below 10 −2 . In certain processes like double charm production, the DPD for one unpolarized and one linearly polarized gluon is relevant. Except at high x, it increases with Q 2 and the corresponding degree of polarization decreases rather gently. Our quantitative results in the gluon sector depend strongly on the initial scale of evolution and on the gluon densities used in the initial conditions, which entails a much stronger model dependence than for quarks. Broadly speaking, we find that almost all polarization effects become small for x ≤ 10 −2 and Q 2 ≥ 10 4 GeV 2 , whereas for x above a few 10 −2 many polarization correlations remain sizeable even at high scales.
The phase space available for parton radiation in DPDs is reduced compared with the case of single parton densities. As a result, evolution does not conserve the factorization of DPDs into separate functions of the momentum fractions x 1 and x 2 if one assumes this property at a certain scale. Quantitatively, this effect is important only if at least one of the momentum fractions is of order 0.3 or larger.
In summary, we find that the effect of scale evolution on parton correlations is important and should be included in quantitative estimates. The assumption that parton radiation will quickly wash out correlations is true in a few cases but cannot serve as a general guideline. How this affects double parton scattering processes remains to be studied in future work.

A Choice of single parton densities
The models we use for the initial conditions of unpolarized DPDs in sections 4, 5 and 6 are constructed from products of ordinary single parton distributions. For this we use LO PDFs, so as to match the leading-order evolution we perform for the DPDs. We need these PDFs at low scales, down to Q 0 = 1 GeV, where different PDF sets significantly deviate from each other. This is clearly seen in figure 24, where we show LO PDFs from Alekhin (a02m lo) [42], CTEQ6 (cteq6ll) [43], GJR (GJR08lo) [39] and MSTW (MSTW2008lo) [33]. We also show the dedicated Monte Carlo PDF set MRSTMCal [44]; the other set of that study (MRST2007lomod) looks similar. All PDF values are generated using the LHAPDF interface [45]. Not included in the figure are the LO PDFs of NNPDF2.1 [46] since they are not available at Q = 1 GeV via LHAPDF. The CTEQ6 gluon distribution turns negative at low x and is hence not suited for our purpose (one of our two models for polarized DPDs in section 5 builds on the positivity of parton distributions). For the same reason we discard the dedicated Monte Carlo PDFs of CT09 [47] (not shown here), although they are set to zero below some value of x instead of going negative. The LO gluon distribution of HERAPDF1.5 [48] is very close to the one of CTEQ6 at Q = 1 GeV (and hence not shown in the figure). It also turns negative at low x.  Figure 24: Comparison of recent LO PDF sets at scale Q = 1 GeV for gluons (a) and for u quarks (b). Note the different y ranges in the two panels.
Among the positive gluon PDFs shown in the figure, the sets of GJR08 and MSTW2008 represent extremes in the sense of having a very steep or a very flat behavior over a wide x range. We chose these two sets for our investigations of DPDs, expecting that results obtained with different PDFs should approximately lie within the range covered by the two representatives we have selected. In the right panel of the figure we see that the spread of u quark distributions in the different PDF sets is notable, but not as large as for the gluons.
For the evolution of DPDs, we adjusted the values of α s and of the quark masses to those used by the two PDF sets we have selected. This is

B Double parton distributions from perturbative splitting
For small interparton distances y, or more precisely in the limit yΛ ≪ 1, where Λ is a typical hadronic scale, the dominant contribution to DPDs is given by the short-distance splitting of a single parton into two [12]. To leading order in α s , one can then express the DPD as the product of a usual PDF with an expression for the perturbative splitting; at higher orders the product turns into a convolution. From the results of section 5.2 in [12] we can readily extract the corresponding expressions for the collinear color-singlet DPDs we are studying in the present work. For unpolarized or doubly polarized DPDs, we have at leading order in α s fractions. All other unpolarized or doubly polarized distributions, including f t p 1 p 2 (x 1 , x 2 , y), do not receive any contribution from perturbative splitting at this accuracy.
DPDs with one polarized and one unpolarized parton arise by perturbative splitting only for linearly polarized gluons, Because f g δg and f q δg are multiplied by two vectors y in their definitions, the associated distributions F jj ′ g δg and F jj ′ q δg diverge like 1/y 2 for y → 0, just as the distributions in (B.1). The distributions f g δq and fq δq do not receive contributions from perturbative splitting due to the chiral invariance of massless QCD.
The above distributions (together with their counterparts that are zero at the considered accuracy) saturate several of the positivity bounds discussed in [23]. Indeed, with these distributions one obtains two vanishing eigenvalues in each of the spin-density matrices ρ for quark-antiquark, quark-gluon and gluon-gluon DPDs, given in eq. (3.2) to (3.6) of [23].