On sum rules for double and triple parton distribution functions and Pythia's model of multiple parton interactions

Multi-parton distributions in a proton, the nonperturbative quantities needed to make predictions for multiple scattering rates, are poorly constrained from theory and data and must be modelled. All Monte Carlo event generators that simulate multiple parton interactions (e.g. Pythia) contain such a model of multi-parton PDFs. One important theoretical constraint for the case of double parton distributions is provided by the so-called number and momentum sum rules. In this paper we investigate to what extent the double parton distribution functions used in the Pythia event generator obey these sum rules. We also derive the number and momentum sum rules for the triple parton distribution functions and discuss how one can use the Pythia code to construct triple parton distribution functions which approximately satisfy these sum rules.

Double parton distribution functions (dPDFs) are crucial ingredients of theoretical predictions for DPS, encoding the correlations between the two partons arising from a given hadron 1 . Unfortunately, due to the absence of experimental fits, the dPDFs remain unknown. Therefore, in phenomenological studies of the DPS phenomenon one still has to rely on certain model dependent assumptions about the shape of dPDFs and their relation to standard single parton distribution functions (sPDFs). There are different approaches to modelling dPDFs. To name some we mention the Light-Front approach [39][40][41][42][43][44], the "bag model" of a proton [45], the AdS/QCD method [46,47], chiral quark models [44] and the models based on the solution of so-called "double" DGLAP evolution equations [24,26,[48][49][50]. First-principles calculations of the lowest Mellin moments of double parton distributions have also now been performed in lattice QCD [51][52][53].
Apart from the aforementioned approaches to the dPDFs there exist a broad class of Monte Carlo (MC) models of multiple partonic interactions (MPI) as implemented in general purpose event generators, e.g. PYTHIA [54][55][56], HERWIG7 [57] and SHERPA [58]. The detailed description of MPI models being used in PYTHIA and HERWIG7 can be found in corresponding publications [59][60][61][62][63][64][65]. The MPI model of SHERPA is based upon the old MPI model of PYTHIA [60] with some modifications that are described in Ref. [58] 2 . These MPI models are also widely used to simulate DPS processes.
Despite the variety of the models of dPDFs available on the market, there are common criteria dictated by first principles QCD which, to some extent, have to be preserved by each model of dPDFs. One such criteria is that the dPDFs in which the two partons are probed at the same scale ('equal-scale dPDFs') should satisfy the sum rules originally proposed by Gaunt and Stirling (GS sum rules) [26,37,67,68]. In this paper we demonstrate how one can use the MPI model of PYTHIA to construct sets of equal-scale dPDFs which approximately obey the GS sum rules, and quantify the extent to which these sum rules are satisfied. The PYTHIA model always considers MPI as an ordered sequence of interactions, and thus naturally yields dPDFs that are not symmetric in the first and second parton arguments, even in the equal scale case. We find that these 'asymmetric' dPDFs satisfy the GS sum rules when integrating over the x fraction of the 'second struck' parton at the level of a few percent when the other x argument is in between 10 −6 and 0.8. Since the equal-scale dPDFs should in fact be symmetric in the first and second parton arguments, we then generate symmetrised versions of these dPDFs via a simple procedure. We find that these dPDFs obey the GS sum rules (when integrating over either x argument) at about 20% accuracy level except for x fractions bigger than 0.4. As a benchmark we use the GS09 [26] set of dPDFs which represents solutions of the inhomogeneous "double" DGLAP evolution equations for an input which was constructed to approximately respect the GS sum rules.
Whereas the dPDFs and the corresponding sum rules have been extensively studied in the literature, the case of triple parton distribution functions (tPDFs) has received rather less attention. The evolution of the triple parton densities was studied in Ref. [69] and theoretical aspects of correlations in MPI have been studied in Refs. [30,70,71], however, the available phenomenological studies of triple parton scattering (TPS) [72][73][74][75][76][77] rely on factorization of tPDFs into products of three sPDFs which essentially neglects correlations between partons 3 . Moreover, to the best of our knowledge, no set of tPDFs satisfying the tPDF evolution equations has so far been constructed. Therefore, in the absence of such "first principles" sets of tPDFs, it is useful to have a method to build tPDFs that account for correlations in x-space related to momentum and valence number constraints in an approximate way. In this work we demonstrate how one can use the PYTHIA code to accomplish this task. This paper is organized as follows: in Section 2 we provide a short description of the GS09 model and the model being used in the PYTHIA event generator and in Sections 3-4 we provide our findings. Namely, in Section 3 we demonstrate how one can use PYTHIA code to construct asymmetric and symmetric dPDFs and study how well they satisfy the GS sum rules as well as perform some toy phenomenological studies for the double Drell-Yan process. Then, in Section 4 we derive number and momentum sum rules for the case of TPS, and demonstrate how one can use the PYTHIA code to create tPDFs which approximately obey such sum rules. Finally, in Section 5 we give our conclusions and some closing remarks.

Double parton distribution functions and the GS09 dPDFs
In this subsection we introduce the theoretical definition of dPDFs, and discuss their key properties including the GS sum rules. We will also briefly describe the GS09 dPDF set [26], and discuss simplifying assumptions that have been used in the past in order to utilise dPDFs directly in DPS cross section predictions.
First of all, consider the total DPS cross section for the production of final states A and B in the collision of hadrons h A and h B : 4 σ DPS AB = 1 1 + δ AB j 1 , j 2 , j 3 ,j 4 |y|>max(Q −1 Following [50,81,82], we shall refer to the objects Γ j1 j2/h (x 1 , x 2 , y, Q 1 , Q 2 ) as the double parton distributions (DPDs) (they are also referred to as "two parton generalized parton distributions" or 2pGPDs in the literature [35,86,87]). Loosely speaking, these can be interpreted as a probability to find two partons j 1 and j 2 with longitudinal momentum fractions x 1 and x 2 separated by transverse distance y in a hadron h. A lower cut-off on the transverse separation |y| is included since the region of very small |y| values is more appropriately described by single parton scattering (for more details see Ref. [50]). Let us define the dPDFs as the integral of the DPD functions over y, as follows: The dPDFs in Eq. (2.2) taken at equal factorization scales (Q 1 = Q 2 = Q Λ QCD ) obey so called inhomogeneous double DGLAP (dDGLAP) evolution equations [88][89][90] which, at the leading logarithmic level, have the form dD j 1 j 2 (x 1 , x 2 , t) dt = α s (t) 2π 3) where f j (x, t) are "standard" collinear sPDFs and we have dropped the subscript h A in D j 1 j 2 /h A in order not to overload the notation 5 . Here we use a dimensionless evolution parameter t = log(Q 2 /Q 2 0 ) where Q 0 is an arbitrary reference scale. The first two terms on the right hand side of Eq. (2.3) involve standard DGLAP leading order splitting functions P j 1 →j 1 (x) and, despite the presence of two Bjorken-x'es, essentially have the same structure as the standard DGLAP splitting kernels. The last term, however, is new and is linked to the reduction of the cut-off in Eq. (2.2) when Q is increased. It is an inhomogeneous term that involves the standard sPDFs and depends on "1 → 2" splitting functions P j →j 1 j 2 (x) which are defined at the leading order in α s via where κ(i, j) is the only parton that can be emitted at leading order when parton i splits to j (e.g. if i = g and j = q then κ(i, j) =q) and P R i→j (x) are given by the real-emission parts of the DGLAP splitting functions, see Ref. [26]. Therefore, in the following, we refer to the third term on the right hand side of Eq. (2.3) as the "splitting" or 1 → 2 term. As we have mentioned before, the dPDFs in Eq. (2.3) are taken at the same values of the factorization scales. The unequal scale dPDFs can be obtained from the equal scale dPDFs via a single DGLAP evolution step as described in Ref. [26].
The dPDFs obey the following sum rules: These dPDF sum rules were first introduced by Gaunt and Stirling in Ref. [26], who also established that Eqs. (2.5), (2.6) are preserved by dDGLAP evolution provided that the dPDFs obey the sum rules at the initial evolution scale Q 0 (see also Ref. [37]). Thus, we refer to Eq. (2.5) and Eq. (2.6) as the GS sum rules. The fundamental character of the GS sum rules was recently demonstrated by Diehl, Plößl and Schäfer in Ref. [68] where these relations were proved in all orders of perturbation theory for bare dPDFs and for the dPDFs renormalized according to the MS scheme (see also Ref. [93]). The MS-renormalized dPDFs coincide with the dPDFs in Eq. (2.2) up to corrections suppressed by O(α s ) or O(Λ 2 QCD /Q 2 ) [68], such that the dPDFs in Eq. (2.2) also satisfy the GS sum rules up to these corrections. In Ref. [26], a set of dPDFs was constructed which approximately satisfies the GS sum rules at the initial scale Q 0 = 1 GeV. These starting conditions were then evolved according to Eq. (2.3) to higher scales, up to Q 2 = 10 9 GeV 2 and for x 1,2 ≥ 10 −6 . The dPDFs generated by this procedure are referred to as the GS09 set 6 .
If one assumes that Γ j1 j2/h (x 1 , x 2 , y, Q 1 , Q 2 ) can be written as a product of a longitudinal and transverse-dependent pieces, with the transverse piece being some smoothly varying function with a width on the order of the proton radius, then one can write the DPD in terms of the dPDF as follows: provided that F (y) is normalised as d 2 yF (y) = 1. Under Eq. (2.7), the DPS cross section in Eq. (2.1) may be simplified as follows: where σ eff = d 2 y F 2 (y) −1 and has a dimension of area.
Nowadays, it is known that the factorization of the DPDs into the dPDFs and a transverse-dependent piece cannot hold [37,50,86,87,[97][98][99][100]. Actually, the discussion above already indicates this -if Eq. (2.7) held, then the cut-off in Eq. (2.2) would only have a power-suppressed effect, and there would be no inhomogeneous term in Eq. (2.3) 7 . Nevertheless, when using the GS09 dPDF set in cross section calculations below we will use Eq. (2.8) -this is because we would like to make the most direct comparison against the predictions of PYTHIA, which uses a framework for DPS/MPI that is most closely aligned with Eq. (2.8).
It is important to remark that although dPDFs satisfying the GS sum rules do not directly appear in the DPS cross section Eq. (2.1), they are linked to the DPDs that do, via Eq. (2.2). Any procedure that yields dPDFs that (approximately) satisfy the sum rules is thus of interest and utility even under the theoretically rigorous approach associated with Eq. (2.1). It is also worth remarking that it is possible to adapt techniques used to construct dPDFs satisfying the sum rules to construct DPDs which have both the appropriate transverse dependence and which integrate to dPDFs satisfying the sum rules -in Ref. [102] 6 For some phenomenological applications of the GS09 dPDFs see, e.g. Refs. [94][95][96]. 7 A further indication of the inadequacy of Eq. (2.7) comes from experimental measurements of DPS. If we take F (y) in Eq. (2.7) to have a width of order of the proton radius, and further assume that the dPDFs in Eq. (2.8) factorise into a product of single PDFs, then we obtain DPS cross sections which are roughly half as large as those obtained experimentally (see e.g. [101]). Inclusion of perturbative transverse correlation effects yields DPS cross sections in better agreement with data [37]. the techniques used to construct the GS09 dPDFs were adapted, extended and improved to yield such a set of DPDs.

The PYTHIA model of multiple parton distribution functions
Now let us briefly sketch the way the PYTHIA event generator models n-parton distribution functions (nPDFs) and, more specifically, dPDFs. We remind the reader that these dPDFs (and nPDFs) are applied to compute cross sections according to a framework that resembles Eq. (2.8); the transverse dependence is assumed to essentially factor out. The generation of the first hard interaction is performed in the usual way, for example, for the case of massless 2 → 2 scattering we have dσ dp 2 After the successful generation of the first hard interaction, one has to generate other subsequent interactions by taking into account MPIs already generated in a given pp collision.
In the general MPI model of PYTHIA these interactions are generated at successively lower scales. There is also an option to generate DPS events specifically (controlled by the option "SecondHard") -here both orderings of the two processes are simulated (see the beginning of Section 3.3 for more details). To generate the additional interaction(s) PYTHIA dynamically modifies sPDFs in Eq. (2.9) according to a history of previous interactions. In the following, while describing PYTHIA's framework, we use subscripts "r" and "m" to distinguish between unmodified ("raw") and modified sPDFs. In the first model of MPI [59,60] PYTHIA could only preserve overall energy and momentum by 'squeezing' the PDFs for the nth interaction such that instead of lying between 0 and 1 in the momentum fraction x n , they lie in between 0 and X n ≡ x i is the total longitudinal momentum taken by previous interactions. This squeezing is achieved by rescaling the x-argument in the PDFs to be Since the x-argument in a PDF cannot exceed 1, this implies x n < X n < 1. Aside from this argument rescaling, one also multiplies the PDF by 1/X n . In the case n = 2 we thus have a dPDF defined by: where f r j 1 (x 1 , Q) is the unmodified sPDF and is the sPDF modified according to the first interaction involving parton j 1 with momentum fraction x 1 . Since the x-argument of f r j (x, Q) cannot exceed 1, this construction implies that the dPDF is only nonzero when x 1 + x 2 ≤ 1.
Note that for the second interaction we have: That is, the momentum of all partons after the first interaction is 1 − x 1 , as desired. In general, for the nth interaction we will have: (2.14) This approach, however, does not take into account changes in the quark content of a proton. For example, if one probes a "valence" d-quark in the first hard interaction then for the second interaction the number of "valence" d-quarks in a beam remnant should be reduced to zero. Furthermore, if one probes a d-quark sPDF at low-x it means that with a high probability the d-quark originates from a perturbative 1 → 2 splitting of a gluon into a quark-antiquark pair, which, in turn, implies that there is a left-over "companion" antiquark which should be included in a beam remnant, see Fig. 1. These effects, obviously, apply to all quark flavours.
Figure 1: Two different ways to generate (d g → d g) ⊗ d g →d g DPS process in the PYTHIA event generator: a) the initial state d-quark comes from a perturbative g → dd splitting process, b) the initial state d-quark comes from a raw d-quark sPDF.
In 2004 the MPI model of PYTHIA was significantly improved: n-parton PDFs of arbitrarily high order n were constructed such that changes in quark content due to higherscale interactions were taken into account. For the first interaction a standard sPDF is assumed to be valid, as fitted to data. The sPDF for the second interaction thereafter depends on the properties of the first interaction, and so on [61]. The starting point of this approach is to split sPDFs into valence and sea parts as where the flavour-antiflavour symmetry in the sea is assumed. If we denote the raw valence sPDF for flavour f by q f v0 (x, Q), then the valence sPDF after n interactions is given by where N f v0 is an original number of valence quarks of a flavour f , q f vn and N f vn are distribution and number of valence quarks of the flavour f after n interactions correspondingly.
One can see that modification in Eq. (2.16) changes the standard number rule Equation (2.16) allows to account for changes in the number of valence quarks. However, it cannot be used on its own since a subtraction of a valence quark also removes an average momentum fraction carried by q f vn (x, Q) 19) and this will cause the algorithm to violate the momentum rule Eq. (2.14) (which would be preserved if we only rescaled the momentum fraction in Eq. (2.16)). We will describe how this problem is resolved shortly; however, let us first describe how the 1 → 2 splitting effects described earlier are accounted for in the PYTHIA code. Let us say that in a particular interaction we find a sea quark q s in the proton. This must have been produced in conjunction with a companion (anti-)quark, q c : Let us denote the momentum fraction of the parent gluon by x g , that of the sea quark by x s and that of the companion by x = x g − x s . Then, according to Ref. [61], we introduce an sPDF for the companion quark, which is given by 8 where P g→qq (z) = 1 2 z 2 + (1 − z) 2 and C is a normalization constant which is fixed by Note that Eq. (2.21) is defined assuming no momentum loss due to the previous hard interactions and hence has to be rescaled accordingly. Using Eq. (2.21) one can compute an average momentum fraction introduced by a companion quark to a beam remnant X n x f c0 (x s ) defined in analogy with Eq. (2.19) as where index f is used to distinguish between different flavours and q f cn is linked to q f c0 via the momentum rescaling: q f cn (x, x s ) ≡ q f c0 (x/X n , x s )/X n . The average momentum introduced by a companion quark together with the average momenta subtracted due to the removal of valence quarks X n x f v0 (Q) leads to the violation of Eq. (2.14). The way to circumvent this issue is to let the sea and gluon distributions fluctuate in such a way that both contributions are compensated and the momentum rule Eq. (2.14) is preserved together with the number rules where index j is used to distinguish between different companion quarks of the flavour f . The corresponding modifications of the sea and gluon PDFs are given by which allows one to rewrite the momentum rule as 3 Comparison between the GS09 and PYTHIA models of dPDFs

Asymmetric double parton distribution functions
Now let us study in detail how well dPDFs constructed out of the sPDFs being used in the MPI model of PYTHIA satisfy the GS sum rules. We shall note that one can use raw and modified sPDFs to construct asymmetric PYTHIA dPDFs as in Eq. (2.11), or one can construct symmetric dPDFs according to where D j 1 j 2 (x 1 , x 2 , Q) and D j 2 j 1 (x 2 , x 1 , Q) are given by Eq. (2.11). The asymmetric dPDFs have a strict ordering of the arguments since the second sPDF f m←j 1 ,x 1 j 2 (x 2 , Q) is always modified according to the first hard interaction involving a parton j 1 carrying a momentum fraction x 1 . Such ordering appears because in the MPI model of PYTHIA the MPIs are considered as semi-hard corrections to the first hard process. However, one should keep in mind that the dPDFs should be symmetric under simultaneous interchange of flavours, Bjorken-x'es and factorization scales, and hence the GS sum rules should be satisfied when either x fraction is integrated over for the equal scale dPDFs. Nonetheless, it is interesting to see to what extent the asymmetric dPDFs satisfy the sum rules when integrating over the momentum fraction of the second "modified" parton, as this is the scenario where, by design, the PYTHIA dPDFs should best satisfy the sum rules. Therefore, we start the discussion with the asymmetric dPDFs as in Eq. (2.11) and then consider effects of symmetrization. Since the sum rules apply to the case where the scales in the dPDF are equal, we set Q 1 = Q 2 and set this common scale to 91 GeV.
First of all, let us check how the ansatz given by Eq. (2.11) satisfies the momentum sum rule given by Eq. (2.5) (when integrating over the second parton only). By substituting Eq. (2.11) into Eq. (2.5) we get where we use the symbol ? = to denote that this relation should be true for the sum rule to be satisfied. Since j 2 in j 2 runs over all parton species one can write Eq. (3.3) as Looking at Eq. (2.14) with n = 2, X = 1 − x 1 , we see that this relation does indeed hold. That is, the asymmetric PYTHIA dPDFs satisfy the momentum sum rule, when integrating over the second parton, by design. Later in this section we will investigate how well Eq. (2.5) is satisfied by the asymmetric PYTHIA dPDFs in practice. However, before doing that, let us first consider the number sum rule given by Eq. (2.6). It is convenient to express Eq. (2.6) in terms of a response function R j 1 j 2 defined as . (3.5) The meaning of R j 1 j 2 (x 1 , x 2 , Q) becomes clear if one neglects momentum and number conservation and substitutes D j 1 j 2 (x 1 , x 2 , Q) = f r j 1 (x 1 , Q)f r j 2 (x 2 , Q) in Eq. (3.5). The function R j 1 j 2 (x 1 , x 2 , Q) then turns simply into a raw valence quark sPDF multiplied by a corresponding Bjorken-x: (3.6) Using R j 1 j 2 (x 1 , x 2 , Q) one can write Eq. (2.6) as Since in the PYTHIA model where one can consider the difference f m←j 1 ,x 1 as the response of the sPDF of a valence quark of a flavour j 2 to the first hard interaction involving a parton of the flavour f 1 probed at the Bjorken-x equal to x 1 .
In order to be able to perform numerical integration in Eq. (3.9) one has to be able to access sPDFs modified by PYTHIA to simulate a second hard interaction. Normally, it is not possible since the corresponding part of the MPI model is hidden from a standard PYTHIA user. However, it can be achieved by instantiation of objects which are members of the BeamParticle class as explained in Appendix A. Each call of the routine second_xPDF from Appendix A returns the value of the modified sPDF being used for the second interaction which allows us to perform numerical integration in Eq. (3.9). One important comment has to be made: the aforementioned method to get the value of modified sPDF is essentially probabilistic which implies that the modified sPDF being plotted as functions of x 2 will have strong oscillations which make numerical integration in Eq. (3.9) unstable. Such oscillations appear due to the different treatment of valence and sea sPDFs in the DPS (MPI) machinery of PYTHIA. Therefore, depending on the value of x 1 PYTHIA decides in a probabilistic manner if f m←j 1 ,x 1 j 2 (x 2 , Q) should be evaluated assuming f r j 1 (x 1 , Q) to be either valence or sea sPDF. The value of f m←j 1 ,x 1 j 2 (x 2 , Q) will differ for each case which leads to the oscillations of f m←j 1 ,x 1 j 2 (x 2 , Q) as soon as x 1 approaches the region where the valence and sea sPDFs overlap. However, the problem can be solved by considering the values of f m←j 1 ,x 1 j 2 (x 2 , Q) averaged over a sufficiently high number of calls of second_xPDF routine at given values of x 1 and x 2 . Therefore, in the following, we always average f m←j 1 ,x 1 j 2 (x 2 , Q) before performing numerical integration (we use 10 6 and 10 4 averaging calls for the checks of the number and momentum rules correspondingly).
Apart from the dPDFs constructed out of the PYTHIA sPDFs we also consider the GS09 dPDFs. In order to check how well GS09 dPDFs satisfy the number sum rule we use the publicly available GS09 grids with 600 grid points [103]. Since the GS09 dPDFs are based upon the MSTW2008 LO PDF set [104] we use it in all our computations involving sPDFs.
Additionally, as a baseline, we use so-called "naive" model of dPDFs where one replaces a dPDF by a product of two sPDFs and a θ-function to enforce the kinematic constraint x 1 + x 2 ≤ 1: This approach neglects correlations in x-space and violates the GS sum rules. Moreover, as was shown in Refs. [88][89][90], this ansatz does not satisfy the dDGLAP evolution equations. Now let us check how well PYTHIA and GS09 models satisfy the number rule given by Eq. (2.6) taking as an example the R uu (x 1 , x 2 , Q) and R dd (x 1 , x 2 , Q) response functions. In Fig. 2 we compare response functions evaluated with the PYTHIA and GS09 dPDFs against raw valence u-and d-quark sPDFs. For our checks we consider six different fixed values of x 1 , namely 10 −6 , 10 −3 , 10 −1 , 0.2, 0.4 and 0.8. One can see how the 1 → 2 splitting, simulated according to Eq. (2.21), leads to negative values of the response functions at small values of x 2 . For example, if in the first interaction we probe a u-quark at x 1 = 10 −6 it is most likely a sea quark which was created due to a perturbative splitting of a gluon into a uū pair. Therefore, the remaining companionū-quark has to be included into ā u-quark sPDF which means that f r Fig. 2 one can see how the minimum of R uu (x 1 , x 2 , Q) propagates with the value of x 1 (see x 1 = 10 −6 and x 1 = 10 −3 cases). The same considerations, obviously, apply to the g → dd splitting. By contrast, when x 1 ≥ 10 −1 it becomes more likely to pick a valence quark and, therefore, PYTHIA starts to use Eq. (2.16) more frequently than Eq. (2.21), which results in a decrease of the peak of the valence distributions. We also would like to note that in the case of the R uu response function we observe almost a factor of two reduction in the height of the peak of the modified valence quark distribution comparing to the raw valence distribution if x 1 = 0.1 and that in case of the R dd response function the peak of the modified valence quark distribution almost disappears completely (as one would naively expect).
in PYTHIA is equal to 10 −10 ). Comparing results of numerical integration obtained with PYTHIA and GS09 dPDFs for the x 1 = 10 −6 case shown in Tables 1 and 2 we see that by setting a lower integration limit to 10 −8 in the PYTHIA code one can avoid violation of the number sum rule due to the boundary effects. Now let us check numerically how well the GS09 and asymmetric PYTHIA dPDFs satisfy  Table 2: Same as in Table 1 but for the R dd (x 1 , x 2 , Q) response function. In the ideal situation when the GS sum rules are perfectly satisfied N d v1 = 0.
the momentum sum rule. In order to do that we write Eq. (2.5) and Eq. (3.3) as and check how well the LHS of Eq. (3.11) and Eq. (3.12) match unity when we insert the GS09 dPDFs and the modified PYTHIA sPDFs respectively. The results of the numerical checks are given in Table 3. We see that, as expected from Eq. (3.4), PYTHIA dPDFs obey the momentum rule at the permille accuracy level and that GS09 dPDFs satisfy the momentum rule at a few percent accuracy level. We also note that the naive dPDFs obey the momentum sum rule at the few permille accuracy level if x 1 is in the range 10 −6 − 10 −3 , but show a strong disagreement if x 1 ≥ 0.1. The reason for this becomes clear if one substitutes Eq. (3.10) in Eq. (2.5) which leads to which holds only in the x 1 → 0 limit.  Table 3: Check of the momentum sum rule for the PYTHIA, GS09 and naive models of dPDFs. In the ideal situation when the momentum sum rule is perfectly satisfied the number in each cell should be equal to unity.
Similarly to the case of the number sum rule, it is interesting to visualize the behaviour of the integrands in Eqs. (3.11), (3.12). In Fig. 3 we plot the response function for the momentum sum rule: as a function of x 2 for six fixed values of x 1 , and for the case in which j 1 = d. In Fig. 3 a) we compare results obtained with PYTHIA dPDFs (solid lines) and GS09 dPDFs (dotted lines) whereas in Fig. 3 b) we compare PYTHIA predictions against results obtained with naive dPDFs (dashed lines). The blue lines in Fig. 3 correspond to the sum over all flavours in Eq. (3.14) whereas red, magenta and orange lines stand for partial contributions from gluon, d-andd-quarks correspondingly. Finally, the green lines are given by contribution from all quarks except d-andd-flavours. By comparing results obtained with PYTHIA and GS09 dPDFs shown in Fig. 3 a) we see that both models give similar predictions for all partial contributions except in the x 1 = 10 −6 and x 1 = 0.8 cases where the GS09 curves lie consistently below the PYTHIA ones. It is also noteworthy that, for all x 1 values probed, the bulk of area under the curves lies above x 2 = 10 −6 (in contrast to the corresponding plots for the number sum rule in Fig. 2, where we find that for x 1 = 10 −6 an important part of the curve lies below x 2 = 10 −6 ). This means that the value of 0.931 at x 1 = 10 −6 in Table 3 cannot be ascribed to the fact that the GS09 grids are limited to x i > 10 −6 , but is instead linked to the GS09 curve being somewhat "too low" at larger x 2 values ∼ 10 −2 − 10 −1 for x 1 = 10 −6 . This deficiency in the extent to which the GS09 set satisfies the momentum sum rule at the lowest x 1 values only appears after evolution -when one constructs Table 3 at Q = 1 GeV, one finds a value of 0.997 at x 1 = 10 −6 , and the sum rule integrand curve for GS09 lies almost exactly on top of the PYTHIA one for this x 1 value. We were not able to isolate a simple reason why this deficiency appears after evolution -it is presumably due to deficiencies at larger x i values migrating downwards, but the region x 1 = 10 −6 , x 2 = 10 −2 − 10 −1 is fed by various x values and flavours during evolution, making it very difficult to identify what features of the initial conditions are driving this effect. Since the main goal of this paper is in any case to study dPDFs and tPDFs generated using the PYTHIA model, and the GS09 results are used as a baseline, we defer further investigation of this phenomenon to future studies. Now let us discuss the momentum response functions for the naive dPDFs, as shown in Fig. 3 b). We see that if x 1 is equal to 10 −6 or 10 −3 there is no visible difference between the naive and PYTHIA response functions (solid and dashed lines practically coincide) which is in agreement with the results of the numerical integration shown in Table 3. We also see that starting from x 1 = 10 −1 the naive model starts to consistently overestimate the sum rule integrand. One would expect such an overestimate from the fact that the naive momentum sum rule curves do not go to zero at x 2 = 1 − x 1 as they should, but their trajectories are such that they would go to zero at x 2 = 1, if there were not the θ-function cutting them off. This effect is hardly noticeable for very small x 1 10 −1 , but becomes more and more significant at large x 1 -see e.g. the x 1 = 0.8 curve of Fig. 3 b). This increasing overshoot results in the naive model momentum sum rule integral deviating more and more from 1 as x 1 increases, as we observed in Table 3.
Before moving on, let us briefly comment on the case of unequal scale dPDFs in PYTHIA. Using the machinery discussed in Appendix A, we can access asymmetric unequal-scale dPDFs, with either Q 1 < Q 2 or Q 1 > Q 2 . The sPDF of the "second" parton is always modified as in Section 2.2. For such dPDFs it is clear that, when integrating over the second parton, the momentum and number and sum rules will hold for the unequal scale PYTHIA dPDFs to a very similar extent as they held for the equal scale PYTHIA dPDFs just discussed 9 . It is worth remarking that it is not necessarily expected theoretically that the dPDFs should satisfy the GS09 sum rules for the case of Q 1 = Q 2 . If one follows the procedure discussed in Section 2.1 to obtain the unequal scale dPDFs from the equal scale dPDFs (and it is the case that those equal scale dPDFs satisfy the sum rules) then the unequal scale dPDFs will satisfy the sum rules when integrating over the parton with higher Q, but not when integrating over the parton with lower Q.

Symmetrized double parton distribution functions
In Section 3.1 we demonstrated that the asymmetric dPDFs constructed as in Eq. (2.11) satisfy the GS sum rules in the second argument at a few percent accuracy level in the wide range of Bjorken-x's from 10 −6 up to 0.8. We also found that both PYTHIA and GS09 dPDFs have similar behaviour of the response functions as shown in Fig. 2 and Fig. 3. However, as we have mentioned before, QCD requires dPDFs to be symmetric under simultaneous interchange of parton flavours, Bjorken-x'es and factorization scales. Obviously, the dPDFs defined by Eq. (2.11) do not satisfy this requirement; in this section we check how well the symmetrized dPDFs defined in as Eq. (3.1) satisfy the GS sum rules. We start with the momentum sum rule. By substituting Eq. (3.1) in Eq. (2.5) we get where, similarly to Eq. (3.14), we defined . (3.16) In Table 4 we check Eq. (3.15) for the symmetrized PYTHIA dPDFs for the case where j 1 = u. As a reference we also add the check of Eq. (3.3) for asymmetric PYTHIA dPDFs as well as for the GS09 and naive dPDFs. We see that if x 1 = 10 −6 − 10 −1 the symmetric, asymmetric, and naive dPDFs satisfy the GS momentum sum rule at the 2% accuracy level; the GS09 dPDFs satisfy the sum rules to a similar extent except at x 1 = 10 −6 , as discussed in the previous section. We also note that the difference between the considered models increases with the value of x 1 . For example, if x 1 = 0.8 the asymmetric PYTHIA dPDFs and symmetric GS09 dPDFs satisfy the momentum sum rule at the 1% and 13% accuracy level correspondingly, whereas the symmetrized PYTHIA dPDFs and naive dPDFs are demonstrating strong violations of momentum conservation.  The additional terms introduced by the symmetrization procedure also lead to violation of the number sum rule given by Eq. (2.6). Let us illustrate this statement by considering uū andūu flavour combinations as the most prominent examples. Similarly to Eq. (3.5) we define a response function for the symmetrized PYTHIA dPDFs as which can be written as where, omitting factorization scale labels, we defined The second term on the right hand side of Eq. (3. 19) is new and leads to the violation of the number sum rule by the symmetrized PYTHIA dPDFs. Let us consider first the uū flavour combination. The R * uū (x 1 , x 2 , Q) response functions are given in Fig. 4 a) and the corresponding checks of the number sum rule are provided in Table 5. The checks for thē uu flavour combinations are given in Fig. 4 b) and Table 6. We see that in the uū case the deviation of the number sum rule integral from the expected value of Nū v1 ≡ −N u v1 = −1 is about 10 − 25% depending on the value of x 1 . However, in theūu case the deviations become stronger. In particular, the sum rule integral shows a deviation from the expected value of 3 by around 16% for moderate ([10 −3 , 10 −1 ]) values of x 1 and a large deviation for the valence region (about 130% deviation if x 1 = 0.8).  Table 5: Numerical integration with respect to x 2 over R * uū (x 1 , x 2 , Q) response function at fixed x 1 . In the ideal situation when the GS sum rules are perfectly satisfied Nū v1 = −1.  Table 6: Numerical integration with respect to x 2 over Rū u (x 1 , x 2 , Q) response function at fixed x 1 . In the ideal situation when the GS sum rules are perfectly satisfied N u v1 = 3.
In order to find out the origin of the peak of R * uu at x 1 = 0.8 that is responsible for this large deviation, we study in more detail the terms on the right hand side of Eq. (3.18). Let us start with the uū flavour combination. The corresponding terms are given by , (3.20) .
10 −5 10 −4 10 −3 10 −2 10 −1 10 −5 10 −4 10 −3 10 −2 10 −1 10 −5 10 −4 10 −3 10 −2 10 −1 In Fig. 5 a) we plot the four different terms given in Eqs. (3.20) -(3.21) (red and blue curves). The difference between Eq. (3.20) and Eq. (3.21) (which gives us R * uū ) is given by the green dash-dotted line. We see that all four terms have rather similar behaviour. Now let us consider theūu case which yields . (3.23) In Fig. 5 b) we plot the four different terms given in Eqs. (3.22)-(3.23). Unlike the uū case we see that now different terms have a different behaviour. Namely, we see that the term proportional to f r u (x 2 )f m←u,x 2 u (x 1 )/f r u (x 1 ) (red dotted line) has the biggest contribution and is largely responsible for the large peak given by the green dash-dotted line. Let us study the aforementioned ratio in more details. In Fig. 6 a) we plot a decomposition of the f r u (x 2 )f m←ū,x 2 u (x 1 )/f r u (x 1 ) term and in Fig. 6 b) we plot its analogue for theūu flavour . By comparing the middle plots in Fig. 6 we see that f m←ū,x 2 u (x 1 ) rapidly drops down in the region between 10 −2 and 10 −1 whereas f m←u,x 2 u (x 1 ) becomes peaked. The reason for such a different behaviour is as follows. In the case of f m←ū,x 2 u (x 1 = 0.8) we are looking at how finding aū quark with momentum x 2 affects the distribution of u quarks with momentum x 1 = 0.8. Since there is a sizeable distribution of u quarks at large x (arising from the valence contribution), the boost to this distribution from the companion quark effects is negligible, and the only visible effect from theū is just the "squeezing" caused by momentum conservation. That is, according to Section 2.2, for the middle plot in Fig. 6 a) at a given factorization scale we have However, in case of the middle plot in Fig. 6 b) the situation changes. Theū-quark sPDF, having no valence component, is much smaller at large x than the u-quark sPDF 10 . In this case, the companion quark effects are noticeable, particularly when x 2 is close in magnitude with x 1 . This results in a peak in the middle plot in Fig. 6 b). Additionally, the comparison between upper plots in Fig. 6 shows us that the x 2 f r u (x 2 ) term stays roughly the same in the region of interest (changes from 0.75 to 0.5) whereas x 2 f r u (x 2 ) drops from 0.6 down to zero. Combining upper and middle plots in Fig. 6 together and dividing by 2 f r u (x 1 = 0.8) or 2 f r u (x 1 = 0.8) we get a smoothly decreasing function in the lower plot in Fig. 6 a) and a peaked function in the lower plot in Fig. 6 b) correspondingly. Therefore, we conclude that the large peak in the response functions for the symmetrized dPDFs shown in Fig. 4  b) is due to the companion quark contribution to theū sPDF f m←u,x 2 u (x 1 ). We also shall note that similar large peaks occur for other flavour combinations where a sea quark sPDF receives a large contribution from a g → qq splitting, e.g. for the ss flavour combination as shown in Fig. 7.
Thus, we see that the symmetric dPDFs constructed using the MPI model of PYTHIA do not satisfy the sum rules perfectly. Nevertheless, we shall note that the largest violation of the GS sum rules by dPDFs constructed as in Eq. (3.1) occur at very large values of Bjorken-x'es (if x ≥ 0.4) where the absolute values of dPDFs are strongly suppressed. For the rest of the considered Bjorken-x'es the symmetric dPDFs obey the GS sum rules at 10 − 25% accuracy level. It would be interesting to modify the PYTHIA algorithm to suppress the large peaks in Figs. 4 and 7 such that we obtain symmetric dPDFs which obey the GS sum rules at a better accuracy level (for example, by suppressing the probability to find a companion quark accompanying a quark, when both partons have large x fractions). However, such a study is beyond the scope of this paper.

Double Drell-Yan production
In Section 3.1 and Section 3.2 we studied how well asymmetric and symmetric PYTHIA dPDFs satisfy the GS sum rules, and compared these to the GS09 dPDFs. Notably, we saw that the response functions defined by Eqs. (3.5) and (3.14) match rather closely between the GS09 and PYTHIA dPDFs.
Given this close match, an interesting question is to what extent the PYTHIA and GS09 dPDFs are similar overall. Here we study this question via phenomenological MC simulations using both sets of dPDFs. We also include predictions obtained using the naive dPDFs defined in Eq. (3.10), as a benchmark. It is most interesting to study a DPS process where the effects of the GS sum rules can potentially play an important role -i.e. one which receives contributions associated with two quarks of the same (anti)flavour belonging to the   same hadron. Therefore, we concentrate on a four-lepton production through the double Drell-Yan (dDY) process [80,[105][106][107][108][109][110], defined as (qq → γ * /Z * → 2l) ⊗ (qq → γ * /Z * → 2l), where l labels either an electron or muon final state. Before presenting our results, let us note that the older versions of PYTHIA (version ≤ 8.235) were using asymmetric DPS luminosities to simulate DPS production, whereas starting from version 8.240 PYTHIA uses the symmetric ones 11 Eq. (3.26) lead to the same results, since both equivalent ways to produce the A B final state (j 1 j 3 → A ⊗ j 2 j 4 → B and j 2 j 4 → A ⊗ j 1 j 3 → B) will be used approximately equal number of times. To check this statement we simulated the DPS distributions presented in this section with PYTHIA 8.235 and PYTHIA 8.240, and found basically no difference between both programs. Therefore, throughout this section we use the symmetrized DPS luminosities Eq. (3.26) only and omit the "sym" labels to keep our notation simple 12 .
To obtain the predictions with the GS09 and naive dPDFs, we change the weight of each PYTHIA event by multiplying it by a corresponding ratio of DPS luminosities: 13 Before discussing our results we would also like to note that according to the PYTHIA model the total DPS cross section is given by where σ ND is a total non-diffractive cross section [61]. This differs by a factor of σ eff /σ ND from the usual "pocket formula" result: σ DPS = σ 2 1 /2σ eff . The transition from Eq. (3.30) to the pocket formula expression in PYTHIA, however, is used only in the MPI machinery and is not performed if the flag SecondHard is turned on. This implies that the DPS cross sections in our simulations are evaluated according to Eq. (3.30) with σ ND instead of σ eff . However, since we are interested only in the relative difference between distributions of DPS events produced with different models of dPDFs, the aforementioned difference in the value of the DPS cross section does not affect our analysis.
The differential distributions generated with the naive, PYTHIA and GS09 dPDFs are given in Fig. 8. We perform binning in terms of the maximal rapidity difference 12 Bear in mind that the usage of the symmetic DPS luminosities as in Eq. (3.26) does not correspond to working with symmetric dPDFs. The luminosities obtained using symmetric dPDFs are given by: x2, x4, x3). 13 To perform the reweighting procedure one has to either use the method we describe in Appendix A to compute L Pythia or to be able to access DPS luminosities used by PYTHIA for a given DPS event. The latter can be achieved, for example, by adding new members to the Event class to store the information about the sPDFs being used to simulate the second interaction. This additional information has to be added to the event record in the ProcessLevel::nextOne(Event& process) method of the ProcessLevel class after the successful generation of the second hard interaction. The information assigned in such a way is then accessible together with the other standard information stored in the event record and can be used to evaluate ω GS09 and ω Naive during the generation procedure. Also note that by default in pp collisions PYTHIA assign a weight equal to one to each event. We see that the ratio dσ Naive /dσ Pythia is fairly flat in ∆Y, and remains rather close to 1 for the whole range of ∆Y (except at the 3 TeV collision energy where we see some increase at small and large values of ∆Y), whereas dσ GS09 /dσ Pythia decreases as one approaches the largest values of ∆Y. We see also that the ratio dσ GS09 /dσ Pythia remains more or less constant at small and moderate ∆Y as √ S is varied, whilst there is a noticeable decrease in this ratio at large ∆Y as √ S decreases. Let us now try to explain all of these features, starting with the very large ∆Y region (and in particular the last bin with ∆Y ∈ [9.5, 10]). We expect this ∆Y bin to probe a very asymmetric configuration of x fractions in the dPDF, with the values of x increasing as √ S is lowered. This is confirmed in the plots of Fig. 9, which also shows the values of x probed at each value of √ S. Why does the ratio dσ GS09 /dσ Pythia decrease as we decrease √ S and probe larger x values, whilst dσ Naive /dσ Pythia stays fairly close to 1? To try to understand this, we use a simplified setup, in which we consider only one partonic contribution; the one in which a u and aū from one proton collide with aū and a u in the other. This sort of contribution should be dominant at large ∆Y, as we can have a right-moving valence quark in the first subprocess, and a left-moving valence quark in the second. We fix the x values in the dPDFs to correspond to the production of two heavy states of mass M = √ x 1 x 2 S = 13 GeV, with the first being produced with rapidity Y = log(x 1 /x 2 )/2 = 4.5, and the other being produced with rapidity Y = −4.5 (yielding ∆Y = 9 overall), and for simplicity we consider the luminosity rather than the cross section. The corresponding results are given in Fig. 10 where we plot ratios of DPS luminosities as a function of collision energy √ S. As we mentioned before, Fig. 8 is generated using the particular choice of factorization scale Q = 91 GeV. However, in Fig. 10 we choose to generate the ratios for a range of factorization scales in between 1 and 180 GeV, as this helps us explain the observed behaviour at Q = 91 GeV.
Let us first study the Q-dependence of the L Pythia /L Naive ratio which is driven by the modification of the sPDFs affecting the L Pythia luminosity. We see in the upper plot of Fig. 10 that at Q = 91 GeV the ratio L Pythia /L Naive changes if √ S ∈ [3, 9] TeV and becomes almost constant if √ S > 9 TeV, which is consistent with our results shown in Fig. 8 14 . We also see that the ratio L Pythia /L Naive demonstrates a considerable change as the factorization scale increases from 1 to 180 GeV (compare orange and green curves). The Q-dependence becomes milder as Q increases, and is very gradual for Q ≥ 45 GeV. To investigate this behaviour let us note that for a given configuration (x 1 = x 3 , x 2 = x 4 , j 1 = j 3 , j 2 = j 4 ) L Pythia /L Naive is given by . (3.31) The kinematics we consider implies that x 1 ∈ [10 −2 , 1] whereas x 2 is well between 10 −4 and 10 −6 (for instance, if √ S = 3 TeV x 1 ≈ 0.4 and x 2 ≈ 4.81 × 10 −5 ). We found that for the Bjorken-x'es corresponding to Fig. 10 the second factor on the right hand side of Eq. (3.31) is very close to unity whereas the first one is responsible for the observed differences. Since the first factor on the right hand side of Eq. (3.31) implies that in the first interaction we probe a valence u-quark and a seaū-quark in the second one, then, according to the PYTHIA model, we deal with two modification mechanisms: first we squeeze all sPDFs according to Eqs. (2.12), (2.16) and then rescale the sea sPDFs by multiplying them by the a-factor as in Eqs. (2.26), (2.27) to preserve overall momentum conservation. We found that the rescaling of theū sPDFs according to Eq. (2.26) increases the value of f m←u,x 1 u whereas the squeezing according to Eq. (2.12) tends to make it smaller. The effect of the squeezing depends on the rate at which theū PDF grows at small x -if theū PDF at small x is proportional to x −b , then taking only the squeezing effect into account we would have For the MSTW set we use, b > 1 and grows rapidly at small Q ∼ 1 − 2 GeV before reaching an approximately stable value of b ∼ 1.5 at larger Q values. Thus, the squeezing effect corresponds to a suppression of the ratio f m←u,x 1 u /f r u , which increases at low Q but stabilises at larger Q values. This suppression also increases as x 1 increases, or √ S decreases. The "a-rescaling" effect is initially quite large -we find a ∼ 1.2 − 1.3 at Q = 1 GeV, which explains the (perhaps unexpected!) fact that L Pythia /L Naive significantly exceeds 1 at this Q value. This enhancement rapidly reduces to ∼ 1.1 at Q ∼ 5 GeV, and thereafter changes rather slowly. At Q = 1 GeV, there is a noticeable decrease in a as √ S increases, but at higher Q values the behaviour of a is rather flat in √ S. Thus, we see that the reduction of L Pythia /L Naive as Q increases is caused both by the squeezing suppression getting stronger and the a-rescaling enhancement getting weaker. At large Q values such as Q = 91 GeV the downward trend in the ratio as √ S decreases is mainly driven by the squeezing effect getting stronger, as the a-rescaling enhancement is roughly constant with √ S. Now let us consider the ratio L GS09 /L Naive shown in the middle plot of Fig. 10 15 . At Q = 1 GeV this ratio is very close to unity, which is essentially by design -when the x value of one parton is very small (as we have here), the GS09 inputs are constructed to match closely with the naive dPDFs, as this will ensure the momentum sum rule is well-satisfied when integrating over the other parton (c.f. Eq. (3.13)). At Q = 1 GeV the differences between the GS09 and naive dPDFs predominantly manifest themselves in the region where both x values in the dPDF are large (see Ref. [26]). However, as Q starts to increase we observe some important deviations of the ratio L GS09 /L Naive from 1 in Fig. 10 caused by the dDGLAP evolution. We see that if √ S is approximately in the [3, 10] TeV interval the ratio L GS09 /L Naive is smaller than one, whereas at large √ S we observe the opposite behaviour. This observation is consistent with the picture of the dDGLAP evolution from Ref. [26]; at the starting value of the evolution scale Q = 1 GeV the differences between the GS09 and naive dPDFs are located in the region where both x i are large, however, as Q grows the differences propagate down to smaller values of x 1 , x 2 leading to a decrease of the GS09 D uū dPDFs relative to the naive ones. Taking into account the fact that in our analysis the decrease in √ S leads to an increase of both x 1 and x 2 (see the distributions in Fig. 9) we conclude that if √ S ∈ [3, 10] TeV the x i values enter the region significantly affected by the aforementioned propagation of the differences between two models down to the smaller √ S ∆Y ∈ [0, 0.5] ∆Y ∈ [9.5, 10] √ S = 3 TeV, q i q j 53.3% 19.2% √ S = 7 TeV, q i q j 51.3% 29.9% √ S = 13 TeV, q i q j 50.7% 38.5% √ S = 20 TeV, q i q j 50.1% 42.7% Table 7: Flavour composition of the dPDFs belonging to the same proton and contributing to the first and last ∆Y bins (q i q j + q iqj = 100%).
values of momentum fractions. However, as √ S increases both x 1 and x 2 become smaller which lead us to the region of small x'es where the 1 → 2 splitting contributions due to the last term on the right hand side of Eq. (2.3) become important. Since the 1 → 2 feeding term causes an increase of the parton densities we conclude that the change of the behaviour of the L GS09 /L Naive ratio at large values of √ S is due to the positive contribution from 1 → 2 splittings.
Having now studied and understood the behaviour of the L Pythia /L Naive and L GS09 /L Naive ratios, we can now understand the behaviour of L GS09 /L Pythia as this can simply be constructed by taking the quotient of the ratios already studied. In particular, we note that at very low Q's the ratio L GS09 /L Pythia lies significantly below 1 predominantly due to the a-rescaling enhancement in PYTHIA. At large Q, the trend in √ S is mainly driven by momentum/number sum rule suppression in the GS09 dPDFs, which has started to propagate downwards into the region where one x i is large and the other small. As one can see the behaviour of the L GS09 /L Pythia at Q = 91 GeV is consistent with the distributions shown in Fig. 8. Now let us focus on the small ∆Y region in Fig. 8, and initially restrict our attention to the behaviour of dσ Naive /dσ Pythia . We saw in Fig. 8 that this ratio stays within 5% of 1 for the whole range of ∆Y at 7, 13 and 20 TeV collision energies. However, if √ S = 3 TeV we observe about 10% deviation of the aforementioned ratio from 1 at small ∆Y ∈ [0, 0.5].
The ∆Y ∈ [0, 0.5] bin should probe dPDFs where the two x i values are of the same magnitude -this is confirmed in Fig. 11 where we plot the distribution of Bjorken-x'es for the PYTHIA events contributing to this bin. We also see that the distribution of events in this common x value is approximately flat with the bulk of events having x ∈ [10 −5 , 10 −1 ]. Furthermore, unlike the ∆Y ∈ [9.5, 10] bin which predominantly probes dPDFs with flavour structure q iqj (i, j = {u, d, s, c, b}), the events in the ∆Y ∈ [0, 0.5] bin arise from an approximately even mix of events where a q i q j in one proton collide with aq iqj in the other, and events where we have a q iqj vs.q i q j collision -see Table 7. Finally, we found that amongst the five (anti-)quark flavours we consider, the combinations ofū,d, u and d flavours occurs most frequently.
In order to get a better understanding of dσ Naive /dσ Pythia in the region ∆Y ∈ [0, 0.5], we plot the ratio L Naive /L Pythia where we set all Bjorken-x'es to the same value in Fig. 12. In this figure we consider combinations of u and d flavours having different behaviour, namely L uuūū , L uūūu , L udūd and L udūd (here we skip dddd and dddd flavour combinations since they have behaviour similar to the analogous combinations of u andū quarks). We see that for all considered flavour combinations, except for the uūūu case, the ratio is bigger than unity for the whole range of x. The decrease of the ratio L Naive /L Pythia for the uūūu flavour combination below 1 is due to the companion quark contribution from g → uū splitting in the PYTHIA case. We also see that the ratio corresponding to the uuūū flavour combination increases faster compared to the udūd and udūd cases. This behaviour can be explained by the fact that for the uu dPDF we invoke not only momentum squeezing mechanism in PYTHIA, as in Eq. (2.10), but also have a suppression effect due to the change in remaining number of valence quarks, as in Eq. (2.16). It is also interesting to note that the effects due to the momentum and number sum rule constraints can penetrate down to rather small Bjorken-x values ∼ 10 −2 − 10 −1 . Since at √ S = 3 TeV collision energy we probe Bjorken-x'es at larger values than at 7 -20 TeV (see blue distributions in Fig. 11 a)), we conclude that the increase in the ratio of naive to PYTHIA distributions at small values of ∆Y at 3 TeV collision energy, as shown in Fig. 8, is due to the number and momentum sum rule effects which manifest themselves at the larger values of Bjorken-x'es accessible at the 3 TeV collision energy. x Figure 12: Ratio of the DPS luminosities L Naive /L Pythia as a function of x where Let us finally turn to considering the ratio dσ GS09 /dσ Pythia at small to moderate values of ∆Y ∈ [0, 5]. We see in Fig. 8 that this ratio is approximately constant for ∆Y ∈ [0, 5], and is ∼ 1.1 across all values of √ S. This enhancement of GS09 over PYTHIA is linked to 1 → 2 splitting contributions to the GS09 dPDFs during dDGLAP evolution, see the the last term on the right hand side of Eq. (2.3). To illustrate this, in Fig. 13 we plot the ratio of the ∆Y distributions obtained by solving dDGLAP evolution equation with and without the 1 → 2 splitting contribution. As one can see the evolution effects due to the 1 → 2 splittings increase the ∆Y distributions reasonably uniformly in the whole range of binning but somewhat more at the small values of ∆Y than at the large values of ∆Y. One notes that the enhancement factor when including the 1 → 2 splitting effects is ∼ 1.30 − 1.35 at small ∆Y, which is rather larger than the factor ∼ 1.1 enhancement of GS09 over PYTHIA that we see in Fig. 8. This implies that there is some suppression in GS09 from momentum/number sum rule suppression effects migrating to lower x values during evolution, which partially compensates the 1 → 2 splitting enhancement.
To finish this section, we would like to comment that despite the similarities between the response functions computed with PYTHIA and GS09 dPDFs, as shown in Figs. 2 and 3, the ∆Y distributions generated by these models of dPDFs differ from one another, especially at small and large values of rapidity separation. This implies some important differences (qq → γ * /Z * → 2l) ⊗ (qq → γ * /Z * → 2l) Figure 13: Comparison between ∆Y = max|y i − y j | DPS distributions generated with the GS09 dPDFs produced by solving dDGLAP evolution equation with and without 1 → 2 splitting. The collision energy √ S is equal to 3, 7, 13 and 20 TeV.
between the dPDF sets themselves. We also note that the ∆Y distributions generated with naive and PYTHIA dPDFs also show some important differences not only at large values of rapidity separation (something one would expect, since large values of ∆Y allow to probe dPDFs with one large x value, where momentum sum rule effects play a significant role) but also at small values of ∆Y at the 3 and 7 TeV collision energies.

Sum rules for triple PDFs
In the previous sections we studied how well the asymmetric and symmetric dPDFs constructed with PYTHIA code obey the GS sum rules. Now we would like to extend the analysis performed in Sections 3.1-3.2 to the case of triple PDFs. However, since the sum rules for tPDFs do not exist in the literature, we first need to establish what these sum rules actually are; in Section 4.1 and Section 4.2 we derive the GS sum rules for the bare and renormalized tPDFs correspondingly. Then, in Section 4.3, we investigate how well the PYTHIA tPDFs satisfy these sum rules.

Sum rule for bare tPDFs
We start by writing down the GS sum rules generalized to the case of bare tPDFs. The momentum sum rule is given by where T B j 1 j 2 j 3 and D B j 1 j 2 (x 1 , x 2 ) are the bare triple and double parton distribution functions respectively. To prove this relation we apply the same method as was used to show that GS sum rules hold for bare dPDFs in Ref. [93]. First of all let us write down the light-cone representation [112][113][114][115][116] of dPDFs and tPDFs: By substituting Eq. (4.3) in Eq. (4.1) we get The last sum in Eq. (4.6) can be written as where we used the fact that N m z m = 1 and that z k = x 1 and z l = x 2 . Therefore, Eq. (4.6) turns into Now let us concentrate on the generalization of the number rule Eq. (2.6) for the case of bare tPDFs. First of all let us note that and, therefore, (4.10) The two sums in the parentheses in Eq. (4.10) can be written as The difference between N j 3 |{β i } and Nj 3 |{β i } is the number of valence quarks, which is independent of {β i }: (4.12) Inserting Eqs. (4.11) and (4.12) into Eq. (4.10), we obtain which is Eq. (2.6) generalized for the case of bare tPDFs.

Sum rules for renormalized tPDFs
We have shown that Eq. (4.1), (4.13), hold for the bare distributions. Now let us sketch the all order proof for the renormalized tPDFs following the method similar to the one used in Ref. [68]. Using notation of Refs. [68,82], the bare sPDF, dPDF and tPDF are defined as , (4.14) where we use the light-cone representation for four-vectors: . The quantity n j is 0 if j is a quark or antiquark, and 1 if j is a gluon. The twist-two operators O in Eqs. (4.14)-(4.16) are bilinear in partonic fields and are defined as Momentum sum rule. To prove the momentum sum rule for the renormalized tPDFs, let us consider: For the momentum sum rule to hold, ∆ j 1 j 2 mtm must be zero. The first term on the right hand side of Eq. (4.29), after inserting the expression for T j 1 j 2 j 3 from Eq. (4.20), becomes: (4.30) We now simplify each of the terms in this expression. For the first term in Eq. (4.30), we have: (4.31) Going from the first line to the second line, we have used Eq. (64) from Ref. [68], where A and B are arbitrary single argument functions (specifically we have used the latter relation). Then, going from the second to the third line, we have used the momentum sum rule for the bare tPDFs, plus the momentum sum rule for the Z ij : j 3 XZ j 3 j 3 = 1 (see Eq. (97) in Ref. [68]). For the second term in Eq. (4.30), we have: (4.33) Once again we use Eq. (64) from Ref. [68] to go from the first to the second line, and then we use the momentum sum rule for the Z ij along with the momentum sum rule for the bare dPDFs (as proven in Ref. [68]). In one of the terms denoted by "(permutations)", we have: (4.34) Here we use Eq. (66) from Ref. [68] to go from the first to the second line, with D an arbitrary function of two arguments and A a function of one. Then to get to the third line, we use the momentum sum rule for the Z ij,k factors as given in Eq. (103) of Ref. [68]. The other term in the "(permutations)" has the same structure but with 1 ↔ 2.
Let us now consider the second term in Eq. (4.29). Expressing D j 1 j 2 in terms of bare quantities, using Eq. (70) from Ref. [68], we have: (4.36) Now many of the terms on the right hand side of the expression for ∆ mtm cancel. The terms in the , , boxes cancel against each other directly, whilst the terms in the boxes cancel against the other "(permutations)" term that we have not written out explicitly. All the terms containing D B on the right hand side then cancel out, and we are left with: (4.37) To reach the last line we have used: and (4.39) The quantity ∆ j 1 j 2 mtm must be finite, as it is expressed in terms of renormalized quantities. However, in Eq. (4.37) we see that it can be expressed in terms of a pure pole factor convolved with a finite quantity (f k ). The only way this can be true is if both left and right hand side are zero -thus ∆ j 1 j 2 mtm = 0 and the momentum sum rule for tPDFs holds at the renormalized level.
Given that ∆ j 1 j 2 mtm = 0, we obtain from Eq. (4.37) a sum rule for the Z ijk,l factors: Number sum rule. The proof of the number sum rule for renormalized tPDFs proceeds in an analogous fashion. Consider: If this is zero, the number sum rule holds. We consider 3 T j 1 j 2 j 3v , and express it in terms of bare quantities using Eq. (4.20): (4.42) Using Eq. (64) from Ref. [68], the number sum rule for the Z ij (see Eq. (90) from Ref. [68]), and the number sum rule for the bare tPDFs, we may express the first term of Eq. (4.42) as: (4.43) For the second term, we use Eq. (64) from Ref. [68], the number sum rule for the Z ij , and the number sum rule for the bare dPDFs, to obtain: (4. 44) In one of the terms in the "(permutations)", we can use Eq. (66) from Ref. [68] and the number sum rule for the Z ij,k (as given in Eq. (94) of Ref. [68]): (4. 45) The other term in the "(permutations)" has the same structure, but with 1 ↔ 2.
We express the second term in Eq. (4.41) in terms of bare quantities using Eq. (70) from Ref. [68]: (4.46) As for the momentum sum rule, a number of cancellations now occur. The terms in the , , , , boxes cancel against each other directly, whilst the terms in the boxes cancel against the other "(permutations)" term that we have not written out explicitly. Once again, the terms with D B cancel out, and we obtain (using Eq. (4.38)): According to a similar argument as was used for ∆ mtm , we find that ∆ num = 0 -that is, the number sum rule holds for the renormalized tPDFs. We also obtain a number sum rule for the Z ijk,l factors: (4. 48) 4.3 Check of the GS sum rules for tPDFs constructed with PYTHIA code In Sections 4.1 -4.2 we proved the GS sum rules generalized for the case of tPDFs. Now let us demonstrate how one can build tPDFs which approximately obey Eq. (4.1) and Eq. (4.13). As in Section 3.1 we begin the discussion with the case of the asymmetric tPDFs which are given by As for the dPDFs, we only check the sum rules when integrating over the final parton (parton 3 here) -this is the case in which we expect the sum rules to be satisfied the best. To construct tPDFs we use the PYTHIA code as it is explained in Appendix A. In this section we always set all scales in the tPDF to be equal (this is the case for which the tPDF sum rules should hold), and set the value of this common scale to 91 GeV. First of all let us check the momentum sum rule given by Eq. (4.1). Inserting Eq. (4.49) into Eq. (4.1), we get: (4.50) The results of the numerical check of Eq. (4.50) are given in Table 8. The problem of oscillations for the PYTHIA dPDFs mentioned in Section 3 becomes stronger for the case of tPDFs. Therefore, we increase a number of averaging calls by one order of magnitude comparing to the case of dPDFs. For our check we choose j 1 = j 2 = u, and vary x 1 between 10 −6 and 0.8 intervals while keeping x 2 equal to 10 −4 . As a baseline for our computation we use naive tPDFs defined in analogy with Eq. (3.10) as Now let us check the number rule given by Eq. (4.13). In order to do that, similarly to Eq. (3.5), we define which can be seen as a response of the valence sPDF f j 3v (x 3 , Q) to the two interactions involving partons j 1 and j 2 with momentum fractions x 1 and x 2 scales respectively. Due to the large number of different flavour combinations we restrict ourselves to a few illustrative cases. Let us first look at the uuu and uūu combinations. The results of our numerical checks are given in Tables 9 -10 and the corresponding response functions are given by blue lines in Fig. 14. We see that the aforementioned 1 → 2 splitting mechanism now give rise to a non-trivial modification of the valence u-quark sPDF after the first two interactions. As a consequence, the modified valence u-quark sPDF can possess several minima and maxima. According to Eq. (4.13) the integrations over the uuu and uūu response functions should   Table 9: Integration over R uuu response function with respect to x 3 at fixed x 1 , x 2 . In the ideal situation when the GS sum rules are perfectly satisfied N u v2 = 0.  Table 10: Integration over R uūu response function with respect to x 3 at fixed x 1 , x 2 . In the ideal situation when the GS sum rules are perfectly satisfied N u v2 = 2.
(4. 54) We see that for most of the cases the results of the numerical integration given in Table 9 and Table 10 agree with Eq. (4.53) and Eq. (4.54) up to a few per mille. Very similar results are obtained for the other flavour combinations. For example, for the sss case the number sum rule reads as follows: The response functions for this case are plotted in Fig. 15, and numerically we find that the integral on the left hand side of Eq. (4.55) does give zero up to deviations of order 0.001, for the same x values as were probed in the uuu and uūu cases. In Fig. 15 we see once again that the companion quarks produced via the g → qq splitting introduce additional strange and anti-strange quarks to the beam remnant, leading to peaks and troughs appearing in the response function.  Now let us move to discuss symmetrized tPDFs -we recall that the tPDFs should be symmetric when the flavour, x fraction and scale arguments for parton i are exchanged with those of parton j (where i, j ∈ {1, 3}), and that the (equal-scale) tPDFs should satisfy the sum rules when integrating over any of the partons, not just the final one. Similarly to Eq. (3.1) we define symmetrized PYTHIA tPDFs as where the sum on the right hand side of Eq. (4.56) runs over all partonic permutations. For the sake of simplicity we omit the factorization scale label in what follows (recall we take Q = 91 GeV in the numerics). The momentum rule Eq. (4.1), therefore, becomes (4.57) The numerical check of Eq. (4.57) for the j 1 = j 2 = u case is given in Table 11. We see that when both x 1 and x 2 are small the sum rule is pretty well satisfied (with deviations of  Table 11: Test of the momentum sum rule for the tPDFs. In the ideal situation when the momentum sum rule is perfectly satisfied the number in each cell should be equal to unity.  Table 12: Same as in Table 9 but for the symmetrized R * uuu response function. In the ideal situation when the GS sum rules are perfectly satisfied N u v2 = 0. the order of only 4%), but the violations become larger when one x fraction approaches 1, reaching the 70% level when x 1 = 0.8, x 2 = 10 −4 .
Similarly to the case of the symmetrized dPDFs we define the symmetrized response functions as , (4.58) where D sym j 1 j 2 are symmetrized PYTHIA dPDFs defined by Eq. (3.1). The numerical checks for R * uuu and R * uūu are given in Tables 12 and 13 and in Fig.14. We see that for R * uuu and R * uūu , the number sum rules are preserved reasonably well -at most of the x 1 , x 2 points tested, we have deviations of order 0.2 − 0.3. At certain points we see larger deviations up to ∼ 0.5 − 0.6 -for the uuu case, this occurs when one of the x values gets close to 1, but interestingly, for the uūu case, this occurs when the x fractions are at the smallest values tested (x 1 = 10 −6 , x 2 = 10 −4 ).
Whereas the symmetrized tPDFs in Fig. 14 show relatively small violations of the GS sum rules the problem of the large peaks discussed in Section 3 also shows for certain tPDFs, for instance for the sss combination in Fig. 15. Here we see a large violation of the number sum rule when x 1 is very large, x 1 = 0.8, similar to what we saw for the dPDF case in Fig. 7 Table 13: Same as in Table 10 but for symmetrized R * uūu response function. In the ideal situation when the GS sum rules are perfectly satisfied N u v2 = 2.

Summary and conclusions
We have demonstrated how one can use the MPI model of the PYTHIA event generator to construct asymmetric and symmetric dPDFs which approximately obey the Gaunt-Stirling (GS) number and momentum sum rules. The asymmetric dPDFs constructed in this way obey the GS sum rules at a few percent accuracy level in one of the arguments. We found that the symmetrized dPDFs obey the GS sum rules at about 20% accuracy level for x fractions smaller than 0.4. This is not perfect, but we remind the reader that the task of constructing dPDFs that even approximately satisfy the sum rules in this way is highly nontrivial -for example, in the GS09 paper [26], the sum rules are only satisfied "to better than 25% precision (for x < 0.8)". The concept of a "response function" was introduced for both the momentum and number sum rules, a function of x 1 and x 2 constructed out of dPDFs and sPDFs that, when integrated over x 2 , should give a constant value according the sum rules. The precise definitions of these response functions are given in Eqs. (3.5) and (3.14). An interesting observation we made is that even though only the "areas" under the response function curves in the x 2 direction are fixed by the sum rules, the overall response functions computed with PYTHIA and GS09 dPDFs show very similar behaviour.
We also compared the GS09 and PYTHIA dPDFs in terms of their predictions for the double Drell-Yan cross section, focussing in particular on the distribution in the maximal rapidity separation between the leptons (∆Y). We observed differences in shape between the GS09 and PYTHIA ∆Y distributions, with discrepancies up to the 10% level. This indicates that although the response function curves are rather similar between the two models, there are notable differences between the dPDFs themselves.
Additionally, we generalized the GS sum rules to the case of tPDFs. We demonstrated that these hold at the bare level using the light-cone wave function approach, and showed that they also hold at the renormalized level following the approach of Ref. [68]. We demonstrated how one can use the PYTHIA code to construct the asymmetric and symmetric tPDFs which approximately obey the corresponding version of the GS sum rules. We argue that in the absence of available sets of tPDFs one can use PYTHIA tPDFs as a good first approximation for phenomenological studies of the TPS phenomena.
Finally, we shall note that both GS09 and PYTHIA approaches do not take into account transverse dependence (y-dependence) between partons, and that the cross section formulae for DPS or TPS actually involve the quantities that include this transverse dependence, the double parton distributions (DPDs) or triple parton distributions (TPDs). Theoretical studies of DPS phenomena [37,50,86,[98][99][100]117] have demonstrated that proper account of the y-dependence results in contributions with 1 → 2 partonic splittings being numerically more important than in approaches neglecting transverse dependence of dPDFs, e.g. GS09.
Even though dPDFs and tPDFs do not appear directly in the rigorous treatment of the DPS and TPS cross section formulae, they are nonetheless linked to the DPDs and TPDs which do -the former can be obtained from the latter via an appropriate integral over transverse separation parameter y. One could, presumably, extend the PYTHIA approach to construct DPDs and TPDs approximately satisfying sum rules, following a similar approach to Ref. [102]. The DPDs and TPDs constructed in a such way would be a useful low-scale input for new tools (e.g. dShower [118,119]) that follow the theoretical rigorous approach to multiple scattering.
it takes the reference &beam to the object of the BeamParticle type, the Particle Data Group ID numbers, Bjorken-x'es x 1 , x 2 and factorization scales Q 1 , Q 2 . The reference to the beam object can be obtained by setting BeamParticle &beam = pythia.beamA; in the analysis code after the initialization of the PYTHIA object. The dPDF, therefore, is given by Eq. (2.11). return r e s ; } Listing 1: An example of a C++ function to access sPDFs being used for the second hard interaction.
On can easily extend this approach to get sPDFs after n subsequent MPI events. For example, if we need to get the sPDF being used to simulate the third interaction we need to append the lines given in Listing 2 to the function defined in Listing 1. The tPDFs are then given by Eq. (4.49).

B PYTHIA setup we use.
Here we provide the setup we have used for our simulations of the double Drell-Yan production.    [16] D0 collaboration, Evidence for simultaneous production of J/ψ and Υ mesons, Phys.