Sum rule improved double parton distributions in position space

Models for double parton distributions that are realistic and consistent with theoretical constraints are crucial for a reliable description of double parton scattering. We show how an ansatz that has the correct behaviour in the limit of small transverse distance between the partons can be improved step by step, such as to fulfil the sum rules for double parton distributions with an accuracy around 10%.


Introduction
To analyse data taken at the Large Hadron Collider in the best possible way, it is of great importance to have sound theoretical control over the QCD dynamics of proton-proton collisions. The mechanism of double parton scattering (DPS), in which two partons in each proton participate in a hard-scattering process, can give important contributions to particular final states and in particular kinematic regions. A prominent example is the production of two W bosons with the same charge [1][2][3][4][5][6], a channel that is also a background in searches for new physics (see e.g. [7][8][9]). A variety of DPS processes have been studied experimentally at the LHC [6,[10][11][12][13][14][15][16][17][18][19][20] and at lower energies [21][22][23][24][25][26][27][28][29][30][31] (see e.g. figure 4 of [18] and figure 15 of [32] for overviews). Recent years have seen significant progress in the QCD description of double parton scattering, see e.g. [33][34][35][36][37][38][39][40][41] and the brief overview in [42]. In particular, the formalism developed in [38,[43][44][45][46] extends the factorisation proofs for single Drell-Yan production [47][48][49] to double parton scattering with colourless final-state particles and achieves a consistent combination of single and double parton scattering contributions to a given final state. The non-perturbative quantities in DPS factorisation formulae are double parton distributions (DPDs), which specify the joint distribution of two partons in a proton. In the formalism just mentioned, these distributions depend in particular on the spatial separation y of the two partons in the plane transverse to the proton momentum. Alternatively, one may work with the transverse momentum ∆ that is Fourier conjugate to y. Given the complexity of measuring and computing DPS cross sections, a largely modelindependent fit of DPDs to experimental data, akin to what is done for single parton distributions (PDFs), will not be possible for a considerable time. It is hence essential to develop realistic models for DPDs. Considerable efforts have been made to compute them in quark models [50][51][52][53][54][55][56][57][58][59][60][61], and lattice calculations of Mellin moments of DPDs are underway [62,63]. In addition, there are important theoretical constraints on DPDs. On the one hand, there is the perturbative splitting of one parton into two [33-38, 40, 41, 44, 56, 64-75], which determines the behaviour of DPDs at small y and likewise puts constraints on DPDs depending on ∆. On the other hand, there are sum rules [76,77], which involve DPDs integrated over y (or evaluated at ∆ = 0) and express the conservation of momentum and quark number. So far, only a small number of studies [2,76,[78][79][80] have used these sum rule to constrain DPDs, and it is the goal of the present paper to continue this line of work. Whereas the DPD models in [2,76,78,79] are formulated for DPDs at ∆ = 0, we work with DPDs in y space, because these are the quantities required for computing DPS cross sections in the formalism of [44]. This paper is organised as follows. In section 2, we recall the theory underlying our model construction, highlighting in particular the nontrivial relation between DPDs depending on y and those depending on ∆. The starting point for our DPD model, taken from [44], is described in section 3. In section 4 we give a few technical details about our numerical calculations. In section 5, we make a series of changes to our model DPDs, improving at each step their agreement with the sum rules. The scale dependence of our results is studied in section 6, before we conclude in section 7.

Theory
The model analysis in this paper is based on the theory for double parton distributions developed in [44]. Let us briefly present the most important results of that work for our context. Consider the distribution function F a 1 a 2 (x 1 , x 2 , y; µ) for finding two partons a 1 and a 2 in the proton. The momentum fractions of the partons are x 1 and x 2 , and y denotes their spatial separation in the transverse plane. At leading order (LO) in α s , the scale dependence of DPDs is given by evolution equations dF a 1 a 2 (x 1 , x 2 , y; µ) d log µ 2 = α s (µ) 2π with the same DGLAP splitting functions P ab (v) that govern the evolution of ordinary PDFs at LO. For simplicity, we take a common factorisation scale µ for both partons in the present work, but it is straightforward to use different scales µ 1 and µ 2 . The behaviour of F a 1 a 2 (x 1 , x 2 , y; µ) at small y = |y| is dominated by the perturbative splitting of a single parton a 0 into the observed partons a 1 and a 2 . Evaluating the splitting mechanism at LO in α s , one obtains F a 1 a 2 ,spl,pt (x 1 , x 2 , y; µ) in D = 4 − 2 dimensions, where f a 0 is the PDF for parton a 0 . The function P a 1 a 2 ,a 0 (v, ) is equal to the ordinary DGLAP splitting function P a 1 a 0 (v) for = 0, and its form for nonzero may be found in [81]. Instead of the DPDs F (x 1 , x 2 , y; µ) in transverse-position space, one may also consider distributions depending on the transverse momentum ∆ that is Fourier conjugate to y. Since according to (2) the distribution F (x 1 , x 2 , y; µ) behaves like 1/y 2−4 at short distances y, its Fourier transform w.r.t. y requires an additional renormalisation in the ultraviolet. One way to achieve this is to perform the Fourier transform in 2 − 2 transverse dimensions. This gives rise to a 1/ ultraviolet pole that can be renormalised using standard MS subtraction, after which one can set to zero. Owing to this additional renormalisation, the evolution equations of the resulting momentum space distributions F (x 1 , x 2 , ∆; µ) differ from those of F (x 1 , x 2 , y; µ) by an inhomogeneous term that can readily be deduced from (2). This inhomogeneous equation has long been known and discussed in the literature [64-67, 70, 76]. An alternative is to start from the distributions F (x 1 , x 2 , y; µ) in D = 4 physical dimensions and to cut off their 1/y 2 singularity at short distances in the Fourier transform: Here ν is a scale with dimension of mass, and Φ(u) is a suitable function, which may be taken as a hard cutoff where γ is the Euler-Mascheroni constant. This choice of b 0 is such that certain analytical expressions simplify, see [44,81].
Since the distributions in (3) differ from those defined with MS subtraction only by the treatment of the ultraviolet region, one can use the small y expression (2) to derive a perturbative matching equation between the two types of DPD: where we have abbreviated P (v, ) = ∂P (v, )/∂ and v = x 1 /(x 1 + x 2 ). Here Λ denotes a non-perturbative scale. It is understood that one should take ν ∼ µ to avoid logarithmically enhanced higher-order corrections . Under this condition, the ν dependence cancels between the first and second line of (5) within the stated accuracy. We will investigate this numerically in section 6.2. We remark in passing that the previous discussion can be extended beyond LO. The higherorder forms of (2) and (5) involve convolutions instead of ordinary products, and the NLO kernels for unpolarised partons have been computed in [81]. The distributions F (x 1 , x 2 , ∆; µ) are of particular interest because at the point ∆ = 0 they fulfil the sum rules formulated in [76].
and express the conservation of quark number and of momentum, respectively. Here F a 1 qv = F a 1 q −F a 1q denotes the valence combination for quark flavor q, and N qv is the number of valence quarks with flavour q in the target. Equivalent sum rules can be written down for DPDs integrated over x 1 , given the trivial symmetry relation F a 1 a 2 (x 1 , x 2 ; µ) = F a 2 a 1 (x 2 , x 1 ; µ).
Note that naively F (x 1 , x 2 ; µ) just corresponds to the integral of F (x 1 , x 2 , y; µ) over all y, as one would expect for a sum rule. As discussed above, this simple correspondence is however invalidated by the singular short-distance behaviour of the y dependent distributions. As shown in [77], it is indeed the momentum space DPDs defined with MS renormalisation and taken at ∆ = 0 that appear in the above sum rules (together with MS renormalised PDFs). Already in [76] it was pointed out that the inhomogeneous term in the evolution equations for momentum space DPDs is essential for ensuring that (6) and (7) are valid at all µ. The matching relation (5) allows us to devise models for the position space distributions F (x 1 , x 2 , y; µ), which are the primary quantities needed to compute cross sections in the formalism of [44] and at the same time to use the DPD sum rules (6) and (7) as constraints for these models. In practice, the sum rules will then only be fulfilled approximately and in a particular range of momentum fractions. This is the strategy adopted in the present work. One might think of a different procedure and start with a model for the momentum space DPDs F (x 1 , x 2 , ∆; µ), constructed such that the sum rules are satisfied exactly. Using the extension of (5) to arbitrary values of ∆, given in [81], one can then compute the functions F Φ (x 1 , x 2 , ∆; µ, ν). The latter can be used instead of F (x 1 , x 2 , y, µ) to compute the double parton scattering cross section, as shown in section 8 of [44]. This possibility shall not be pursued here. We note that it has proven to be difficult to devise a general ansatz for distributions F (x 1 , x 2 ; µ) that satisfy the sum rules exactly, with the only consistent solution so far being limited to the pure gluon sector [79]. Until further progress is made in that direction, the best one can achieve with either momentum or position space models is that the sum rules are satisfied approximately to a degree one deems satisfactory.

Initial model
As starting point of our work, we take the DPD model introduced in [44]. Let us briefly recall its features and motivation. We require that the DPDs have the small y behaviour given by the perturbative splitting mechanism at LO. This is achieved by using a two-component ansatz where F spl tends to the perturbative splitting form at small y, whilst F int remains finite in that limit. The µ dependence of both components is required to follow the evolution equations (2). The physical idea behind the separation (8) is that in F a 1 a 2 ,int the partons a 1 and a 2 originate from the "intrinsic" part of the proton wave function, whilst in F a 1 a 2 ,spl they are obtained from a parton a 0 in the proton by perturbative splitting. It should be borne in mind that this is meant to be a heuristic picture, rather than a distinction that could be formulated in a field theoretically rigorous way. For the intrinsic part of the DPD, we make an ansatz at the scale µ 0 = 1 GeV. It consists of the product of two PDFs with a factor for the y dependence and a "phase space factor" ρ a 1 a 2 suppressing the distributions close to the kinematic boundary x 1 + x 2 = 1, with ρ a 1 a 2 (x 1 , Apart from the factor ρ a 1 a 2 , this form is obtained if one assumes that the two partons a 1 and a 2 in the proton are completely uncorrelated. Under that assumption, one can express a DPD as a convolution of two impact-parameter dependent PDFs f a (x, b), cf. [82] and section 2.1 of [38]. If one furthermore assumes that the impact-parameter dependent PDFs can be expressed in terms of ordinary PDFs and a Gaussian impact parameter profile, i.e.
then the convolution integral in (11) yields a Gaussian with a width that is the sum of the single-particle widths, i.e. h a 1 a 2 = h a 1 + h a 2 . For the single-particle widths we use the values whose physical motivation is discussed in [44]. The phase space factor ρ a 1 a 2 ensures that the distributions go to zero when approaching the kinematical boundary x 1 + x 2 = 1. The first or second power of (1 − x 1 − x 2 ) is frequently used in the literature, but as observed in [76], this results in a strong violation of the sum rules in the region x 1 1. A much better agreement is obtained with a phase space factor that does not yield any suppression in that limit. This is achieved by dividing ( For the "splitting part" of the DPD, we make the ansatz where F a 1 a 2 ,spl,pt (x 1 , x 2 , y; µ y ) = 1 πy 2 is the splitting form (2) in D = 4 dimensions. As required by theory, the ansatz (14) tends to the perturbative result for small y, with power corrections of order y 2 /h a 1 a 2 . At large y, the 1/y 2 falloff of the perturbative result is dampened by the Gaussian factor in (14). For lack of better guidance, we take the same parameters h a 1 a 2 in this factor as in the intrinsic part (9). The splitting form (14) is evaluated at the scale with y max = 0.5 GeV −1 . In the perturbative regime y y max this corresponds to the natural choice µ ∼ 1/y, which avoids logarithmically enhanced corrections from higher orders. For large y, the scale µ y approaches a limiting value b 0 /y max ≈ 2.25 GeV, which ensures that neither α s nor the PDFs on the r.h.s. of (14) are evaluated at too small scales. For the parton densities appearing in both (9) and (15), we take the MSTW2008 LO distributions [83] with the small modification described in section 3.2 of [76]. The latter ensures that thed and thes PDFs are positive and thus admit a probability interpretation. For the strong coupling, we use the starting value α s (µ 0 ) = 0.682 adopted in the MSTW2008 LO analysis. Throughout this work, we fix the number of active quark flavours to n f = 3.

Technical implementation
With the general prescription (5) and the two-component model (8), the DPDs entering the sum rules are given by where the matching term follows from (5). In (17) we have used that the position space DPDs depend on y only via y. Whilst evaluating F match is straightforward, the numerical computation of the intrinsic and splitting terms is more demanding. In the following paragraphs, we give some details about our numerical implementation. A reader not interested in these technicalities may skip forward to section 5.
DPD evolution and grids. To evolve F int and F spl from their respective starting scales in (9) and (14) to the scale µ at which the sum rules are to be evaluated, we use a modified version of the code employed in the study [44], which was itself a modification of the original code described in [76]. With this code, we compute position space DPDs on grids in the momentum fractions x 1 and x 2 , the interparton distance y, and the renormalisation scale µ.
The momentum fraction grids are equidistant in the variables u i = log(x i /(1−x i )). We use 89 grid points in each x i direction, with the smallest and largest x i values being x min = 5 × 10 −5 and x max = 1 − x min . For the factorisation scale, we use 51 points on an equidistant grid in log µ 2 , with largest scale µ max = 172 GeV. For each grid point µ i , we define a grid point in y i such that µ i = µ y i with the function µ y given in (16). This is convenient for evaluating F spl at its starting scale. The smallest value µ min on the µ grid thus corresponds to the largest value on the y grid and is just slightly larger than the limiting value b 0 /y max ≈ 2.25 GeV of µ y for infinitely large y. It turns out that for evaluating the integrals in (17), the y grid just described is not quite dense enough at small y values, and that for F spl we also need additional points at large y.
Extending the y grid appropriately, we end up with 60 points for the intrinsic part and 90 points for the splitting part of the DPD.
Integration. At the starting scale µ 0 , the y dependence of the intrinsic part F int is given by a simple Gaussian factor. This does not remain true at other scales µ, because quark and gluon distributions mix under evolution and have different Gaussian widths in our model. However, we find that at the µ values we consider, the y dependence of the evolved distributions F int is reasonably well approximated by a linear superposition of Gaussians with widths h qq , h qg , and h gg . Determining the appropriate superposition by a fit for each pair (x 1 , x 2 ) on our grid, we can evaluate the first y integral in (17) analytically. This strategy does not work for the splitting part F spl , for which we perform the y integral numerically, using the values of the distribution on the grid in y. Finally, the integral over x 2 in the sum rules is evaluated numerically, using the values of F a 1 a 2 (x 1 , x 2 ; µ) on the x 2 grid. For both the y and the x 2 integrals, integration rules for equidistant grids with errors of order 1/N 4 are used if N > 4, where N is the number of grid points in the relevant integration interval.

Refining the model
In this section, we describe how the initial model described in section 3 is modified so as to fulfil the DPD sum rules to a good approximation over a wide range in x 1 . The modifications are performed in several steps, after each of which we quantify the degree to which the sum rules are satisfied. To this end, we follow [76] and consider the "sum rule ratios" with a 1 being a quark, an antiquark, or a gluon. Note that a number sum rule ratio cannot be defined for F ddv in this way, because the denominator of R a 1 qv is zero in that case. The same holds for F a 1 sv unless a 1 = s or a 1 =s.
Postponing the discussion of F ddv to the end of this section, we now take a closer look at the number sum rules involving s v . We first observe that the PDFs underlying our DPD model satisfy the relation f s (x) = fs(x), which is of course stable under LO evolution. As a consequence, our initial DPD model satisfies at all scales µ. This will remain true with the modifications made in the present section. One thus obtains R ssv = Rs sv and needs to consider only one of these ratios. Furthermore, the number sum rules for F a 1 sv with a 1 = s,s are satisfied exactly. To prove the relations (21), we first note that they hold separately for F int (x 1 , x 2 , y, µ 0 ) and for F spl (x 1 , x 2 , y, µ y ) in the model specified in section 3. It is easy to see that they are stable under LO evolution. Since they also hold for the matching term in (17), they are valid for the distributions F a 1 sv (x 1 , x 2 ; µ) entering the sum rules. We will separately evaluate the contributions of the three terms in (17) to the numerators of R a 1 qv and R a 1 , so as to see which part of the DPD model requires adjustment to improve a specific sum rule. We will show plots for selected sum rules that are representative of the general situation, or -when there are large differences between sum rules -show the best and worst cases. Throughout this section, we evaluate the distributions (17) for ν = µ = µ min , where µ min = 2.25 GeV is the smallest value on the grid described in the previous subsection. Other scale choices will be explored in section 6.

Zeroth iteration: Initial ansatz
Let us start with the DPD model described in section 3 and consider the momentum sum rules. They turn out to be satisfied surprisingly well, as is illustrated in figure 1. Notice that there is a rather large contribution from F spl to R g . This is readily explained by identifying which parton combinations can be produced by perturbative splitting at LO, namely qq, qg,qg, and gg, as well as channels obtained from those by interchanging the two partons. All 2n f + 1 DPDs appearing in the g momentum sum rule thus receive a sizeable splitting contribution. By contrast, for the u momentum sum rule shown in figure 1(a) there are just two parton combinations with a large splitting contribution, namely uū and ug. Another noteworthy feature is the relatively small size of the matching contribution, which is a consequence of our choice ν = µ.
Let us investigate at this point the phase space factor ρ a 1 a 2 in (9). In some of the earlier works on DPDs, a simple factor (1 − x 1 − x 2 ) has been suggested [84][85][86][87], whilst the more  (17), as well as the full result. The ±10% deviations from unity are indicated by a light grey band. Not shown is the separate contribution from the matching term F match , which is negligible in this case. The remaining plots in this section will follow the same conventions unless explicitly stated otherwise.
recent study [88] argued that a factor (1 − x 1 − x 2 ) 2 is more appropriate. Even higher powers n would be obtained if one were to generalise the Brodsky-Farrar quark counting rules [89,90] from PDFs to DPDs. Each of these variants leads to a very strong suppression of DPDs in the region where x 1 ≈ 1 and x 2 ≈ 0 (or vice versa), since in that region the suppression from the phase space factor comes on top of the suppression of the corresponding PDF. As discussed in section 3, it is more appropriate to divide (1 − x 1 − x 2 ) n by (1 − x 1 ) n (1 − x 2 ) n for a given n in order to remove the phase space suppression in the regions x 1 ≈ 0 and x 2 ≈ 0. Including this division, we have investigated the momentum sum rules for different values of n and find that best agreement is achieved for n = 2, as is illustrated by the comparison of figure 2 with figure 1(b).
Turning to the number sum rules, we find that these are violated quite strongly in the initial    model, as is illustrated in figure 3. The agreement does not improve with other choices of the power n just discussed. The adjustments discussed in the following will improve the situation considerably. Let us at this point note that the number sum rules for equal quark flavours (such as those in the lower row of figure 3) can receive a substantial contribution from g → qq splitting at the starting scale µ y of F spl . This contribution is negative for a 1 = q and positive for a 1 =q, given that F a 1 qv = F a 1 q − F a 1q .

First iteration: number effects and modified phase space factor
In the first iteration of our model, we implement the same two adjustments that were already made in [76]. To describe these adjustments, it is convenient to specify the ansatz (9) for F a 1 a 2 ,int with a 1 and a 2 taking the values q v ,q, g instead of q,q, g (with q v denoting the linear combination q −q). This switch from quarks and antiquarks to "valence" and "sea" quarks is familiar from the parametrisation of ordinary PDFs. Following the argumentation in [76], it is natural to change the ansatz for distributions with two valence quark labels so as to take into account "number effects", i.e. the fact that we have a finite number of valence quarks in the proton, two u and one d. To do this, we set F uvuv,int to half the value given by (9) and set F dvdv,int to zero. The latter corresponds to the simple intuition that the probability to find two "valence d quarks" in the proton is nil.
The second adjustment argued for in [76] is to modify the phase space factor from the parton independent form in (10) to with α(a) = 0.5 for a = q v 0 for a =q, g Whilst the original form in (10) satisfies 0 ≤ ρ a 1 a 2 ≤ 1, the phase space factor in (22) becomes greater than 1 when the momentum fraction of a valence parton tends to 0 and the momentum fraction of the other parton (valence or sea) tends to 1. Due to the PDFs in the ansatz (9), the intrinsic part of the DPD still goes to zero in that limit. With these modifications, we find that the momentum sum rules are further improved, such that for most of the x 1 range, relative deviations are less then 10%. This is illustrated in figure 4, which is to be compared with figure 1 for the initial model.  2) and the modified phase space factor given by (22) and (23). The corresponding plots for the original model are shown in figure 1.
A more significant improvement is obtained for the number sum rules, as can be seen from the comparison of figure 5 with figure 3. The modified phase space factor yields a weaker suppression for valence partons at large momentum fractions of the other parton. This largely mitigates the steep decrease of the sum rule ratios with x 1 in the initial model. Taking into account number effects strongly reduces the value of R uuv at low x 1 , which is much too high in figure 3(c).

Second iteration: parameter scan for the phase space factor
Given that there is no strong motivation to take the particular value 0.5 for α(u v ) and α(d v ) in (23), it is natural to explore whether tuning these parameters can improve the sum rule ratios further. We have therefore performed a parameter scan over these two powers. To quantify the degree to which the sum rules are fulfilled, we introduce  as a quality measure for each sum rule ratio R, where x min = 5 × 10 −5 . A global quality measure is then the sum δ gl of these measures over all sum rules, excluding of course the cases for which R a 1 qv cannot be defined, as specified below (20). Notice that in (24) we have taken an upper integration limit of x 1 = 0.8. This is because for very high x 1 , we consider even large relative deviations from the DPD sum rules to be acceptable: DPDs in this region are expected to be very small and should hence not play any role in cross sections that are of measurable size. The values of δ gl obtained in our parameter scan over α(u v ) and α(d v ) are shown in figure 6. A minimum is reached at which we take as the second iteration of our model. As illustrated in figure 7(a), the momentum sum rules are not strongly affected by this change of parameters. The same holds for number sum rules that do not involve u quarks, which is not surprising because α(u v ) has significantly changed whereas α(d v ) has not. Furthermore, we see in figure 7(b) that the change in R udv is very small. By contrast, all number sum rules for u v are significantly improved in the range x 1 ≤ 0.8, as illustrated in the lower plots of figure 7.
One may wonder whether tuning other parameters in our model can lead to further improvements. Candidates for such an endeavour are the parameter y max in the starting scale µ y of F spl , as well as the widths h a 1 a 2 of the Gaussian damping factor, which appears in both F int    (23) and (25) in the phase space factor. and F spl . We find, however, that changing these parameters does not lead to a significant decrease of δ gl , and that the minimum of δ gl is achieved for parameter values very close to those specified in section 3. We hence leave these parameters at their initial values. Notice that the Gaussian factor in the intrinsic part (9) of the DPD is normalised such that its integral over all y gives unity. Restricting this integral to y ≥ b 0 /ν has little effect, which explains why a change of h a 1 a 2 has almost no impact on the contribution of F int to the sum rules.

Third iteration: modifying the splitting part at large distances
After several modifications to the intrinsic part F int of our DPD model, we now turn to the splitting part F spl . While the latter can be computed for perturbatively small y, its form at large distance y needs to be modelled. We now modify the initial ansatz (14) and multiply F spl,pt by the superposition of two Gaussians in y, with a relative weight depending on the momentum fractions: The factor multiplying F spl,pt can be rewritten as the sum of two Gaussians, one multiplied with 1 − g a 1 a 2 (x 1 + x 2 ) and the other multiplied with g a 1 a 2 (x 1 + x 2 ). For the new width parameters h * a 1 a 2 we make the same ansatz as we did for h a 1 a 2 , i.e. we set h * a 1 a 2 = h * a 1 + h * a 2 . We take values such that the Gaussian factor exp −y 2 /(4h a 1 a 2 ) + y 2 /(4h * a 1 a 2 ) multiplying g a 1 a 2 is approximately the same for all parton combinations. Admittedly, the form (26) is rather special among all possible functions that have the correct limit at small y. Clearly, the requirement of fulfilling the sum rules is not nearly enough to determine the functional form of DPDs at large y, so that a particular ansatz must be made. Our choice has the feature of introducing a nontrivial interplay between the dependence on y and on the parton momentum fractions, controlled by a one-variable function g a 1 a 2 (x 1 + x 2 ) for each LO splitting process a 0 → a 1 a 2 . We will find that this is an adequate degree of complexity, in the sense that the sum rule constraints are sufficient to determine this function. Whilst strict positivity ofF spl requires g a 1 a 2 (x 1 + x 2 ) > 0, the procedure described below yields negative values of this function in some cases. We checked that the resulting full DPDs F int +F spl are still positive in the range of x 1 , x 2 and y covered by our DPD grids. Let us first consider the splitting g → qq, where q takes one of the values u, d, s. This splitting feeds into the number sum rules for equal quark flavours, which at this stage are least well satisfied. Judging the impact of the function g a 1 a 2 (x 1 + x 2 ) is complicated by the fact that the ansatz (26) forF spl is made at the y dependent scale µ y and needs to be evolved to the scale µ min where we evaluate the sum rules. For definiteness, we consider the sum rule (N qv + 1) fq(x 1 , µ min ) = 1−x 1 0 dx 2 d 2 y Fq qv,int (x 1 , x 2 , y, µ min ) +Fq qv,spl (x 1 , x 2 , y, µ min ) where here and in the following it is understood that the integrals over y are restricted to y ≥ b 0 /ν = b 0 /µ min . To simplify the determination of g a 1 a 2 (x 1 + x 2 ), we make two approximations. Firstly, we use that for small y the initial and modified splitting model do not differ significantly, i.e.
Combining both approximations gives where we use (29) below y sep and (30) above. Taking y sep = 1 GeV −1 ensures that (30) is rather well fulfilled, as µ min and µ y differ by at most 12%. We will find that |g qq | < 12, which corresponds to a relative discrepancy below 30% between the l.h.s. and the r.h.s. of (29). While this may not seem to be very precise, it will turn out to be sufficient for improving the sum rules significantly.
Using (26) and (31), the sum rule (28) can be approximated as where Fq qv denotes the second iteration of our model and we have abbreviated Here we used that at the scale µ y one has Fq qv,spl = Fq q,spl = F qq,spl and a corresponding relation forF a 1 a 2 ,spl . Shifting the integration variable on the r.h.s. of (32) from with We recognise in (34) a Volterra equation of the first kind [91]. We discretise this equation by taking both x 1 and x on the grid for DPDs discussed in section 4. The integral over x is turned into a sum using a simple trapezoidal rule in the variable u = log(x/ (1 − x)). The result is a linear system of equations with an upper diagonal matrix K ij , which is readily solved using Gauss-Jordan elimination. In order to have an analytic formulation for our model, we fit the obtained discrete values of g a 1 a 2 (x) to the form for each of the splittings g → uū, g → dd, and g → ss. This reproduces the general shape of the numerical results rather well, except for some deviations at very large x. The resulting functions are shown in figure 8(a) to 8(c), and the fitted parameters are given in table 1.
With these modified g → qq splittings, the agreement of the model with theqq v number sum rules improves significantly, as can be seen in figure 9(b) to 9(d). Remarkably, the modification of the g → uū splitting improves not only the sum rule forūu v but also one for uu v , as seen in figure 9(a). x (d) ggg Figure 8: Modification functions g a 1 a 2 (x) for the g → qq and g → gg splittings. For each channel we display the fit to the form (37) and the direct solution of the discretised Volterra equation (36). The direct solution is shown as a dashed curve with linear interpolation between each data point.  Table 1: Parameters of the modification functions g a 1 a 2 defined by (26) and (37). At this point, we recall that the ratio R a 1 qv is undefined for F ddv . In order to quantify how well the number sum rule for this distribution is satisfied, we introduce the modified ratiõ in which the zero prefactor in the denominator of (19) has been replaced with unity. The ratiõ R ddv should be close to zero. We see in figure 10 that this is indeed the case: the modification of the g → dd splitting improves not only the sum rule for Fd dv but also the one for F ddv . Altogether, we have reached a satisfactory agreement of our model with all number sum rules. The modification of the g → qq splitting also affects the quark momentum sum rules, as illustrated in figure 11. In the cases shown in the figure, the agreement of the momentum sum rule becomes slightly worse, whereas the changes in the remaining cases are insignificant. One could improve Rū and Rd by modifying the g → gū and g → gd splittings, but this would also affect the number sum rules ratios R guv and R gdv . We refrain from such an exercise, considering that the agreement shown in figure 11 is still satisfactory. The sum rule ratio that is farthest away from 1 after these improvements is the one for the gluon momentum sum rule. This can be adjusted by modifying the g → gg splitting at large y in the same way as discussed for g → qq. The parameters of the modification function g gg (x) are given in table 1, and the function itself is shown in figure 8(d). The resulting improvement of the sum rule can be seen in figure 12, and we have checked that none of the other sum rule ratios is adversely affected by this final modification of our model. (a) g momentum sum rule Figure 12: Change of the gluon momentum sum rule due to the modification of the g → gg splitting at large y.

Scale dependence
So far, we have evaluated the sum rules for DPDs and PDFs at the scale µ = µ min = 2.25 GeV, and with the matching between position and momentum space DPDs computed for a cutoff scale ν = µ. In this section, we investigate how the sum rules change if these scales are chosen differently.

Renormalisation scale
As shown in [76], the DPD sum rules are preserved under LO evolution. If they are approximately valid at some scale, one may expect that they are still approximately valid when the DPDs and PDFs are evolved to a different scale. We verified that this is indeed the case for the DPD model developed in the previous section. This is illustrated in figure 13 for momentum sum rules and in figure 14 for number sum rules. We evolved the distributions from µ min to µ = 144.6 GeV, which is a point on our µ grid. The DPD matching at the high scale is evaluated with ν = µ. In the case of the g momentum and the gu v number sum rule, we notice that the individual contributions from F int and F spl to the sum rule ratios change considerably under evolution, while the sum of all contributions remains nearly the same. This highlights the relevance of the perturbative splitting mechanism for ensuring the scale independence of the DPD sum rules, which was pointed out in a number of different studies [69,70,77].
In the figures for theūu v sum rule, we observe that the oscillatory behaviour of Rū uv , which is a consequence of the modified splitting term in our model, is less pronounced after evolution to µ = 144.6 GeV. This is a typical feature of scale evolution, which tends to "wash out" details of distributions when going from low to high scales.

Cutoff scale
The matching relation given in (5) is only accurate up to higher orders in α s and up to power corrections in Λ/ν. The higher order analysis in [81] reveals that the term of order α n s in the matching relation is accompanied by up to n powers of log(µ 2 /ν 2 ). Varying ν around its "natural value" µ thus provides an estimate of higher order and power corrections in the matching relation. Following a widespread practice for scale variations, we vary ν between µ/2 and 2µ, taking again µ = µ min . The resulting variation of the sum rule ratios for our final DPD model is illustrated in figure 15. We find the ν dependence to be moderate, with changes of 10% or less in the sum rule ratios in almost all cases. These variations are hence of the same order as the agreement of the sum rule ratios with 1. The theoretical uncertainties reflected by the ν variation also suggest that it is of limited value to tune the sum rule ratios obtained for ν = µ much further than we have done. The only sum rule ratio with a larger ν dependence is Rs sv , shown in figure 15(d), which varies up to 20%. To understand this, we note that the ν dependence of the splitting and matching terms in (17) is stronger than the ν dependence of the intrinsic term. The latter gives an important contribution to all sum rule ratios, except for Rs sv , where within our model it is strictly zero. One might wonder whether a change of the scale ν could systematically improve the agreement of our initial model with the sum rules. The examples in figure 16 show that this is not the case: the ν variation is not able to bring the ratios Rū or Rū uv close to 1 for all x 1 ≤ 0.8. We also note that the change of the sum rules with ν is roughly of the same size in our initial and final models. This justifies our choice of ν = µ for the tuning of the model described in the previous section.

Conclusions
The number and momentum sum rules for DPDs put important constraints on DPD parametrisations. We have shown that one can construct physically plausible models for DPDs in position space that approximately fulfil these constraints. Our starting point was the DPD ansatz used in [44], the construction of which ensures the correct small y limit given by LO perturbation theory, but does not take into account DPD sum rule constraints at all. That ansatz was then sequentially modified: we started by adapting the modifications discussed in [76] to our case and furthermore tuned some model parameters, using parameter scans and a measure that quantifies how well the sum rules are globally satisfied. In the last step, we modified the form of the parton splitting term at large y, where perturbation theory is not  applicable and this term has to be regarded as part of the non-perturbative model. Whilst the specific form of that modification was motivated more by practical considerations than by physical intuition, our exercise shows that one can adapt position space DPDs up to the point where all momentum and number sum rules are satisfied within about 10% accuracy. An exception to this statement is the region of parton momentum fractions x > 0.8, where even ordinary PDFs are poorly known and where double parton scattering processes will have tiny cross sections. We verified that the approximate validity of the sum rules remains stable under evolution from low to high scales. Furthermore, we find that the sum rules are robust under variation of the cutoff scale ν, which appears when converting DPDs from position to momentum space. The largest ν variation is observed for the number sum rule that involves only strange quarks, where we see effects of up to 20%. Since the ν variation reflects in particular the size of uncomputed higher orders in the parton splitting, and since we vary ν around a central value of 2.25 GeV, we find a scale variation of this size not too surprising. One can expect that the inclusion of perturbative splitting terms at NLO, which have been computed in [81], will improve the situation. Whilst perturbative calculations for double parton scattering have been pushed to higher orders in recent years, the construction of more reliable DPD models remains an outstanding task. The present work shows that two major theoretical constraints on DPDs, namely the small y limit and the sum rules (where y is integrated over) can be satisfied simultaneously at least in an approximate way. Of course, this theoretical input alone is not sufficient to pin down the DPDs, and ultimately, the predictions obtained with any DPD model should be compared with experiment. This will be a huge endeavour and must be left to future work.