Double hard scattering without double counting

Double parton scattering in proton-proton collisions includes kinematic regions in which two partons inside a proton originate from the perturbative splitting of a single parton. This leads to a double counting problem between single and double hard scattering. We present a solution to this problem, which allows for the definition of double parton distributions as operator matrix elements in a proton, and which can be used at higher orders in perturbation theory. We show how the evaluation of double hard scattering in this framework can provide a rough estimate for the size of the higher-order contributions to single hard scattering that are affected by double counting. In a numeric study, we identify situations in which these higher-order contributions must be explicitly calculated and included if one wants to attain an accuracy at which double hard scattering becomes relevant, and other situations where such contributions may be neglected.


JHEP06(2017)083
1 Introduction The precise description of high-energy proton-proton collisions in QCD is imperative for maximising the physics potential of the LHC and of possible future hadron colliders. An important issue in this context is to understand the mechanism of double parton scattering (DPS), in which two pairs of partons undergo a hard scattering in one and the same proton-proton collision. In many situations DPS is suppressed compared with single parton scattering (SPS), but this suppression generically becomes weaker with increasing collision energy. For specific kinematics or specific final states, DPS can become comparable to or even larger than SPS. An overview of recent experimental and theoretical activities in this area can for instance be found in [1,2]. Consider a DPS process pp → Y 1 + Y 2 + X, where Y 1 and Y 2 are observed particles or groups of particles produced in two separate hard scattering processes, whilst X denotes all unobserved particles in the final state. A cross section formula has been put forward long ago in the framework of collinear factorisation, where the transverse momenta q 1 and q 2 of Y 1 and Y 2 are integrated over [3]. A corresponding expression has been given in [4,5] for transverse momentum dependent (TMD) factorisation, where q 1 and q 2 are small compared with the hard scales in Y 1 and Y 2 .
Inside a proton, the two partons that take part in the two hard scatters can originate from the perturbative splitting of one parton. The relevance of this splitting mechanism for the evolution equations of double parton distributions (DPDs) has been realised long ago [6,7] and studied more recently in [8][9][10]. However, it was only noted in [4,5] that the same mechanism dominates DPDs in the limit of small transverse distance between the two partons, and that the splitting contribution leads to infinities when inserted into the DPS cross section formula. These infinities are closely connected with double counting between DPS and SPS in particular Feynman graphs, a problem that had been pointed out earlier in the context of multi-jet production [11].
Different ways of dealing with this issue have been proposed in [12][13][14][15] and [16]. As discussed later in this paper, we find that these proposals have shortcomings either of theoretical or of practical nature. In the present work we present an alternative scheme for computing the cross section in a consistent way, including both DPS and SPS (as well as other contributions). Our scheme allows for a nonperturbative definition of DPDs in terms of operator matrix elements, and it is suitable for pushing the limit of theoretical accuracy to higher orders in α s .
Our paper is structured as follows. In section 2 we recall the theoretical framework for describing DPS and specify the theoretical problems mentioned above. The short-distance behaviour of DPDs is discussed in section 3, since it is essential for the scheme we propose in section 4. We first show how this scheme works at leading order (LO) in α s , before giving examples of its application at next-to-leading order (NLO) in section 5. Collinear DPDs evolve according to DGLAP equations, and in section 6 we discuss several consequences of this scale evolution. Our scheme is naturally formulated with DPDs that depend on the transverse distance y between the two partons, but we show in section 7 how one may instead use DPDs depending on the transverse momentum conjugate to y. This allows JHEP06(2017)083 us to compare our results with those of [12,13] and [14,15], which we do in section 8. Whilst the focus of the present work is theoretical rather than phenomenological, we give in section 9 some quantitative illustrations of our scheme, obtained with a relatively simple ansatz for the DPDs. We summarise our findings in section 10. Some Fourier integrals required in the main text are given in an appendix.

Setting the scene
In this section we recall theoretical issues originating from the perturbative splitting mechanism in double parton distributions, namely the appearance of ultraviolet divergences in the naive cross section for double parton scattering, the problem of double counting between DPS and single parton scattering, and the treatment of the so-called 2v1 (two versus one) contributions to DPS. 1 We also give some basic definitions and results. Throughout this work we use light-cone coordinates v ± = (v 0 + v 3 )/ √ 2 for any four-vector v. We write v = (v 1 , v 2 ) for the transverse components and v = √ v 2 for the length of the transverse vector.
Since the perturbative splitting mechanism in DPDs leads to issues in the ultraviolet (UV) region, renormalisation plays a crucial role in our context. We work in dimensional regularisation and extend the definitions of [5] to D = 4 − 2ǫ dimensions, changing phase space factors (2π) 4 into (2π) D . Bare two-parton distributions are then given by with twist-two operators for quark distributions, where q denotes the bare field. W (z) is a Wilson line from z to infinity along a prescribed path, which we do not recall here. The Dirac matrices Γ i select different quark polarisations. Analogous definitions hold for distributions involving gluons, with quark fields replaced by gluon field strength operators. For ease of writing, we omit colour labels on the operators and distributions throughout this work, bearing in mind that different colour couplings are possible for the four parton fields in (2.1).
In the process of deriving factorisation, one finds that the proton matrix element in (2.1) needs to be multiplied by a combination of so-called soft factors, which are vacuum expectation values of products of Wilson lines. More information for the case of single parton distributions can be found in chapter 13.7 of [18] and in the recent overview [19]. A brief account for DPDs is given in section 2.1 of [20], and a more detailed discussion will be provided in [21]. The product F bare × {soft factors} depends on a parameter ζ which regulates rapidity divergences. A scheme in which a soft factor appears explicitly in the cross section formula was presented in [22]. Since the treatment of soft gluons does not JHEP06(2017)083 affect parton splitting in any special way, we will not discuss it further in the present paper. Correspondingly, we will suppress the argument ζ in all DPDs.
As a final step one performs UV renormalisation, which we assume to be done in the MS scheme. The DPDs obtained from (2.1) are appropriate for TMD factorisation, with the transverse positions z i being Fourier conjugate to the transverse parton momenta that determine the final state kinematics. These distributions renormalise multiplicatively, with one renormalisation factor for each product of a quark field (or gluon field strength) with the Wilson line at the same transverse position and one factor for each pair of Wilson lines at the same transverse position in the soft factors. Denoting the product of these factors with Z, one obtains the final DPD as F = lim ǫ→0 Z × F bare × {soft factors} . It obeys a renormalisation group equation which is a straightforward generalisation of the one for TMDs (given e.g. in [19]).
The DPDs needed for collinear factorisation are obtained by setting z i = 0 in F bare and in the associated soft factors, before renormalisation. Setting z i = 0 introduces ultraviolet divergences in the operators O 1 and O 2 , and in the associated soft factors. The renormalised DPDs are then obtained as F = lim ǫ→0 Z 1 ⊗ Z 2 ⊗ F bare × {soft factors} , where Z i renormalises the operators associated with parton i and where the convolution products are in the momentum fractions x i . In the colour singlet channel, where both operators O i in (2.1) are colour singlets, the soft factors reduce to unity and one obtains the renormalised twist-two operators that appear in single parton densities.
Since the operators associated with partons 1 and 2 renormalise independently (both for the TMD and the collinear case) one may choose different renormalisation scales µ 1 and µ 2 in each of them. This is useful when the two hard subprocesses in double parton scattering have widely different hard scales. In particular, one can then approach the kinematics of the so-called underlying event, with a very hard scattering at scale µ 1 and additional jet production at a much lower scale µ 2 (of course µ 2 needs to remain in the perturbative region for our factorisation approach to be justified).
With different scales µ 1 , µ 2 in the collinear DPDs, we have a homogeneous evolution equation d d log µ 2 1 F a 1 a 2 (x 1 , x 2 , y; µ 1 , µ 2 ) in µ 1 and its analogue for µ 2 . For colour singlet DPDs, the kernels P a 1 b 1 on the r.h.s. are the usual DGLAP splitting functions for single parton densities. In colour non-singlet channels, both the DPDs and the splitting kernels have an additional dependence on the rapidity parameter ζ [21]. Following the notation of [5], the labels a i , b i in (2.3) denote both the species and the polarisation of the partons. The relevant labels are q, ∆q, δq for unpolarised, longitudinally polarised and transversely polarised quarks, likewiseq, ∆q, δq for antiquarks, and g, ∆g, δg for unpolarised, longitudinally polarised and linearly polarised gluons, respectively. Note that the polarisations of the two partons can be correlated with each other, even in an unpolarised proton.

JHEP06(2017)083
To simplify our presentation, we will consider the production pp → V 1 + V 2 + X of two electroweak gauge bosons V i = γ * , Z, W . Our results readily generalise to other processes for which DPS factorisation can be established; in the case of TMD factorisation this requires that the produced particles are colour singlets. We denote the four-momenta of the two bosons by q i , their squared invariant masses by Q 2 i = q µ i q iµ and their rapidities by Y i = 1 2 log(q + i /q − i ). We work in the proton-proton centre-of-mass frame, taking the proton with momentum p (p) to move in the positive (negative) 3 direction. Furthermore we define with s = (p +p) 2 . For the phase space of each gauge boson we have The "naive" cross section formulae (not taking into account the UV problems discussed below) read for TMD factorisation and for collinear factorisation. The combinatorial factor C is 1 if the observed final states of the hard scatters are different and 2 if they are identical. For simplicity we will consider the case C = 1 throughout this paper, unless mentioned otherwise. As explained in section 2.2.1 of [5], there are further contributions with DPDs that describe the interference of different parton species. They can be discussed in full analogy to the contributions given in (2.6) or (2.7), and we do not treat them explicitly in the present work for ease of notation. The hard scattering cross sectionsσ i in (2.6) are for the exclusive final state V i with transverse momentum q i , which must satisfy q i ≪ Q i . By contrast,σ i in (2.7) is integrated over q i and inclusive for final states V i + X. At leading order in α s it contains δ functions that enforce x ′ i = x i andx ′ i =x i . The subtractions for collinear and soft regions inσ i are different in the two factorisation frameworks, but in both cases they lead to a dependence on the factorisation scale µ i that cancels against the µ i dependence of the DPDs, up to powers in α s beyond the accuracy of the calculation. This happens separately for the two hard scatters (i = 1, 2) and by construction works exactly as in the case of single hard scattering.

JHEP06(2017)083
As was pointed out in [5], the framework discussed so far suffers from problems in the region of small transverse distances between the two partons in a DPD. The leading behaviour of the collinear distributions F (x i , y) at small y can be computed from the perturbative splitting of one parton into two and gives a behaviour like y −2 up to logarithmic corrections. When inserted in the factorisation formula (2.7) this gives a quadratically divergent integral at small y, which clearly signals an inappropriate treatment of the ultraviolet region. As we will review in section 3.1, the short-distance behaviour of the distributions F (x i , z i , y) is less singular but still leads to logarithmic divergences in the TMD factorisation formula (2.6).
Instead of using DPDs depending on transverse positions, one may Fourier transform them to transverse momentum space, integrating 1 (2π) 4−4ǫ d 2−2ǫ z 1 d 2−2ǫ z 2 e −i(z 1 k 1 +z 2 k 2 ) (2.8) for TMDs and d 2−2ǫ y e iy∆ (2.9) for TMDs and collinear distributions. For collinear distributions, this transformation must be made before subtracting UV divergences and setting ǫ → 0: with F (x i , y) ∼ y −2 in D = 4 dimensions, the Fourier integral (2.9) would be logarithmically divergent at y = 0. In D = 4 − 2ǫ dimensions this singularity turns into a pole in 1/ǫ, which requires an additional subtraction as we will review in section 7. Rather than being associated with the individual operators O 1 and O 2 , this subtraction is related to the singularity arising when the transverse distance y between the two operators goes to zero. It leads to an inhomogeneous evolution equation for the DPD F (x i , ∆) in momentum space, which has been extensively discussed in the literature [6][7][8][9][10] for the case ∆ = 0. Notice, however, that this extra µ dependence does not cancel in the cross section when (2.7) is rewritten in transverse momentum space. Moreover, the additional UV renormalisation of F (x i , ∆) does not remove all UV divergences at the cross section level. The singularity of F (x i , y) at small y translates into a behaviour F (x i , ∆) ∼ log(µ 2 /∆ 2 ) at large ∆ (see section 7), which gives a quadratic divergence for the ∆ integration in the DPS cross section. It is easy to identify the origin of the UV divergences just discussed. Both in the y and ∆ representations, one has integrated over the full range of the integration variable and thus left the region in which the approximations leading to the DPS cross section formulae are valid, namely the region where ∆ ≪ Q i or, equivalently, y ≫ 1/Q i (see section 2.1.2 in [5]). Outside this region, the DPS approximations are not only unjustified, but they give divergent integrals in the cross section.
This points to another problem, namely that of double counting contributions between SPS and DPS. To see this, let us analyse the graph in figure 1a. Since the transverse boson momenta q 1 and q 2 are approximately back to back (up to effects from the transverse momenta of the incoming gluons) it is convenient to introduce the combination For q ≪ Q i the graph gives a leading contribution to dσ/dq 2 if the transverse quark momenta in the loops are all of order Q i . This carries the quark lines far off shell, so that this contribution is naturally associated with SPS, with g + g → V 1 V 2 as the hard subprocess. A leading contribution is also obtained from the region where all transverse quark momenta are much smaller than Q i . This region is naturally described as DPS, with two disconnected hard scattering processes qq → V 1 and qq → V 2 and double parton distributions with perturbative g → qq splittings, indicated by the boxes in the graph. We denote this as a 1v1 (1 versus 1) contribution to DPS, emphasising its close relation to SPS. A double counting problem for this graph obviously arises if one takes the loop integrals in the SPS cross section over all transverse quark momenta (including the DPS region), and likewise if one integrates the DPDs cross section over the full range of transverse positions, which is equivalent to integrating over all transverse momenta in the quark loops. The problem persists if one integrates the cross section over q.
Let us now turn to the graph in figure 1b and consider the cross section integrated over q 1 and q 2 . In the region of large q ∼ Q i , the quark lines at the top of the graph have transverse momenta of order q and are far off shell. The proper description of this region is in terms of a hard scattering qq + g → V 1 V 2 , convoluted with a collinear twisttwo distribution (i.e. an ordinary PDF) at the top and a collinear twist-four distribution at the bottom of the graph. For brevity we refer to this as the "twist-four contribution" henceforth. In the region q ≪ Q i , the g → qq splittings are near collinear and the approximations for DPS are appropriate. We call this the 2v1 contribution to DPS, recalling that there is a qq + g subprocess in the graph. Both small and large q give leading contributions to the integrated cross section, and in a naive calculation adding up the twist-four term and the DPS term has again a double counting problem, as well as divergences in each contribution. The naive DPS cross section has a logarithmic divergence at small y, which is seen by inserting the 1/y 2 splitting behaviour of only one DPD in (2.7). In turn, the hard scattering cross section in the twist-four term contains a collinear divergence in the form of an integral behaving like dq 2 /q 2 at q → 0, as we will show in section 4.1.2.

JHEP06(2017)083
Clearly, one needs a consistent scheme for computing the overall cross section, without double counting and without divergences in individual terms. An intuitive approach for evaluating DPS is to separate the "perturbative splitting" part of a DPD from its "intrinsic" nonperturbative part. 2 This has been pursued independently by Blok et al. [12,13] and by Ryskin and Snigirev [14,15]. Taking the intrinsic DPD for each proton, one obtains the 2v2 part of DPS, which does not contain any perturbative splitting and is shown in figure 1c. The splitting part of the DPD is explicitly computed in terms of a single parton distribution function (PDF) and a perturbative kernel. This is multiplied with an intrinsic DPD to compute the 2v1 term. Finally, the product of two splitting DPDs is used to compute the 1v1 contribution in the approach of [14,15], where an ultraviolet cutoff must be imposed to regulate the quadratic divergence we mentioned earlier. By contrast, the authors of [12,13] advocate to omit this term and replace it entirely with the SPS contribution to the cross section.
We are, however, not able to give a field theoretic definition of the "intrinsic" or "nonperturbative" part of a DPD. The consideration of Feynman graphs in the preceding arguments is instructive, but a satisfactory definition should only appeal to perturbation theory in regions where it is applicable. We regard a nonperturbative definition of DPDs as indispensable for a systematic theory approach, for instance for deriving evolution equations and other general properties.
The setup we propose in this work defines DPDs as operator matrix elements as described above, containing both splitting and intrinsic contributions. UV divergences in the DPS cross section are avoided by introducing (smooth or hard) cutoffs in the integrations over transverse distances. The double counting problems are treated within the subtraction formalism used in standard factorisation theorems, described in detail in sections 10.1 and 10.7 of [18] and briefly recalled in section 4 here. The subtraction terms that avoid double counting also remove the above mentioned collinear divergence in the twist-four term. A distinction between "splitting" and "intrinsic" contributions to a DPD will be made in the limit of small transverse distances, where it can be formulated in terms of an operator product expansion (see section 3.3), and when making a model ansatz for DPDs at large distances, which is of course necessary for phenomenology.

Contributions to the cross section: power behaviour and logarithms
In preparation for later sections, we now recall some results for the power behaviour of different contributions to the cross section, referring to section 2.4 of [5] for a derivation. We also recall which logarithms appear in the lowest order 1v1 and 1v2 graphs. As already stated, we take the process pp → V 1 + V 2 + X as a concrete example, but the discussion readily generalises to other cases.
The differential cross section dσ/(d 2 q 1 d 2 q 2 ) in the region q i ≪ Q i can be computed using TMD factorisation. Here and in the following we write Q i to denote the generic size of Q 1 and Q 2 , and likewise for q i . The transverse momenta q i may be of nonperturbative JHEP06(2017)083 size Λ or much larger. In the latter case, further simplifications are possible by expressing transverse momentum dependent distributions in terms of collinear ones [21], but we shall not discuss this here. The leading power behaviour of the cross section is When q i goes to zero, it should be replaced by Λ. Three types of mechanisms contribute to the leading behaviour, namely DPS, SPS, and the interference between SPS and DPS. Corresponding graphs are shown in figures 1c, 2a and 2b, respectively. As discussed in the previous section, certain graphs contribute both to DPS and to SPS, depending on the kinematics of their internal lines. The 1v1 graph in figure 1 also has leading regions in which one of the loops is hard and the other is collinear. These regions contribute to the SPS/DPS interference, as shown in figure 2c. The double counting problem thus concerns both SPS, DPS and their interference. Note that the SPS graph in figure 2a contributes to the SPS/DPS interference but not to DPS.
Both the amplitude and its conjugate in the 1v1 graph contains a loop integral that behaves like d 2 k/k 2 in the region Λ, q i ≪ k ≪ Q i . When integrated over the full phase space, each loop thus builds up a so-called DPS logarithm, namely log(Q i /q i ) when q i > ∼ Λ and log(Q i /Λ) when q i < ∼ Λ. Whether these logarithms reside in SPS, DPS or their interference depends on how exactly one handles the double counting problem. We come back to this issue in section 4.2.
Let us now turn to the cross section integrated over q 1 and q 2 , which can be described using collinear factorisation. The leading power behaviour of the cross section, σ ∼ 1/Q 2 i , is given by the SPS mechanism alone. Several mechanisms contribute at the power suppressed level σ ∼ Λ 2 /Q 4 i . These are 1. DPS, which is suppressed because it can only populate the region q i ≪ Q i rather than the full phase space up to q i ∼ Q i , Figure 3. (a) A graph with a twist-four distribution for one proton and a twist-two distribution for the other. (b) A graph with twist-three distributions for both protons.
2. the interference between SPS and DPS, which is suppressed for the same reason, 3. hard scattering with a twist-four distribution for one proton and a twist-two distribution for the other. Example graphs are figure 3a, as well as figure 1b with the box removed.
4. hard scattering with twist-three distributions for both protons. An example graph is figure 3b.
The rationale for considering such contributions is that -whilst being power suppressed compared with SPS -they may be enhanced by higher parton luminosities at small momentum fractions x, or by coupling constants in the relevant hard scattering subprocesses. Let us emphasise that a complete calculation of the cross section at the level of Λ 2 /Q 2 i corrections would be a formidable task, and it is not even established whether factorisation (in particular the cancellation of Glauber gluons) holds at that level. Notice that in collinear factorisation, the SPS/DPS interference term involves collinear twist three distributions for both protons, because the SPS mechanism forces the two partons in the interfering DPS amplitude to be at same transverse position (see section 2.4.1 in [5]). In this sense, mechanism 2 in the above list may be regarded as a special case of mechanism 4, with a disconnected hard scattering in the amplitude or its conjugate (see figure 2b).
A full treatment of contributions with twist-three or twist-four distributions is beyond the scope of this paper. We remark however that twist-n operators contain n or less than n parton fields, and that different operators are related by the equations of motion. For a detailed discussion we refer to [31][32][33]. Twist-n operators with n parton fields were called "'quasipartonic" in [34] and involve only the "good" parton fields in the parlance of lightcone quantisation [35]. These are exactly the fields appearing in the definitions of multiparton distributions, so that graphs with a double counting issue between higher-twist hard scattering and DPS (or the SPS/DPS interference) involve only quasipartonic operators.
The matrix elements of quasipartonic twist-three operators in an unpolarised target satisfy the important selection rule that the helicities carried by the parton lines must balance. This excludes three-gluon operators since three helicities ±1 cannot add up to zero. For quark-antiquark-gluon operators it forces the quark and antiquark fields to have JHEP06(2017)083 opposite chirality, i.e. one only has the operator combinationqσ +j q, where the transverse index j is contracted with the polarisation index of the gluon. As for non-quasipartonic twist-three distributions in an unpolarised target, one finds that they are absent in the pure gluon sector [36], whereas the corresponding quark-antiquark distributions are again chiral-odd [37]. Since chiral-odd distributions cannot be generated by gluon ladder graphs, they lack the small x enhancement that is one of the motivations to keep higher twist contributions in the cross section. We will therefore not discuss them further in this work. Note that corresponding selection rules do not hold for TMD correlators, where an imbalance in the helicities of the parton fields can be compensated by orbital angular momentum.
Let us finally recall the appearance of DPS logarithms in collinear factorisation. The 2v1 graph (figure 1b) has a behaviour dσ/dq 2 ∼ 1/q 2 in the region Λ ≪ q ≪ Q i , which gives a log(Q i /Λ) when integrated over the full phase space. Depending on how the double counting between DPS and the twist-four mechanism is resolved, this logarithm can appear in different contributions to the cross section. We will discuss this in section 4.1.2.

Short-distance limit of DPDs
In this section, we analyse the behaviour of DPDs in the limit where the transverse distance between partons becomes small compared with the scale of nonperturbative interactions. In this region, the splitting of one parton into two becomes dominant. Generalising results in [5] we give expressions in D = 4 − 2ǫ dimensions, which are necessary in intermediate steps when constructing a factorisation formula for the cross section.

TMDs
A useful choice of position variables for describing the parton splitting mechanism is with Fourier conjugate momenta 3 The relation between DPDs in position and momentum space reads in D = 4 − 2ǫ dimensions. As seen in figure 4, one can identify y + (y − ) as the transverse distance between the two partons on the left (right) hand side of the final state cut in the DPD. Correspondingly, the transverse momentum difference between the partons on the left (right) hand side of the cut is k − (k + ). The splitting singularities of the DPDs thus occur at y ± → 0 or k ± → ∞.
JHEP06(2017)083 Figure 4. The perturbative splitting mechanism for a DPD, with momentum and position assignments. Here and in the following, the line for the final-state cut of the spectator partons is not shown for simplicity.
The perturbative splitting contribution F spl,pt to transverse-momentum dependent DPDs in momentum space has been calculated at leading order in section 5.2.2 of [5]. Generalising these results to D = 4 − 2ǫ dimensions, we have where j, j ′ are transverse Lorentz indices and f a 0 (x 1 + x 2 , K) is an unpolarised singleparton TMD. 4 The ellipsis denotes a term that involves a TMD for polarised partons in an unpolarised proton and depends on K but not on k ± . In position space we then get using the Fourier integral (A.2), where the term denoted by an ellipsis depends on Z but not on y ± . It is understood that for transverse quark or linear gluon polarisation, both F a 1 a 2 and the kernel T carry additional transverse Lorentz indices. f a 0 (x 1 + x 2 , Z) is the Fourier transform of f a 0 (x 1 +x 2 , K). The form (3.4) gives the leading behaviour of the DPD for large k ± ≫ Λ, and correspondingly (3.5) gives the leading behaviour for y ± ≪ 1/Λ. If one inserts these results into the cross section formula and sets D = 4, logarithmic divergences appear at y + = 0 and y − = 0. To make them explicit we transform variables to with q defined in (2.10). Performing the angular integration in

JHEP06(2017)083
where y stands for y + or y − , we see that the integral with J 0 is divergent at y = 0. Given the range of validity of (3.5) one should impose y ≪ 1/Λ in (3.7), although the integrals are finite for y → ∞ due to the oscillations of the Bessel functions. The perturbative splitting contribution to DPDs at higher order in α s involves graphs with additional partons radiated into the final state as shown in figure 5a, as well as virtual corrections. It is natural to expect that it will again be singular at y ± = 0. A calculation of this contribution is outside the scope of the present work, so that we will limit our analysis of TMD factorisation to perturbative splitting at LO.
To compute the DPD cross section, we must also consider the case where only one of the distances y + or y − is small, whereas the other one remains large. In this case, one has a perturbative splitting only on one side of the final-state cut, as illustrated in figure 5b. We will not discuss the detailed expression of the DPD in this regime, but give its general structure. Setting D = 4 for simplicity, we have were U α 0 →α 1 α 2 is a kernel for the splitting α 0 → α 1 α 2 in the amplitude (hence its complex conjugate appears in (3.8)). D α 1 α 2 |α 0 is the position space version of a transversemomentum dependent twist-three distribution, constructed from the hadronic matrix element where φ i is a "good" field for parton α i (cf. section 2.1). Distributions D α 0 |α 1 α 2 (x i , y − , Z) where α 0 belongs to the amplitude and α 1 , α 2 to the complex conjugate amplitude are defined in analogy. In the second step of (3.9) we have used translation invariance and shifted the parton fields to the same position as in the corresponding DPD (see figure 4). The labels α i denote the parton species; it is understood that in (3.8) and the following equations parton helicities are taken fixed on the l.h.s. and must be appropriately summed over on the r.h.s. Note the difference between this notation and the labels a i , which denote parton species and polarisation (none, longitudinal, transverse or linear) and thus refer to a pair of parton legs. The notation with a i is hence not suitable for distributions with three parton fields. If y + is small, then D α 1 α 2 |α 0 (x i , y + , Z) itself can be generated by perturbative splitting, as shown in figure 5c. We have Notice that a quark and antiquark produced by perturbative splitting have opposite helicities, so that the corresponding quark-antiquark operator φ 2 φ 1 in D qq |g must be chiral-even.

Collinear DPDs: splitting contribution
We now turn to collinear DPDs, i.e. to the case where z 1 = z 2 = 0. Let us first consider distributions for two unpolarised or two longitudinally polarised partons, so that the DPDs do not carry any transverse Lorentz indices. The lowest order splitting has been computed in [5]. For 4 − 2ǫ dimensions, one finds the general form F a 1 a 2 ,spl,pt (x 1 , x 2 , y; µ) The kernel for the splitting g → qq reads for instance with a factor f = 1 for the colour singlet and f = −1/ √ N 2 − 1 for the colour octet DPD. In terms of the kernel in (3.4) we have T jj ′ g→qq = δ jj ′ P g→qq . We recognise in P g→qq (u, 0) the usual DGLAP splitting kernel without the terms proportional to δ(1 − u).
Going beyond leading order, one can deduce the general form of the perturbative splitting contribution using dimensional analysis and boost invariance. For colour singlet distributions one finds 1 F a 1 a 2 ,spl,pt (x 1 , x 2 , y; µ) in D = 4 − 2ǫ dimensions. The convolution integral over v is familiar from factorisation formulae for hard scattering processes. Both f and V on the right-hand side are understood to include all necessary subtractions, so that they are finite at ǫ = 0. The splitting kernel V is a double series The µ (and thus on dimensional grounds the y) dependence of V follows from the fact that the mass parameter of dimensional regularisation appears in graphs only via µ 2ǫ α s (µ); terms with n > m in (3.15) are due to the subtractions of ultraviolet or collinear divergences. At lowest order, the hard splitting graphs are disconnected (with no partons across the final state cut), so that (3.14) this gives a form consistent with (3.12). Using that at order α n s the poles of highest order are 1/ǫ n−1 , we find in the physical limit ǫ = 0.
For colour nonsinglet DPDs one must regulate rapidity divergences, which complicates the preceding result. Taking e.g. Wilson lines along non-lightlike paths introduces additional vectors and changes the analysis of boost properties of the kernel. We will not pursue this issue here.
DPDs with transverse quark or linear gluon polarisation carry transverse Lorentz indices. Their perturbative splitting expressions thus have a tensor structure containing additional factors of y j /y compared with the formulae above. At leading order one readily finds from (3.5) and the appropriate splitting kernels that the factor 1/y 2−4ǫ in (3.12) is to be replaced with y j y j ′ /y 4−4ǫ times a tensor constructed only from Kronecker deltas.

Collinear DPDs: all contributions
Let us now study the small y behaviour of collinear DPDs in more general terms. We start by writing the relation between unrenormalised DPDs in position and momentum space as omitting all arguments other than y and ∆. The first term on the r.h.s. is a collinear twist-four distribution, independent of any transverse variable. For small y, the second term is dominated by large ∆, so that one can replace F (∆) by its approximation for large ∆, following the power counting analysis of section 5.2 in [5]. This leads us to write F y→0 (y) = F spl,pt (y) + F tw3,pt (y) + F int,pt (y) , (3.18) where the three terms on the r.h.s. will be described shortly. In D = 4 dimensions, they respectively go like y −2 , y −1 and y 0 , up to logarithmic corrections. Further terms from the perturbative expansion of F (∆) give contributions to F (y) that vanish like y or faster.
JHEP06(2017)083 Figure 6. Graphs for the short-distance behaviour of a DPD that involve a twist-four distribution.
x i and u i denote longitudinal momentum fractions.
One may also derive the expansion (3.18) from the operator product expansion, without taking recourse to the transverse momentum representation (3.17). In the definition of collinear DPDs one has a product O 2 (0, z 2 ) O 1 (y, z 1 ) of operators with z 1 = z 2 = 0 but nonzero y. This can be expanded around y = 0 in terms of light-ray operators where all fields are at transverse position 0. These operators have twist 2, 3, 4 for the first, second and third term in (3.18), respectively.
The spitting contribution F spl,pt is given by graphs as in figures 4 and 5a and has already been discussed in the previous subsection. The term F tw3,pt originates from two types of graphs. The first type involves a single perturbative splitting and a quasipartonic collinear twist-three distribution as shown in figure 5b. The second type has two splittings as in figure 4 and a twist-three distribution with one "good" and one "bad" parton field. Given the helicity constraints discussed in section 2.1, collinear twist-three distributions in an unpolarised proton involve a quark and antiquark with opposite chirality (and possibly an extra gluon). As announced earlier, we discard twist-three terms in the following, since they are expected to become unimportant at small momentum fractions x 1 , x 2 .
Finally, the term F int,pt contains contributions without any perturbative splitting; we hence refer to it as the "intrinsic" part of the DPD. It can be written as F int,pt (x 1 , x 2 , y; µ) = G(x 1 , x 2 , x 2 , x 1 ; µ) + C(· · ·, y; µ) ⊗ G(· · ·; µ) + . . . (3.19) where G is a quasipartonic collinear twist-four distribution and C a perturbative splitting kernel, corresponding to graphs as in figure 6. The convolution ⊗ is in the longitudinal momentum fractions indicated by · · · (cf. figure 6a). The first term in (3.19) corresponds to the first term in (3.17) and is the only contribution that does not involve a hard splitting at all. The ellipsis denotes terms with non-quasipartonic twist-four distributions containing three or two parton fields, together with one or two parton splittings. While having the same power behaviour in y, one may expect that at small x 1 , x 2 these terms become less important than the terms with quasipartonic twist-four distributions, which should roughly grow as fast as the square of two parton densities. We now take a closer look at the second term in (3.19). The kernel C can be determined by computing both sides of (3.19) for a given graph. At order α s only "non-diagonal" interactions, i.e. interactions connecting partons 1 and 2 as in figure 6a and b, contribute to C. The ladder graph in figure 6c is independent of y and thus gives identical contributions to the matrix elements F int,pt and G. As a consequence it does not contribute to C.

JHEP06(2017)083
At this point we can comment on the scale evolution of the different terms in (3.19). The l.h.s. evolves according to the homogeneous double DGLAP equation for DPDs, which describes "diagonal" interactions, either between the partons with final momentum fraction x 1 or between those with final momentum fraction x 2 . By contrast, the evolution of G(x 1 , x 2 , x 2 , x 1 ; µ) contains both diagonal and non-diagonal ladder interactions [34]. The non-diagonal interactions in the evolution must thus be cancelled by the µ dependence of the term C ⊗ G. At leading order in α s , this dependence comes only from the coefficient function C, which indeed contains just non-diagonal interactions as just discussed.
We finally emphasise that an unambiguous decomposition of F (y) into splitting, intrinsic and twist-three parts is only possible in the limit of small y. If y is of hadronic size, neither the operator product expansion nor the notion of perturbative parton splitting make sense. One may however use the short-distance decomposition (3.18) as a starting point for a model parameterisation of DPDs in the full y range. We describe a simple version of this strategy in section 9.

A scheme to regulate DPS and avoid double counting
In this section, we present a scheme that regulates the DPS cross section and solves the double counting problem between DPS and SPS, as well as between DPS and the twist-four contribution (figure 1b). Before doing so, we discuss some general considerations that motivate our scheme.
The following properties are in our opinion desirable for any theoretical setup describing double parton scattering.
1. It should permit a field theoretical definition of DPDs, without recourse to perturbation theory. This is the same standard as for the ordinary parton distributions in SPS processes. In particular, it allows one to derive general properties and to investigate these functions using nonperturbative methods, for instance lattice calculations.
One may object that so far not even ordinary PDFs can be computed to a precision sufficient for phenomenology. However, important progress has been made in the area of lattice computations, and more can be expected for the future. Furthermore, whereas ordinary PDFs are being extracted with increasing precision from experiment, it is hard to imagine a similar scenario for DPDs, because of their sheer number and because DPS processes are much harder to measure and analyse than most processes from which PDFs are extracted. In such a situation, even semi-quantitative guidance from nonperturbative calculations (such as the relevance of correlations of different types) is highly valuable.
As already discussed in section 2, the requirement of a nonperturbative definition prevents us from separating the "perturbative splitting" contribution of a DPD in a controlled way for all distances y (or equivalently for all conjugate momenta ∆).

2.
To pave the way for increased theoretical precision, the scheme should permit a formulation at higher orders in perturbation theory. Furthermore, the complexity of the required higher order calculations should be manageable in practice.

JHEP06(2017)083
3. For collinear factorisation, one would like to use as much as possible existing higherorder results for SPS processes, namely partonic cross sections and splitting functions. This means that the scheme should not modify the collinear subtractions to be made in hard scattering kernels, nor the validity of standard DGLAP evolution for DPDs in the colour singlet channel.
4. For TMD factorisation, it is desirable not to modify Collins-Soper evolution and the handling of rapidity divergences. This again allows one to re-use calculations done for SPS, although rapidity evolution for DPS is necessarily more complicated due to the complexities caused by colour [5,21].
5. One would like to keep procedures as similar as possible for collinear and TMD factorisation. This will in particular facilitate the computation of DPS processes at perturbatively large transverse momenta in terms of collinear DPDs [21], adapting the well known procedure for single Drell-Yan production [38].
In principle one can use dimensional regularisation to handle the UV divergences that are induced in the DPS cross section by the perturbative splitting mechanism, as is done with the UV divergences that arise in simpler situations such as single hard scattering. However, contrary to that case, the UV divergences discussed in section 2 arise not at the level of individual DPDs but only when two DPDs are multiplied together and integrated over y. This means that if one treats these divergences in dimensional regularisation, only the product of two DPDs is defined in D = 4 dimensions but not the DPDs separately. This possibility was explored in [16]. However, DPDs and their products remain nonperturbative functions at large y, which according to present knowledge cannot be reduced to simpler quantities in a model independent way. In practice, one therefore needs to model or parameterise them at some starting scale. This is more involved for the product of DPDs than for DPDs themselves, as is the practical implementation of scale evolution. We will come back to this scheme in section 8.
Ultraviolet regularisation. We define the regularised DPS cross section by multiplying the integrand in the DPS formula (2.6) for measured transverse momenta with Φ(y + ν) Φ(y − ν) and the integrand in the collinear DPS formula (2.7) with Φ 2 (yν). The function Φ(u) goes to 1 for u ≫ 1 and to 0 for u → 0, and we can restrict ourselves to the case where 0 ≤ Φ(u) ≤ 1 for all u. More specific requirements are given below.
Collinear and transverse-momentum dependent DPDs are defined as specified in section 2, without any modifications. Constructed from operator matrix elements, they contain both splitting and non-splitting contributions. They quantify specific properties of the proton and have a simple physical interpretation, with the same limitations as single parton densities. (We recall that a literal density interpretation of PDFs and TMDs is hindered by the presence of Wilson lines and of ultraviolet renormalisation.) Double counting subtraction. To treat the double counting between DPS and other contributions, we adapt the recursive subtraction formalism of Collins, which we JHEP06(2017)083 briefly sketch now (details are given in sections 10.1 and 10.7 of [18]). Consider a graph (or sum of graphs) Γ that receives leading contributions from a set of loop momentum regions R. An approximation for Γ is then given by In each term one integrates over all loop momenta. The operator T R applies approximations designed to work in momentum region R. Subtraction terms avoid double counting the contributions from smaller regions R ′ (regions that are contained in R).
In these subtraction terms one applies the approximations designed for R and those designed for the smaller regions. One can show [18] that C R Γ then provides a valid approximation in the region R and in all smaller regions, and R C R Γ gives a valid approximation of the graph in all relevant regions. All approximations discussed here are valid up to power corrections.
In our context, we have graphs in which a set of collinear partons split into partons that can be either collinear (as in DPS) or hard (as in SPS). A slight adaptation of the above formalism is required since we compute DPS using DPDs and a regulating function Φ that depend on transverse distances y rather than transverse momenta. A collinear splitting region R ′ then corresponds to large y and the corresponding hard region R to small y, but we keep the ordering of regions R ′ < R from momentum space when implementing (4.1). We will show in section 4.3 that our use of subtractions in position space is equivalent to the one in momentum space up to power suppressed effects.
The subtraction terms for the DPS region turn out to have a very simple form. They can be obtained by replacing the DPDs in the UV regularised DPS cross section with their appropriate limits for small y ± in TMD factorisation and for small y in collinear factorisation. Details will be given in sections 4.1, 4.2 and 5.
Criteria 1 and 5 above are obviously satisfied in this scheme. Regarding criterion 2, we note that the higher-order calculations required for the double counting subtraction terms are for the short-distance limit of DPDs, which involve much simpler Feynman graphs than the full scattering process. The introduction of a function Φ in the DPS cross section avoids an explicit modification of the definition of DPDs and thus respects criteria 3 and 4. In particular, the collinear DPDs F (x i , y) in transverse position space follow the homogeneous DGLAP evolution equation (2.3). Since for colour singlet DPDs, the evolution kernels are the familiar DGLAP kernels, the associated scale dependence in the cross section cancels by construction against the one of the hard cross sections computed with the same collinear subtraction as for SPS.
In our scheme, we have introduced an additional momentum scale ν to separate DPS from SPS and the twist-four contribution. In practical calculations one may take ν equal to the UV renormalisation scale µ in DPDs, but we find it useful to keep it separate in the general discussion. As a minimal requirement, ν must be of perturbative size, so JHEP06(2017)083 that the double counting subtraction terms remove all contributions from nonperturbative regions in the hard kernels of the SPS and twist-four cross sections, making their calculation consistent. By construction, the ν dependence in the physical cross section cancels between DPS and the double counting subtraction terms, up to higher orders in α s that are beyond the accuracy of the computation. We will come back to this issue in section 6.2.
Let us now take a closer look at the properties required for the function Φ.
• The introduction of Φ in the DPS cross section must not spoil the correct description of the physics in the region where the DPD approximations work. Specifically, we require that the modifications introduced by Φ in that region should be power suppressed in the large scale. This requires that Φ(u) must approach 1 sufficiently fast when u ≫ 1.
Anticipating the momentum space analysis in section 4.3, we demand that for collinear factorisation the Fourier transformation of 1 − Φ 2 exists in 2 − 2ǫ dim (for positive and negative ǫ). Therefore the integral du u 1−2ǫ 1 − Φ 2 (u) must converge for u → ∞. The corresponding criterion for TMD factorisation is obtained by replacing Φ 2 → Φ.
• Φ(u) at small u must regulate the UV divergences of the naive DPD cross section in 2 dimensions. This means that du u −1 Φ(u) must be integrable at u = 0 in TMD factorisation. For collinear factorisation we have the stronger requirement that with some δ > 0.
• To compute the double counting subtraction terms for SPS and the twist-four contribution, we must perform integrals over y in 2 and 2 − 2ǫ transverse dimensions, respectively, as we shall see in the next sections. We choose Φ such that the required integrals are known analytically. This is especially important for the twist-four contribution, where we must expand the result around ǫ = 0.
A suitable function for both collinear and TMD factorisation is the step function where γ is Euler's constant. This corresponds to a hard lower cutoff y > b 0 /ν in the y integration, where the constant b 0 ≈ 1.12 is taken to simplify certain analytical results. An alternative choice is Important integrals for these functions are given in table 1. Further possible choices for collinear factorisation are Φ 2 (u) = 1 − exp(−u p ) or Φ 2 (u) = u p /(1 + u p ) with p > 2; the corresponding integrals in the third and fourth row of table 1 can be performed and give Euler Γ functions.

Leading order analysis: collinear factorisation
In this section we show how our formalism works in collinear factorisation, concentrating on the leading order in α s . Following the procedure of all-order factorisation proofs (see for instance [18]), we rewrite individual Feynman graphs in a way consistent with the final factorisation formula. At the end of this procedure, the factors associated with longdistance physics in that formula, such as PDFs and DPDs, are expressed in terms of operator matrix elements and thus defined in a nonperturbative way. During our Feynman graph analysis we can separately consider splitting and intrinsic contributions to a DPD, F spl and F int . Likewise, we can consider separate contributions σ 1v1 , σ 2v1 and σ 2v2 to the DPS cross section, which respectively involve the combinations F spl F spl , F int F spl + F spl F int and F int F int , corresponding to graphs a, b and c in figure 1. In the final factorisation formula, only the full DPDs F and the full DPS cross section σ DPS will appear. For brevity, we write σ instead of dσ/(dx 1 dx 2 dx 1 dx 2 ) here and in the following.

Squared box graph
Let us start with the squared box graph in figure 1a, integrated over all transverse momenta in the final state. As discussed in section 2, it receives its leading contribution from the region q ∼ Q. The quark loops in the amplitude and its conjugate are then in the hard momentum region. This gives a part of the SPS cross section σ SPS and is computed as usual in terms of PDFs and the cross section for gg annihilation into two gauge bosons. Other graphs that contribute to SPS, such as the one in figure 2a, have no overlap with DPS. In the present subsection we restrict σ SPS and other terms in the cross section to refer to the double box graph only. The DPS region of the graph has near-collinear splittings in both the amplitude and its conjugate, which is only possible at q ≪ Q. This gives the 1v1 contribution σ 1v1 to the DPS cross section, which is power suppressed compared with σ SPS because we integrate over the transverse boson momenta. Since the double box graph is integrated over all q in σ SPS , we subtract the contribution of the DPS region in order to prevent double counting. (At leading power accuracy this is not strictly necessary, but we will shortly see that it is useful). According to (4.1) this contribution is obtained by applying to the graph the JHEP06(2017)083 approximations for DPS as well as those for SPS. This gives the DPS expression with the g → qq splittings computed perturbatively, the incoming gluons treated as collinear on-shell partons, and quark masses in the loop set to zero. These approximations are just what goes into the perturbative splitting approximation F spl,pt of a DPD, which is given in (3.12). The subtraction term thus reads Note that in σ 1v1,pt we use the short-distance approximation F spl,pt (even for large y, where it is not valid), whereas in σ 1v1 we use the un-approximated splitting part F spl of each DPD.
To ensure that σ 1v1 and σ 1v1,pt receive contributions only from the DPS region, one should take ν ≪ Q.
To compute (4.5) we can use the known form (3.12) for the splitting DPDs, with the kernels collected in [39]. With σ SPS being computed in fixed order perturbation theory, we make the same approximations in σ 1v1,pt , taking the renormalisation and factorisation scales in the splitting kernel and in the collinear PDF fixed. In D = 2 dimensions the splitting DPDs then depend on y like 1/y 2 . The y integration in (4.5) is thus readily performed, This is proportional to ν 2 , which is expected since the unregulated integral has a quadratic divergence in y and ν is the scale regulating this divergence. The relevant integral for different Φ(u) is obtained from the last row of table 1 by setting ǫ = 0. We anticipate that the subtraction term at higher orders in α s involves an integral as in (4.5), with splitting DPDs given by (3.14) and (3.16). Compared with (4.6) the integrand then has additional powers of log(yµ). The integrals can be obtained from those in the last row of table 1 by Taylor expanding around ǫ = 0.
The complete contribution of the double box graph to the cross section is As is well known, an appropriate choice for the factorisation scale µ in collinear factorisation is µ ∼ Q. In the same spirit we take ν ∼ Q rather than ν ≪ Q, although this extends the y integrals in σ 1v1 and σ 1v1,pt to values y ∼ 1/Q where the DPS approximation does not work. Both σ 1v1 and σ 1v1,pt are then of the same order as σ SPS rather than power suppressed. Let us see that this still gives the correct approximation of the overall cross section. Although σ SPS is naturally computed in momentum space, one can perform a Fourier transform w.r.t. the transverse momentum difference ∆ (defined as in figure 4 but for the graph of the cross section). Then σ SPS is a loop integral over the conjugate distance y, just as the two other terms in (4.7).
We now show how the subtraction formalism works in this case. In the region of small y ∼ 1/Q, one has σ 1v1 ≈ σ 1v1,pt , because the perturbative approximation of the DPD is JHEP06(2017)083 designed to work at short distances. The dependence on the cutoff function Φ(yν) cancels between these two terms. One is therefore left with σ SPS , which gives the appropriate description of the graph for y ∼ 1/Q. We note that the DPDs in σ 1v1,pt are computed at fixed order in α s , whereas in σ 1v1 they should be resummed using the double DGLAP equations. The cancellation between the two terms is therefore only up to higher orders in α s , which are beyond the accuracy of the result. We shall investigate this in detail in section 6.2.
Turning to the region of large y ≫ 1/Q, one finds that σ SPS ≈ σ 1v1,pt , because precisely in that region the DPS approximation is designed to work for the graph. We are therefore left with the DPS term σ 1v1 , which is exactly what we want to describe this region. The cutoff function Φ(yν) has no effect here, since we take ν ∼ Q and require Φ(u) ≈ 1 for u ≫ 1. The combination (4.7) thus gives the correct approximation for both small and large y.

2v1 graph
We now derive a description for the 2v1 graph in figure 1b. Integrated over transverse momenta in the final state, the graph receives leading contributions from small and large transverse momenta of the quarks on the splitting side (i.e. the top of the graph), corresponding to small or large transverse momenta of the produced gauge bosons. The region of small transverse momenta corresponds to DPS, whilst the contribution from large transverse momenta is correctly described by the twist-four mechanism introduced in section 2.
Whereas the SPS contribution in the previous subsection gives a finite result when computed as a perturbative two-loop graph, the twist-four contribution has a collinear divergence when integrated over the full region of phase space. To exhibit this divergence, we use dimensional regularisation and write As in the previous subsection, we restrict the meaning of σ tw4 and other terms in the cross section to the contribution from a single graph. The labels a 1 and a 2 indicate parton species and polarisation, and G is a collinear twist-four distribution. H includes overall factors and the squared amplitude of the hard scattering process qq + g → V 1 V 2 . Note that in the hard scattering graph one has q 1 = −q 2 = q. Momentum fractions are as shown in figure 7. The plus-and minus-momentum fractions of the produced bosons are respectively given by with x i andx i defined in (2.4). For the momenta of internal lines in the hard scattering subprocess, minus and transverse components are fixed by the final state kinematics, so that there is only a loop integration over the plus-momentum fractions v and v ′ . Note that by virtue of (4.9) both X i and X i have an implicit dependence on q 2 at given x i ,x i . An upper limit on q 2 follows from the requirements that X 1 + X 2 ≤ 1 and X 1 + X 2 ≤ 1. The expression (4.8) includes a contribution from DPS region, characterised by q i ≪ Q i and v, v ′ ≪ 1. Let us approximate the integrand in that region. We first perform a Fierz transformation for the Dirac indices of the upper quark lines at the gauge boson vertices and retain only the leading Dirac structures. After this, the cross sectionsσ i for qq annihilation appear in the expression. Furthermore we set v = v ′ = 0 in G and extend the integrations over v and v ′ to the full real axis (after which they are easily done using Cauchy's theorem). Approximating X i ≈ x i and X i ≈x i and collecting all factors, we obtain

JHEP06(2017)083
with the splitting kernel T from (3.4). If we drop the restriction q ≪ Q i , then the expression in the second line of (4.10) becomes the perturbative splitting approximation F spl,pt (x i , ∆ = 0) of the transverse-momentum integrated DPD. This DPD is evaluated at ∆ = 0 in order to fulfil momentum conservation at the gauge boson vertices, because the transverse momenta of the right moving partons are set to zero in the overall hard scattering kernel H. The q integral of (4.10) has an infrared divergence in D = 4 dimensions, which must be cancelled by the subtraction term for the DPS region. That term is obtained by applying the DPS approximations to the graph, in addition to those for the twist-four region. This gives again the form (4.10), but with an unrestricted integration over q, and a regulator of the corresponding ultraviolet divergence. To implement the regulator, we Fourier transform the collinear DPD, , and then multiply with Φ 2 (yν) as prescribed by our formalism. This gives

JHEP06(2017)083
Alternatively, one can obtain the subtraction term by starting with the contribution σ 2v1 of graph 1b to the DPS cross section and applying the additional approximations adequate for the twist-four region. Let us translate the latter from transverse momentum (where the twist-four approximations are formulated) to transverse position (where our UV regulator for the DPS cross section is local). On the l.h.s. of the relation one neglects ∆ in the upper part of the graph, replacing F spl (x i , ∆) by its value at y). The twist-four approximation implies computing the g → qq splitting in perturbation theory, with massless quarks and the gluons taken collinear and on shell. This corresponds to replacing F spl with F spl,pt . For the graph under discussion one has F int (x i , y = 0) = G(x 1 , x 2 , x 2 , x 1 ) according to (3.19) and thus obtains (4.11). We see that, as in the case of the 1v1 graph, the subtraction term is obtained from the DPS cross section by replacing the DPDs with their small y approximation at the appropriate order in α s . In section 5 we will find that this also holds at higher orders, where we have to take into account the second term in (3.19). Let us finally perform the y integral in the subtraction term (4.11). Collecting all y and ǫ dependent factors in (3.12), we have where Ω 2n = 2π n /Γ(n) is the surface of a sphere in 2n dimensions. For brevity we do not display the momentum fraction in the splitting kernel P . The y integral is infrared divergent in 2 dimensions and converges for ǫ < 0. Using the third row of table 1, we get in particular The contribution of the graph to the overall cross section is finally given by in analogy to (4.7). By the same argument as in the previous subsection, one finds that this combination reproduces the graph over the full range of y values, with the second term cancelling against the third term for y ∼ 1/Q and against the first term for y ≫ 1/Q. As we have just seen, the product of (µ/ν) 2ǫ with the 1/ǫ pole gives a logarithm log(ν/µ) in the subtraction term σ 2v1,pt . In turn, the q integral in the expression (4.10)

JHEP06(2017)083
of σ tw4 gives log(Q/µ), where Q comes from the upper integration limit and the µ dependence from the regulated infrared divergence. In the combination σ tw4 − σ 2v1,pt this gives log(Q/ν). The ν dependence is cancelled by a log(ν/Λ) from the y integral of the DPS cross section σ 2v1 , where ν regulates the ultraviolet divergence. Combining all terms, one obtains log(Q/Λ) in (4.15). With the choice ν ∼ Q, the large part of the logarithm is contained only in the DPS term.

Combining contributions
We can now add all contributions to the physical cross section, including 1v1, 2v1 and 2v2 graphs in DPS. We include the same regulator Φ(yν) in all DPD contributions, so that the splitting and intrinsic contributions to the individual DPDs can be added up to distributions F = F spl + F int that are defined by operator matrix elements as described in section 2.
2v2 graphs like the one in figure 1c are dominated by the DPS region, where transverse parton momenta are small, whilst regions with large transverse momenta are subleading. In impact parameter space this means that including the Φ regulator to suppress the region of small y does not significantly change σ 2v2 . We will see in section 4.3 that the effect of Φ in the 2v2 term is at the level of power corrections.
The master formula for the cross section with integrated transverse momenta is then where σ DPS = σ 1v1 + σ 2v1 + σ 2v2 includes all contributions to the DPS cross section and involves the full DPDs F . From the discussion in the previous two subsections it follows that the ν dependence of σ DPS cancels against the one of σ 1v1,pt + σ 2v1,pt , up to terms of higher order in α s that are beyond the accuracy of the calculation. As discussed in section 2.1 we neglect contributions involving twist-three functions in (4.16). The terms in (4.16) are now meant to include all contributing graphs, including graphs without a double counting issue. The hard scattering cross sections needed for the DPS term are available for many final states, and σ 1v1,pt can be readily obtained from existing results at leading order in α s . For double gauge boson production, the SPS cross section has been computed at high perturbative order, including the double box graphs as well as the first radiative corrections to them [40]. On the other hand, to the best of our knowledge, no calculations of σ tw4 have been performed so far. To leading logarithmic accuracy, one can however omit the terms σ tw4 − σ 2v1,pt in (4.16), because with the choice ν ∼ Q the large DPS logarithm log(Q/Λ) generated by 2v1 graphs is entirely contained in σ DPS , as discussed in the previous subsection.

Leading order analysis: TMD factorisation
We now show how our formalism works for TMD factorisation. As explained in section 2.1, the cross section for q 1 , q 2 ≪ Q receives leading-power contributions from DPS, SPS and from their interference. The 2v1 contribution is power suppressed in this case and will not be discussed further in this subsection. The discussion of leading momentum regions JHEP06(2017)083 is a bit more involved now, because the two quark loops in the double box graph can be collinear or hard independently of each other.
We start with the SPS/DPS interference, for which a graph is shown in figure 2b. Its contribution to the V 1 V 2 production cross section has the form for DPS in the amplitude or in its complex conjugate. Here D are the twist-three TMDs introduced in section 3.1 and the proportionality is up to kinematic and numerical factors.
We have different hard scattering amplitudes H, with α 0 and β 0 being gluons and the remaining partons quarks or antiquarks. The regions of small y + or y − are regulated by Φ(y + ν) and Φ(y − ν).
Notice that the graph in figure 2b has a leading contribution when the g → qq splittings in the conjugate amplitude become collinear, in which case we are back in the DPS region. From the first expression in (4.17) one should therefore subtract with F y − →0 given in (3.8) and represented in figure 5b. From the second expression in (4.17) one should subtract an analogous term with y + → 0. According to (3.7), the computation of the subtraction terms involves the integrals which are given in the appendix. We see in table 1 that for n = 0 both our choices for Φ give a logarithmic behaviour log(ν/q) when ν ≫ q, as expected for a regulated logarithmically divergent integral.
Let us now turn to the double box graph. It has leading contributions when the loops are both hard (SPS region) or both collinear (DPS region), or when one loop is hard and the other one collinear (region of SPS/DPS interference). The terms to be added to the SPS cross section for removing the contribution from the SPS/DPS interference regions (4.20) The first terms in the square brackets are obtained from (4.17) by taking the perturbative splitting form for all distributions, the one for D α 1 α 2 |α 0 (x i , y + , Z) being given in (3.10). The second terms are obtained from (4.18) and its analogue for y + → 0 by taking both splitting vertices perturbative, using the DPD given in (3.11). Notice that the subtraction terms (4.20) contain subtractions themselves, showing the recursive nature of the construction in (4.1). Finally, one should remove the contribution from the DPS region by adding where again both splitting vertices in each DPD are taken as perturbative. Note that in both (4.20) and (4.21) all proton matrix elements are expressed in terms of ordinary twist-two TMDs, as is the case for the SPS cross section. Adding up contributions and double counting subtractions, we finally obtain where we have again omitted differentials dx 1 dx 2 . . . for brevity. The DPS cross section σ DPS contains contributions with and without splitting in each of the DPDs, which are defined by operator matrix elements as in the case of collinear factorisation. Let us now discuss the large logarithms appearing in the different terms. The hard scattering kernel H α 0 β 0 in (4.17) has a DPS logarithm log(Q/q), and correspondingly there is a log 2 (Q/q) in the unsubtracted SPS cross section. We take ν ∼ Q, so that these logarithms are fully removed by the double counting subtractions. They reappear in the logarithmic integrals over y + and y − in the DPS cross section and in the SPS/DPS interference terms (4.17). In the full cross section (4.22) one then has a squared logarithm log 2 (Q/q) in the DPS term, and a single log(Q/q) in the subtracted SPS/DPS interference in the first line, whereas the subtracted SPS contribution in the last line has no leading logarithm. If q is in the nonperturbative region, one has log(Q/Λ) instead of log(Q/q).
At a practical level, the computation of the SPS/DPS interference requires twist-three TMDs as input, and the DPS term requires transverse-momentum dependent DPDs. In general, the modelling of these functions is very difficult and largely unconstrained. The situation is much better when q 1 , q 2 ≫ Λ: in this case one can compute these functions in terms of collinear ones [21].
On the perturbative side, the calculation of the hard scattering cross sections in the DPS term is easiest, while the SPS and interference terms require the full computation of box graphs. If one is satisfied with keeping only the leading double logarithm in Q/q, then it is sufficient to retain only the DPS contribution with ν ∼ Q.

Subtraction formalism in momentum space
In [18] it was proven that the subtraction formalism leading to (4.1) correctly approximates a graph up to power suppressed terms, provided that each approximant T R correctly reproduces the graph in its design region R, up to power suppressed terms. This proof works with regions in momentum space. However, our UV regulator Φ for DPS singularities is multiplicative in y space, which gives convolutions in transverse momenta that do not appear in the original Feynman graphs and hence are not part of the proof in [18]. We now show that this does not cause problems, so that our setup gives adequate approximations when transformed to momentum space, with corrections that are subleading in powers of Λ/ν. With ν ∼ Q these corrections are of the same order as other approximations in the factorisation formula.
We begin with the case of collinear factorisation, using dimensional regularisation for UV and IR divergences as usual. As mentioned earlier, we choose Φ such that the Fourier transform is finite for all values of ∆/ν. In the following we show that terms that depend on Ψ in the cross section are power suppressed and can hence be discarded when establishing the validity of the subtraction formalism. Let us first discuss σ 2v2 , which is proportional to The productσ 1σ2 of hard scattering cross sections does not affect our argument and is hence omitted. To keep the notation simple, we only display transverse-momentum arguments in this subsection. For simplicity we perform the power counting for D = 4 dimensions; changes in D = 4 − 2ǫ are by fractional powers and do not alter our conclusions. We then have for both ∆ ∼ Λ and ∆ ∼ ν, which are the two scales present in (4.24). 5 The large ∆ behaviour of F int follows from dimensional analysis (Λ 2 in the numerator comes from a collinear twist-four matrix element and 1/∆ 2 from the hard splitting kernel in the Fourier transformed version of (3.19)). The Ψ independent term in (4.24) receives its leading contribution of order Λ 2 from the region ∆ ∼ Λ, which is the design region of the DPS approximation. By contrast, the Ψ dependent term behaves like Λ 4 /ν 2 , with contributions JHEP06(2017)083 from both ∆ ∼ Λ and ∆ ∼ ν. We thus see that our UV regulator (which is not necessary in σ 2v2 as noted earlier) only gives changes suppressed by Λ 2 /ν 2 . For ν ∼ Q this does not degrade the overall accuracy of the calculation. Next we consider the contribution σ 2v1 − σ 2v1,pt + σ tw4 from 2v1 graphs. Here we have a subtraction term, which is obtained by replacing the DPDs in y space with their perturbative approximations. In momentum space this reads The second line is obtained by replacing F with F int in the first term and with F int,pt in the second term on the r.h.s. of (3.17) and then taking the Fourier transform with respect to y. We thus find that σ 2v1 − σ 2v1,pt is proportional to For power counting we use (4.25) and additionally which follows from dimensional analysis of the splitting graph. The perturbative approximations F int,pt and F spl,pt scale like their un-approximated counterparts. The Ψ independent terms in (4.27) are then found to be of order Λ 2 , which after multiplication withσ 1σ2 gives the same power behaviour in the cross section as σ tw4 . To analyse the Ψ dependent terms, we rearrange them as follows: The first line is suppressed by the integration phase space for ∆ ∼ Λ, by Ψ(∆−∆ ′ )−Ψ(∆) for ∆ ∼ ν and ∆ ′ ∼ Λ, and by F int (−∆ ′ ) − F int,pt (−∆ ′ ) for ∆, ∆ ′ ∼ ν. The second line is again suppressed by the integration phase space for ∆ ∼ Λ, and by F spl (∆) − F spl,pt (∆) for ∆ ∼ ν. The sum of all Ψ dependent terms is hence power suppressed compared with the Ψ independent terms. One readily finds that the relative suppression is by Λ 2 /ν 2 , as it was in the 2v2 case. Notice that this suppression is obtained for the combination σ 2v1 − σ 2v1,pt but not for the two terms separately.

JHEP06(2017)083
Finally, we discuss the contribution σ 1v1 − σ 1v1,pt + σ SPS from 1v1 graphs. The combination σ 1v1 − σ 1v1,pt is proportional to (4.30) One finds that both the Ψ dependent and the Ψ independent terms on the r.h.s. scale like Λ 2 . Multiplied withσ 1σ2 this is power suppressed by Λ 2 /Q 2 compared with σ SPS . Notice that the SPS approximation itself has relative corrections of order Λ 2 /Q 2 . This means that the contribution in (4.30) is beyond the accuracy of the calculation of 1v1 graphs. The fact that the Ψ dependent part is of similar size as the Ψ independent one is hence of no concern. We recall from section 2.1 that the rationale to evaluate the 2v2 and 2v1 terms, which are also suppressed by Λ 2 /Q 2 , is that they may be enhanced by other factors. The case of TMD factorisation can be discussed along the same lines. In analogy to (4.23) we define Ψ as the Fourier transform of [ Let us discuss the SPS/DPS interference graph in figure 2b. Its contribution to σ DPS/SPS in (4.17) is proportional to where on the r.h.s. we have twist-three TMDs in momentum representation. An analogous relation holds for the contribution to σ DPS/SPS, y + →0 , which according to (4.22) is to be subtracted from σ DPS/SPS in the overall cross section. In this case one should take the perturbative splitting approximation D spl,pt of the distributions, which in y + representation is given by (3.10). For power counting we take q 1 , q 2 ∼ Λ and use where k − and K can be of order Λ or ν. This can be derived by a general analysis of graphs as in section 5.2 of [5]; the 1/k − dependence is also directly obtained by Fourier transforming (3.10). D spl,pt behaves like D, and of course we have Ψ ∼ 1/ν 2 as before. We see that for k − ∼ K ∼ Λ the Ψ independent term in (4.31) is of order 1/Λ 2 , whilst the Ψ dependent one is of order 1/ν 2 . The only region in which the Ψ dependent term gives a leading contribution is when k − ∼ k ′ − ∼ ν and K ∼ Λ. However, in this JHEP06(2017)083 region one has D ≈ D spl,pt , so that its contribution is suppressed in the combination σ DPS/SPS −σ DPS/SPS, y + →0 . We thus find that in the overall cross section, Ψ dependent terms are power suppressed, as in the case of collinear factorisation. The same holds of course for σ SPS/DPS − σ SPS/DPS, y − →0 . Finally, a corresponding analysis (involving the momenta k + , k − and K) can be given for the combination σ DPS − σ DPS, y − →0 − σ DPS, y + →0 + σ DPS, y ± →0 , which covers all terms in (4.22) that depend on our UV regulator.

Subtraction terms at higher orders
Our construction of the double counting subtractions for DPS regions is not limited to the leading-order graphs discussed so far. In the present section we illustrate how the formalism works at highers in α s , taking as an example the 2v1 mechanism in collinear factorisation. In the master formula (4.1) we encounter nested subtractions in this case. We show for selected graphs how the formalism gives the different contributions to the cross section, namely σ tw4 , σ 2v1 and the subtraction term σ 2v1,pt . In particular we find that the latter can be obtained from σ 2v1 by replacing the DPDs with their short-distance expansions F int,pt and F spl,pt , as found at leading order in section 4.1.2.
We focus here on the general structure and use a schematic notation, in particular suppressing indices and arguments that are not essential for the discussion.
Graph 8a. To introduce this notation, let us briefly review the leading-order graph in figure 8a, which was discussed in detail in section 4.1.2. It has two leading momentum regions, which we denote by C and H. In region C the gauge bosons have small transverse momenta and the g → qq splitting is near collinear. This is the region of DPS, and its JHEP06(2017)083 approximated contribution to the cross section reads is the tree-level cross section for qq → V i , and we write y as a shorthand for d D−2 y. The contribution of this graph to the splitting DPD at the top is denoted by F In region H, the transverse momenta of the gauge bosons are large and thus the qq pair into which the gluon splits is far off shell on both sides of the final state cut. This is described by the twist-four contribution and requires a subtraction for the smaller region C. The subtracted contribution is Here q is shorthand for d D−2 q, the tree-level hard scattering cross section for qq + g → V 1 V 2 is denoted by H (1) , whilst f (0) and G (0) respectively are the zeroth order contributions to the gluon distribution at the top and the twist-four qq distribution at the bottom of the graph. As seen in (4.8), there is a convolution integral over longitudinal momentum fractions in G, denoted here by ⊗. In the subtraction term, V (1) represents the lowest-order perturbative g → qq splitting kernel, including all terms in (3.12) except for the singleparton density f a 0 . In this term, the momentum fractions of G (0) are fixed by external kinematics, so that the convolution ⊗ is absent. The perturbative expressions of the DPDs for small y are at leading order, so that the combined contribution of the C and H regions can be written as As already mentioned, the double counting subtraction in the second term is simply obtained from the DPS term by replacing each DPD with its small-y expansion at the appropriate order in α s .
Graph 8b. We now analyse the NLO graph in figure 8b. The qq splitting can be hard or collinear, and likewise the additional gluon with momentum ℓ. The leading momentum regions thus are CC, HC and HH, where the first letter refers to the qq splitting and the second to the additional gluon. There is also a region where the gluon is hard whereas the JHEP06(2017)083 qq splitting to the right of the cut is collinear. The splitting at the left of the cut is then hard by momentum conservation. One finds that this region is not leading: the loss of phase space for having the gauge bosons at low transverse momenta is stronger than the gain from having one but not both splitting vertices collinear. Let us discuss the leading regions in turn.
In the CC region, the extra gluon is part of the DPD for the right moving particles, and one simply has where the superscript (1) indicates a contribution to F int at order α s . Likewise, in the HC region, the extra gluon becomes part of the twist-four distribution G, so that one has ren . (5.7) Note that the integral over the gluon momentum in the graph giving G ren contains an ultraviolet divergence. It is understood that G (1) ren is renormalised, i.e. includes the appropriate UV counterterm, which is indicated by the subscript "ren".
For the HH region, the recursive character of the subtraction formalism becomes evident. We have a subtraction for the next smallest region HC, which itself contains a subtraction for the region CC. In addition, a CC subtraction has to be made for the overall graph. This gives The term in the first line is for the hard region, with H (2) being a one-loop contribution to the cross section for qq + g → V 1 V 2 . In the second line, the term for the HC region (including its subtraction for the CC region) involves a non-diagonal splitting kernel K ren (0) for the exchange of a gluon between the two parton lines with momentum fractions x 1 and x 2 , which are taken at relative transverse position y = 0 because this is part of the approximation for the HC region. In the last line, which is the overall subtraction for the CC region, the same kernel is to be taken at relative distance y, because in this case the qq splitting is not assumed to be harder than the extra gluon emission. One therefore has to take a convolution in transverse momenta, which turns into a product of two y dependent factors V (1) (y) and K (1) (y) after Fourier transformation.
Notice that in dimensional regularisation the limit y → 0 in K (1) (y) is not smooth. The graph for K (1) ren (0) contains a scaleless momentum integral, which is zero. After renormalisation, i.e. subtraction of the ultraviolet counterterm, one is left with an infrared JHEP06(2017)083 divergence. 6 By contrast, K (1) (y) at finite y does not have an ultraviolet divergence (and hence no counterterm), but it does have the same infrared divergence as K (1) ren (0). Adding all contributions and using (5.3) we obtain spl,pt (y) F int,pt (y) In the first line of (5.9) we have an NLO contribution to the twist-four term, with the NLO kernel H (2) having a subtraction for the region where the gluon becomes collinear, as is standard in collinear factorisation. In the DPS subtraction term we have a contribution to the small-y expansion of F int at NLO, with the three terms in (5.10) corresponding respectively to the three terms in the general expression (3.17). The infrared divergence of the kernel K drops out here, as it must. We see that the DPS subtraction term in (5.9) has the same form as in (5.5), with F int,pt instead of F int,pt . In (5.9) all ultraviolet and infrared subtractions have been carried out, so that one can set D = 4 in dimensional regularisation. The quantities with subscript "ren" then depend on the associated scale µ. In the first line, this dependence cancels between G Graph 8c. We next consider the graph in figure 8c. The contribution from the region CC is int,ren (y) , (5.11) where in contrast to graph 8b the one-loop contribution to F int now has an ultraviolet divergence (since the two partons connected by the gluon are at the same transverse position) and needs renormalisation. The contribution from region HC is analogous to (5.7). Unlike for graph 8b, a leading contribution is now also obtained for the momentum region CH, where the g → qq splittings are collinear, whilst the gluon carries a large transverse momentum (balanced by the vector boson momentum q 2 ). Including the subtraction for the region CC, we have 12) 6 The UV divergent part of K ren(0) gives DGLAP splitting kernels generalised to non-diagonal momentum fractions, known from the evolution of twist-four distributions [34] and of GPDs [41].

JHEP06(2017)083
where σ (1) 2 is the cross section for qq → V 2 + g and K (1) ren is a diagonal splitting kernel including the necessary UV counterterm. Its divergent part involves the contribution of this graph to the usual DGLAP splitting kernel P qq . Finally, the region HH gives where the last three lines are the subtraction terms for the regions HC, CH, CC, respectively. Adding all contributions, we obtain The latter relation holds because the exchanged gluon in graph 8c gives exactly the same contribution to F int (y) and to G and hence cancels between the second and third term in the general small-y expansion (3.17) of F int (y). In each line of (5.14) we have hard scattering cross sections subtracted for the collinear region, as is characteristic for collinear factorisation. The µ dependence cancels between the subtracted cross sections and the renormalised distributions.
Graph 8d. We now discuss figure 8d, where we have an additional gluon in the upper part of the graph. Its leading momentum regions are CC, HC and HH. As in the case of graph 8b, there is no leading region CH. For the CC region we simply have In the region HC, the extra gluon in the graph is collinear and couples to an off-shell quark, since the g → qq splitting is hard. Among the approximations to be made in this case is the Grammer-Yennie approximation (see e.g. [5,20]). After summing over an appropriate set of graphs, one can apply a Ward identity, after which the gluon couples to an eikonal line, as shown in figure 9. We thus can write

JHEP06(2017)083
ℓ + 3 more graphs = ℓ Figure 9. Graphical representation of applying the Grammer-Yennie approximation and a Ward identity to the sum of four graphs. The three other graphs required are graph e in figure 8, and the analogues to graphs d and e in which the gluon couples to the quark line at the very right. It is understood that the g → qq splitting is hard, whereas the gluons are collinear left moving.
where Γ e+ denotes the sum of graph e and of two graphs analogous to graphs d and e, with the gluon coupling to the rightmost quark line. f (1) ren is an NLO contribution to the gluon density, involving a gluon coupling to an eikonal line. The integral over the gluon transverse momentum diverges in the ultraviolet and hence requires renormalisation. For the HH region we have In the subtraction term for the region HC, we have an NLO contribution K (1) ren to the gg splitting kernel with a gluon coupling to an eikonal line, in analogy to f (1) ren above. Again, this requires summation over several graphs as specified on the l.h.s. In the last line of (5.18) we have the subtraction term for the region CC, in which the gluon and the splitting are regarded as collinear, so that no eikonal approximation is made and one obtains an NLO contribution V (2) to the g → qq splitting kernel. Adding all contributions, we obtain where F spl,pt (y) = d and e+ is a contribution to the perturbative splitting DPD at NLO. The scale dependence of f

JHEP06(2017)083
Other NLO graphs. The discussion of graphs where the extra gluon is exchanged between a right-moving and a left-moving line, such as in figure 8e and f, is more complicated because it involves more momentum regions. Depending on whether the g → qq splitting is collinear or hard, the extra gluon can be hard (H), collinear left moving (L), collinear right-moving (R), or soft (S). In the three latter cases, one can apply the Grammer-Yennie approximation and, after summation over an appropriate set of graphs, apply Ward identities, similarly to the example we have just discussed. This leads to the gluon coupling to eikonal lines, which in regions L and R are part of the Wilson lines in the operator defining the DPDs, whereas in region S they are part of the Wilson lines in the soft factor. For more details we refer to [20]. Our formulation of the DPS cross section in transverse position space is well suited for the treatment of these Wilson lines, which are taken at the same transverse positions as the quark or gluon fields in the DPD definitions. In transverse momentum space, one instead has a more complicated structure with convolution integrals.

Double DGLAP evolution
Collinear DPDs F (x i , y; µ i ) evolve with µ i according to the double DGLAP equations. In the present section, we discuss several aspects of this scale dependence.

Resummation of large logarithms
As we have seen in section 4.1.2, the 2v1 graph in figure 8a gives rise to a large logarithm log(Q/Λ) when the cross section is integrated over the transverse boson momenta. This logarithm originates from the g → qq splitting vertices in the kinematic region of double parton scattering and is correctly reproduced in the DPS cross section, provided that we make a suitable choice of scales.
Radiative corrections to figure 8a give rise to additional large logarithms of DGLAP type. One can expect that these logarithms are correctly resummed to all orders by the DGLAP evolution of the DPDs, provided again that one makes suitable scale choices. In this section we show for explicit examples that this is indeed the case. We will allow for different hard scales Q 1 , Q 2 and set the scales as where we define Of course, our conclusions regarding large logarithms will not change if (6.1) is modified by factors of order unity. For definiteness, we will always take a rigid cutoff Φ(u) = Θ(u − b 0 ) in yν. For later use we define the scale DGLAP logarithms are generated by loop integrals in regions of strongly ordered transverse momenta. We will compare the logarithms obtained directly from transverse momentum integrals with the ones obtained from evolved DPDs in y space. Figure 10. Radiative corrections to the hard scattering graph in figure 8a that give rise to DGLAP logarithms.

JHEP06(2017)083
To set the scene, we start with the lowest-order 2v1 graph in figure 8a. In the DPS region, i.e. when the scale of the g → qq splittings is much smaller than Q 1 and Q 2 , the behaviour of a splitting vertex is easily recovered from the computation of the splitting DPD. As can be seen from (3.4), a vertex with quark momentum k gives a factor k j /k 2 , where j is a transverse index. The transverse momentum of the splitting vertices must flow through the gauge boson lines. Thus, the two splitting vertices in figure 8a give a combined 1/q 2 1 behaviour for q 1 ≪ Q min , as is evident in (4.10) after angular integration. The integrated cross section is then proportional to Since the q 1 integral is logarithmic, we must include an infrared cutoff Λ, which limits the integration to the perturbative region where the preceding analysis is valid. We take here the usual procedure for extracting the leading behaviour of a logarithmic integral b a dk Γ(k), approximating the integrand Γ(k) for a ≪ k ≪ b while integrating over the full range a ≤ k ≤ b.
Let us now see how the logarithm (6.4) arises in the DPS cross section formula. According to the perturbative expressions of the DPDs, F spl,pt ∼ 1/y 2 and F int,pt is y independent at leading order in α s , so that we have (6.5) With our choice for the infrared cutoff on y, the expressions (6.4) and (6.5) agree including their overall normalisation. Now consider the NLO graph in figure 10a. We route the gluon momentum ℓ back through the boson line with momentum q 2 , so that ℓ flows through exactly two quark propagators. In the strongly ordered region the graph has the same DPS logarithm as the leading-order result, whilst developing an additional logarithm from the integral dℓ 2 /ℓ 2 , as is familiar from similar gluon emission JHEP06(2017)083 graphs in hard scattering processes. In the case Q 2 < Q 1 this gives whereas for Q 1 < Q 2 we have In both cases, we have an additional logarithm log(Q 2 min /Λ 2 ) compared with the lowestorder result. If Q 1 ≪ Q 2 , then the logarithm of the ratio Q 2 /Q 1 is also large.
Let us compare these results with the 2v1 part of the DPS cross section formula when we take into account evolution of the splitting DPD. Focusing on the dependence on y and µ, we can write at leading order, where the choice of scale µ y ensures that there are no large logarithms of yµ in the higher-order corrections. This should be evolved to F spl,pt (x i , y; µ 1 , µ 2 ) in the cross section. In (6.9) there is a logarithmic dependence on y from the running of α s and from the scale dependence of the gluon density. Graph 10a does not contribute to either of these, so that the corresponding dependence should not be taken into account here. Rather, the DGLAP logarithm of graph 10a corresponds to the evolution of (6.9) from µ y to µ 2 , which gives a factor proportional to α s log(µ 2 2 /µ 2 y ). The relevant y integral in the DPD cross section thus is where we have used the scale choice in (6.1). This agrees with both (6.7) and (6.8). Note the ease with which (6.10) is obtained in y space, compared with the detailed considerations leading to (6.7) and (6.8). We see that with evolved DPDs in σ DPS one includes higher order logarithms while using tree-level partonic cross sections. Treating the 2v1 graphs fully within the collinear twist-four formalism would require the explicit inclusion of the higher-order graphs.
Whereas in the previous example, the gluon is emitted close to the splitting vertices, the opposite is true for the graphs in figure 10b and c. For figure 10b the region generating two large logarithms is

JHEP06(2017)083
To see how these logarithms arise from the evolution of F int,pt , we use the short-distance expansion F int,pt (x 1 , x 2 , y; µ y , µ y ) = G(x 1 , x 2 , x 2 , x 1 ; µ y ) (6.13) at leading order, which we evaluate at scale µ y to avoid large logarithms of yµ in the higher-order corrections, as we did for F spl,pt in (6.9). Graph 10b contributes to the scale evolution of G, where at leading order it gives α s log(µ 2 y /Λ 2 ) when evolution is started from the generic soft scale Λ 2 . The graph also contributes to the evolution of F int,pt from µ y to µ 2 , where it gives α s log(µ 2 2 /µ 2 y ). The two effects add up to a factor α s log(µ 2 /Λ 2 ). In the DPS cross section we thus have in agreement with (6.12). The graph in figure 10c also has a region where DGLAP logarithms are generated, as was already discussed in [17]. We route the loop momentum ℓ through two quark lines at one of the two g → qq splitting vertices. In the strongly ordered region the two lower fermion lines to which the gluon couples have virtualities of order ℓ 2 , whereas both splitting vertices have virtualities of order q 2 1 ≈ q 2 2 . This gives 7 Just as graph 10b, graph c contributes to the evolution of G in (6.13), giving a factor α s log(µ 2 y /Λ 2 ) at leading order. There is however no contribution to the evolution of F int,pt from µ y to µ 2 from this graph, so that in the cross section we have which again reproduces the result (6.16) obtained in momentum space, including the overall factor. A fully analogous discussion can be given for the graph in figure 8b. One finds the same type of DGLAP logarithms as in (6.16), which is again reproduced correctly in the DPS cross section formula. Note that in this case the additional gluon is a virtual correction rather than emitted into the final state -a somewhat unusual situation in the context of DGLAP logarithms.
In conclusion, higher-order graphs with strong ordering of transverse momenta produce DGLAP type logarithms in 2v1 graphs, which can depend on Q min /Λ or Q 2 /Q 1 . Our JHEP06(2017)083 Figure 11. expression for the regularised DPS cross section correctly reproduces these logarithms if we make an adequate choice of DPDs and scales. Specifically, we must • use DPDs that have the correct short-distance behaviour, as discussed in section 3.3. Fixed-order perturbative results for this behaviour can be used reliably at scales of order µ y , as in (6.9) and (6.13).
• take scales µ 1 , µ 2 and ν in σ DPS as specified in (6.1), with possible modifications by factors of order unity.
The situation is quite different for 1v1 graphs. In the DPS region which disappears however when one integrates over q 1 up to Q min . This is inevitable since the integrated 1v1 cross section depends on a single dimensionful quantity Q min and hence cannot have any large logarithm. An infrared cutoff Λ in the q 1 integration does not play any role here since the cross section is dominated by large q 1 ∼ Q min , where strong ordering as in (6.18) is not possible. One finds that the higher-order graph in figure 11b does have a DGLAP logarithm in the strong-ordering regime due to the additional loop integral dℓ 2 /ℓ 2 . But the region (6.20) gives only a power suppressed contribution to the cross section integrated over q 1 , where DPS kinematics is not dominant. 8 There is an exception to the statement just made. If at least one of the momentum fractions x 1 , x 2 is very small and the corresponding scale µ 1 , µ 2 is sufficiently large, then JHEP06(2017)083 the evolution of F spl,pt from µ y to µ 1 , µ 2 can significantly flatten the 1/y 2 behaviour of the initial condition (6.9). The underlying physics is radiation of soft gluons, which strongly enhances parton distributions at low x, as is well known from the evolution of ordinary PDFs. The enhancement by evolution of F spl,pt increases with the distance between µ y and µ 1 , µ 2 , i.e. it increases with y. In the 1v1 (and also the 2v1) contribution to the DPS cross section, this enhances the importance of the region y ≫ 1/Q, where the DPS approximation is valid. In this region, the evolution of the DPDs correctly describes the DGLAP logarithms generated by graphs like the one in figure 11b. In the numerical studies presented in section 9.2, we find that this happens for F gg , F qg and Fq g , as well as for F qq when one momentum fraction is much smaller than the other.

Dependence on the cutoff scale
As already discussed in section 4.1, our formalism relies on the fact that in the region of small y ∼ 1/Q the subtraction terms in the cross section (4.16) cancel against the DPS contribution, up to corrections that are beyond the accuracy of the calculation. Such corrections include higher-order terms in α s , whose parametric size we now investigate.
Let F fo (x i , y; µ) denote the short distance approximation of a DPD, computed in fixedorder perturbation theory. It depends on µ via the scales of α s and of PDFs or twist-four distributions (and at higher orders also via powers of log(yµ) in the splitting kernels). In the cross section, this fixed-order approximation is used in two different ways: 1. For the small-y limit of the DPDs used in the DPS cross section, σ DPS , one takes F fo (x i , y; µ y ) as initial condition and evolves it to the scales µ 1 ∼ Q 1 and µ 2 ∼ Q 2 of the two hard scatters.
2. In the subtraction terms, σ 1v1,pt and σ 2v1,pt , one takes F fo (x i , y; µ h ) evaluated at the same scale µ h that is used in the hard scattering terms, σ SPS or σ tw4 , so as to cancel the contribution from the region y ≫ 1/Q order by order in perturbation theory.
The scale µ h may be a combination of Q 1 or Q 2 or an independent third scale. For now we assume that α s L ≪ 1, where The case where this is not satisfied because the hard scales are widely different is discussed in the following subsection. Consider now a calculation at N k LO (k = 0 for LO), using DPDs evolved with N k LO kernels P and N k LO expressions for the DPS subprocess cross sectionsσ 1,2 and the shortdistance approximations F fo . We have and its analogue forσ 2 . Here ⊗ denotes the usual Mellin convolution for the indicated variable, and we omit labels and sums for parton species. To understand (6.22) we observe

JHEP06(2017)083
that the l.h.s. would be exactly zero ifσ 1 and P were calculated to all orders; the convolution ofσ 1 with the relevant parton densities would then be µ 1 independent. With N k LO truncations on the l.h.s., the first term beyond the accuracy of the calculation is of order α k+1 sσ1 . Furthermore, the highest power of log(µ 1 /Q 1 ) inσ 1 is k according to the general pattern for logarithms related with the renormalisation or factorisation scale. In a similar fashion we can derive (6.23) using that k is the highest power of log(yµ) in the N k LO approximation, according to (3.16). Finally, it is easy to see that for a quantity satisfying where the integration of the differential equation (6.24) between µ a and µ b gives an extra logarithm on the r.h.s. We thus obtain for the small-y contribution to the 1v1 cross section using our initial condition F spl,pt = F spl,fo at scale µ y . The longitudinal momentum arguments and convolutions are as in (2.7). The first line of the last expression in (6.26) is just the double counting subtraction term. We thus find with L defined in (6.21). In the second step we have used that the full SPS cross section is of the same order of magnitude as the DPS subtraction term, regarding its power behaviour in Q 2 as well as the relevant coupling constants and PDFs. In full analogy one can derive

JHEP06(2017)083
Notice that the difference σ tw4 − σ 2v1,pt is by construction dominated by short distances y ∼ 1/Q, whilst the individual terms receive contributions from a wide range of y. We thus find that, for both 1v1 and 2v1 contributions, the DPS term cancels against the double counting subtractions in the region y ∼ 1/Q, up to corrections of relative order α k+1 s . This is beyond the accuracy of an N k LO calculation, as it should be. The only logarithms multiplying α s in this context are logarithms of the ratios between different hard scales, which at this point we must require to be of moderate size.
We can now also quantify the dependence of the overall cross section σ in (4.16) on the scale ν. With two choices ν 1 ∼ ν 2 ∼ Q the difference σ(ν 1 ) − σ(ν 2 ) is just of the size given by (6.27) and (6.28). In analogy to the variation of renormalisation and factorisation scales, one may thus use the variation of ν to estimate uncalculated higher-order corrections.
An important consequence of (6.27) is that where we used that the subtraction term σ 1v1,pt is dominated by y ∼ 1/Q, given the 1/y 2 behaviour of F spl,fo . This means that we can use the ν variation of the 1v1 contribution to σ DPS to estimate the size of the subtraction term σ 1v1,pt and thus of the SPS contribution to the cross section at the corresponding order in α s . This is of great practical value, given the considerable difficulty to compute SPS at this order. We will come back to this point in section 9.2.

Multiscale problems
The arguments leading to (6.27) and (6.28) are based on an expansion in α s and break down when α s L is not small, for instance because there is a large hierarchy Q 1 ≪ Q 2 between the two hard scales. We now show that a modified choice of DPDs in the subtraction terms σ 1v1,pt and σ 2v1,pt can cope with such a situation. Let us for definiteness describe the 1v1 subtraction; the 2v1 case can be treated in full analogy. We use the notation F spl,pt (x i , y; µ b , µ c |µ a ) for the distribution that is initialised to F spl,pt (x i , y; µ a , µ a |µ a ) = F spl,fo (x i , y; µ a ) (6.30) and evolved according to the double DGLAP equations to the scale µ b (µ c ) for the parton with momentum fraction x 1 (x 2 ). The small-y limit of the distributions used in σ DPS is F spl,pt (x i , y; µ 1 , µ 2 |µ y ) and its analogue for F int,pt . Let us now introduce a profile function p(u; µ a , µ b ) that interpolates monotonically between the limiting values µ a at u = 0 and µ b at u → ∞, satisfying We can then take distributions F spl,pt x i , y; p(yν; µ 1 , µ h ), p(yν; µ 2 , µ h ) | p(yν; ν, µ h ) (6.32)

JHEP06(2017)083
in σ 1v1,pt . These distributions tend to F spl,fo (x i , y; µ h ) for y ≫ 1/ν, so that the subtraction term cancels the contribution of σ SPS in this region. For y ∼ 1/ν the distributions tend to F spl,pt (x i , y; µ 1 , µ 2 |ν) ≈ F spl,pt (x i , y; µ 1 , µ 2 |µ y ), so that in that region the subtraction term cancels against σ DPS . Distances y ≪ 1/ν are of course removed in the cross section by the cutoff function Φ(yν). Using the same technique as in section 6.2, one can show that the cancellations just described are up to corrections of order α k+1 s at N k LO, without any large logarithms. In particular, this ensures that the ν dependence of the overall cross section is beyond the accuracy of the computation.
As an alternative, one may replace ν in (6.32) with min(µ 1 , µ 2 ), which requires less computation when ν is varied in the cross section (while µ 1 , µ 2 , µ h are kept fixed). With this choice, the DPDs in the subtraction term reduce to F spl,fo (x i , y; µ h ) in the case of equal hard scales, Profile functions for the renormalisation scale have been used in other contexts, see e.g. [44,45]. An suitable function for our purpose is where u 1 ∼ 1, u 3 ≫ 1 and u 2 = 1 2 (u 1 + u 3 ). Functions that are piecewise polynomials have been used for instance in [45,46]. To estimate theoretical uncertainties, one can vary the profile function p, for instance by varying u 1 and u 3 . For a detailed discussion we refer to [46].

Approximation of the intrinsic distribution
As discussed earlier, the correct small-y limit F int,pt at leading order in α s is obtained by making the ansatz (6.13) at scale µ y and then evolving F int,pt to the scales µ 1 , µ 2 relevant for the DPS cross section. In practical terms, this requires an ansatz or model for the twist-four distribution G at some reference scale µ 0 and subsequent evolution of G to the scale µ y . The technical implementation of twist-four evolution is rather involved, and one may wish to avoid it in situations where the corresponding loss of precision is tolerable.
In the model used later in this paper we indeed follow this path, taking as limiting form for small y at some fixed reference scale µ 0 and then evolving F int,pt to the final scales µ 1 , µ 2 . For the r.h.s. we will take a form proportional to the product f (x 1 , µ 0 ) f (x 2 , µ 0 ) of ordinary PDFs (this model is restricted to DPDs in the colour singlet channel). Multiplying with a suitable function of y, one readily obtains a model for the intrinsic DPD at any value of y. This is indeed standard in the existing literature. As follows from the discussions below (3.19) and below (6.13), evolution of the ansatz (6.34) to scales µ 1 , µ 2 correctly takes into account diagonal interactions as in figure 6c whilst neglecting non-diagonal interactions as in figure 6a and b, which in the correct JHEP06(2017)083 treatment contribute to the evolution of G from µ 0 to µ y . For σ 2v2 , these interactions are irrelevant since the small-y region is power suppressed. Furthermore, we can expect them to be of limited importance for σ 2v1 in the situation described at the end of section 6.1. Evolution then leads to a y dependence of F spl,pt that is significantly flatter than the 1/y 2 behaviour obtained at fixed order. This shifts the dominant region of the y integration in σ 2v1 to larger values, so that on the F int,pt side there remains little room for the evolution of G between µ 0 to µ y . The limited numerical importance of non-diagonal interactions for 2v1 diagrams in such situations was also observed in [17], where calculations were done in transverse momentum space.

Collinear DPDs in momentum space
As we argued in section 4, there are several advantages to work with DPDs in y space when setting up a factorisation formula for the overall cross section. We now consider collinear DPDs depending on the Fourier conjugate transverse momentum ∆, showing that they can be connected with our formalism in a useful way. The motivation for this is twofold.
• Sum rules for momentum space DPDs at ∆ = 0 have been proposed in [9], where it was stated that they are preserved by the leading-order inhomogeneous evolution equations in [6][7][8][9][10]. If it can be shown that they are valid at some starting scale, this will provide valuable (albeit nontrivial) constraints on DPDs. Of course, no such sum rules can hold for d 2 y F (x i , y) since this integral is divergent at small y.
• Previous work on DPS, in particular on the perturbative splitting contributions [12][13][14][15], is formulated in transverse momentum space. The results derived below will facilitate the comparison of that work with ours, which is the subject of section 8.
Both aspects only concern colour singlet DPDs, to which the rest of this section is restricted. Colour non-singlet DPDs have an additional rapidity dependence described by Collins-Soper equations, which are multiplicative in y space (see [21]) and hence more complicated in momentum space. The function Φ that regulates the DPS cross section in our formalism can also be used to introduce momentum space DPDs. We define omitting the superscript indicating the colour singlet channel for brevity. With the conditions on Φ stated in (4.2), the integral converges at y = 0. Throughout this section, we restrict ourselves to the two choices of Φ in table 1 and assume that ν ≫ Λ. The dependence of F Φ on ν is then given by

JHEP06(2017)083
The derivative of Φ(yν) strongly suppresses the nonperturbative region y ∼ 1/Λ in the integral (or removes it altogether if Φ is a step function). This allows us to replace F (x i , y, µ) with its perturbative splitting part on the r.h.s. The contribution from F int,pt (x i , y, µ) is suppressed by Λ 2 /ν 2 . Inserting the expressions (3.12) or (3.14) for F spl,pt (x i , y, µ) at ǫ = 0 one obtains at leading order and For brevity, we have omitted labels and sums for parton species and polarisations. For ∆ = 0 we have I 0 = 1 and recognise in (7.3) the inhomogeneous term of the evolution equation for a conventional momentum space DPD F (x i , ∆, µ) defined by MS subtraction. At leading order, our momentum space DPD F Φ (x i , 0, µ, µ) thus obeys the same evolution equation as F (x i , ∆, µ). At small ∆ the inhomogeneous term is multiplied by I 0 (∆, µ, µ) = 1 + O(∆ 2 /µ 2 ). The higher-order version (7.4) at ∆ = 0 has the same structure as the inhomogeneous term in the evolution equation proposed in [10], although the evolution kernel may not be the same.
Let us now compute the momentum space analogue of F spl,pt (x i , y, µ), both with the cutoff regularisation (7.1) and with MS renormalisation of the splitting singularity. This gives the correct result for the full momentum-space DPD in the limit ∆ ≫ Λ, when the exponent e iy∆ suppresses large y ∼ 1/Λ in the Fourier integral, so that the perturbative splitting approximation can be used.
Restricting ourselves to leading order, we start with (3.12) and use the Fourier integrals given in (A.3) and (A.4). For ∆ ≪ ν we find with corrections of O(∆ 2 /ν 2 ), whereas for ∆ ≫ ν we obtain To compute the usual momentum space DPD, we take the Fourier transform (2.9) of F spl,pt (x i , y) in D = 4 − 2ǫ dimensions, subtract the UV divergence and then set ǫ = 0. This gives where we have abbreviated x = x 1 + x 2 and v = x 1 /(x 1 + x 2 ). Here C ǫ denotes the MS subtraction term for the splitting singularity. Using (A.1) we get where P ′ (v, ǫ) = ∂P (v, ǫ)/∂ǫ. The MS counterterm is equal to so that we obtain We can use the preceding results to compute the matching between the DPD defined with a Φ cutoff and the one obtained by MS renormalisation. Since the integral in (7.1) is finite, it can be smoothly continued to D = 4 − 2ǫ dimensions. We then have The factor [1 − Φ(yν)] strongly suppresses the region y ≫ 1/ν, so that on the r.h.s. we can use the fixed-order perturbative splitting contribution, evaluated in 4 − 2ǫ dimensions. This gives for all ∆, and we have a smooth limit for ∆ → 0. For ∆ ≪ ν we can use (7.6) for Φ(u) = Θ(u − b 0 ) with (7.12) and obtain At ∆ = 0 we can in particular achieve equality of the two DPDs by choosing the scale With the splitting kernel P g→qq from (3.13) this gives for instance a scale varying between ν = µe −1/2 ≈ 0.6µ and ν = µ. These values respectively pertain to v = 1/2 (i.e. x 1 = x 2 ) and to v = 0 or 1 (i.e. x 1 = 0 or x 2 = 0).

JHEP06(2017)083 8 Comparison with other work
Three groups have previously studied the consequences of the perturbative splitting mechanism on DPS cross sections, namely Blok et al. [12,13], Ryskin and Snigirev [14,15], and Manohar and Waalewijn [16]. We refer to them as BDFS, RS and MW, respectively. In this section we compare the results obtained in our formalism with theirs. We restrict ourselves to collinear factorisation. To begin with, we rearrange our master formula (4.16) as and recall that one should take ν ∼ min(Q 1 , Q 2 ) for the cutoff parameter. Let us briefly characterise the different terms in (8.1). In all approaches just mentioned, the SPS cross section σ SPS is defined in the standard way. It is obtained by multiplying two PDFs with a hard scattering cross section and summing over all parton channels.
The term σ tw4 − σ 2v1,pt does not have any equivalent in the work of BDFS and of RS. As discussed in section 4.1.2, this term is lacking the large log(Q/Λ) that is contained in σ DPS and built up by parton splitting taking place over a wide region 1/Q ≪ y ≪ 1/Λ in the 2v1 graph 1b. The approaches of BDFS and RS can thus only be valid at leading logarithmic accuracy. This limitation is explicitly stated in the work of BDFS [12]. The same limitation holds of course in our approach if σ tw4 − σ 2v1,pt is not computed.
The term σ DPS − σ 1v1,pt includes the 2v2 and 2v1 parts of DPS, whereas the 1v1 part has its small y approximation subtracted. As discussed at the end of section 6.1 and illustrated in section 9.2, there are situations where strongly ordered parton emissions in 1v1 graphs enhance the region y ≫ 1/Q for the initial parton splitting. In our approach, this effect is included via the evolution of the DPDs in σ DPS (which thus provides a small x enhancement for σ 2v2 , σ 2v1 and σ 1v1 ).
If such an enhancement does not take place, then the 1v1 part of DPS has an integrand going like 1/y 4 at small y, so that the integral is dominated by the region y ∼ 1/Q. In σ DPS − σ 1v1,pt the contribution from this region then cancels up to higher orders in α s , as shown in section 6.2. In turn, the contribution from y ∼ 1/Λ is suppressed by Λ 2 /Q 2 compared with σ SPS and thus of the same order as other power corrections beyond the accuracy of the calculation. In this situation, the dominant terms in σ DPS − σ 1v1,pt are the 2v2 and 2v1 parts of DPS, provided that they have a small x enhancement compared with generic power corrections to SPS.

The approach of Manohar and Waalewijn
The master formula for the overall cross section in the MW approach is given in eq. (14) of [16]. Adapted to our notation, it reads whereĉ 3 andĉ 4 are hard scattering coefficients and [· · ·] indicates UV renormalisation of the enclosed quantity. This means that the product F (y)F (y) of DPDs is integrated over JHEP06(2017)083 y in 2 − 2ǫ transverse dimensions, and the resulting poles in ǫ are subtracted as usual in dimensional regularisation (we assume MS renormalisation, although this is not explicitly stated in [16]). The terms going withĉ 3 involve the renormalised product of a twist-two distribution in one proton and a twist-four distribution in the other. In the following, the terms withĉ 4 andĉ 3 in (8.2) will be called DPS and twist four terms, respectively. Via the perturbative splitting mechanism in one of the two protons, the operator in the DPS term mixes under renormalisation with the operators in the twist-four terms. This is discussed at LO in [16], but we see no obstacle to extending the analysis to higher orders. In section 4.1.2 we have seen that the computation of the hard scattering coefficient for graph 7 has a divergence from the region where the g → qq splittings become collinear. In our formalism the term σ 2v1,pt in (8.1) subtracts the contribution from this region. In the MW approach, the general subtraction formalism recalled in section 4 will yield an MS subtraction term for the collinear divergence, resulting from the renormalisation procedure described in the previous paragraph. For the choice Φ(u) = Θ(u − b 0 ) we see in (4.14) that σ 2v1,pt is proportional to the MS expression 1/ǫ + log(4π) − γ plus an extra term log(µ 2 /ν 2 ) + P ′ (0)/P (0), where P (ǫ) is the relevant splitting function in 4 − 2ǫ dimensions. This extra term corresponds to the difference between the momentum space DPD F Φ,spl,pt in (7.6) and its MS counterpart in (7.12). Up to this difference, the term σ tw4 − σ 2v1,pt in our master formula (8.1) is thus equivalent to the twist four term in (8.2). The mismatch between the two versions may be understood as a scheme difference.
It remains to compare the DPS term in the MW approach with our term σ DPS −σ 1v1,pt . As the scale evolution of d 2 y F (y)F (y) derived by MW reflects perturbative splitting in one of the two protons, their DPS term should contain the contributions from 2v2 and 2v1 graphs. Our approach also contains the contribution of 1v1 graphs, with a subtraction of their small y approximation. As discussed above, this part is negligible if DPS evolution effects are moderate. Both approaches should then give a valid representation of the overall cross section, with the possibilities to go beyond the leading logarithmic approximation and to systematically include higher orders in α s . On the other hand, strongly ordered emissions can significantly enhance the region y ≫ 1/Q in 1v1 graphs. In our formalism, such emissions are resummed to all orders by the evolution of the DPDs, but we see no counterpart of this in the MW formulation.
MW argue that under renormalisation d 2 y F (y)F (y) does not mix with the product of two PDFs by pointing out that the 1v1 graph in figure 1a gives a quadratically divergent scaleless integral, which is zero in dimensional regularisation. This argument hinges on using the perturbative approximation for the g → qq splitting at all values of y (or equivalently of the Fourier conjugate momentum ∆). This includes the infrared region, where the perturbative approximation does not work. It would be interesting to investigate this point further, but this is beyond the scope of the present paper.

The approaches of Blok et al. and of Ryskin and Snigirev
Since BDFS and RS use the same σ SPS as we do, and since they have no counterpart to our σ tw4 − σ 2v1,pt (which can be neglected at leading logarithmic accuracy), it remains to JHEP06(2017)083 compare our term σ DPS − σ 1v1,pt with their expressions for the DPS cross section. In this comparison, we can take the leading logarithmic approximation of our result.
The approaches of BDFS and of RS are both based on separating the intrinsic and splitting contributions to a DPD at all values of ∆. As discussed in section 4, we avoid such a separation when formulating the factorisation formula for the cross section. When modelling the DPDs -which is of course unavoidable for phenomenology -it is however convenient to use this separation. In y space we thus make the ansatz for all y. We demand that F int and F spl separately follow the homogeneous evolution equation (2.3) and that for y ≪ 1/Λ they tend to the correct short-distance limits F int,pt and F spl,pt . At small y the form (8.3) is hence unambiguous, whereas for large y it represents a model. For F int,pt we make the approximation discussed in section 6.4. Then F int has a finite limit for y = 0 rather than a weak singularity induced by non-diagonal twist-four evolution. Notice that F int and F spl now denote two terms in a DPD model, whereas in previous sections they were associated with different Feynman graphs in the discussion of factorisation. BDFS and RS only consider DPDs for unpolarised partons in the colour singlet sector, neglecting spin and colour correlations in both the intrinsic and the splitting parts. We will therefore do the same when comparing to our formalism. We note that spin correlations tend to be washed out by evolution to high scales [39] and that colour correlations are suppressed by Sudakov logarithms [22,47]. For sufficiently large hard scales, the importance of such correlations will therefore be reduced.
In terms of the momentum space DPDs introduced in section 7, the DPS cross section at leading order in α s is given by times the appropriate subprocess cross sections and a combinatorial factor. As was done in the papers cited above, we allow for different factorisation scales µ 1 , µ 2 . We note that with the behaviour (7.7) of the splitting DPD at large ∆, the integral (8.4) converges in the ultraviolet for both choices of Φ. For definiteness, we take Φ(u) = Θ(u − b 0 ) in the following. We now analyse the 2v2, 2v1 and 1v1 contributions in turn. The ultraviolet cutoff can be neglected in the intrinsic contribution: This is because F int (x i , y) is of generic size 1/Λ 2 for y ∼ 1/Λ and has a smooth behaviour for y ≪ 1/Λ. A cutoff y > b 0 /ν in the integral therefore only removes a contribution of order Λ 2 /ν 2 . If we make the same ansatz for the intrinsic DPD as BDFS or RS, we will thus obtain the same result for σ 2v2 as they do, up to power corrections that are beyond the accuracy of all approaches discussed here.

JHEP06(2017)083
For a typical ansatz, F Φ,int (x i , ∆) strongly decreases when ∆ ≫ Λ. The 2v1 contribution to the integral (8.4) is therefore dominated by the region ∆ ∼ Λ. In this region, the integral receives leading contributions from both small and large y. However, a large logarithm log(ν/∆) is only built up in the region 1/ν ≪ y ≪ 1/∆, where we can use the perturbative form F spl,pt (x i , y) ∼ 1/y 2 . To leading logarithmic accuracy, we can therefore compute the splitting contribution without a nonperturbative model, even for ∆ ∼ 1/Λ, which we now do. We want to improve on our result (7.6) by including evolution effects. To do so, we start from the first order perturbative expression (3.12) at scale µ y = b 0 /y and then evolve using the Green function D(x, µ; µ y ), which solves the single DGLAP equation in the scale µ with the initial condition D(x, µ y ; µ y ) = δ(1 − x). For simplicity we omit labels for and sums over parton species. We then have The function E(x i ; µ i , µ y ) depends weakly on y due to logarithmic corrections from evolution and the running of α s . To make this explicit, we expand the quantities depending on µ y (the Green functions, α s and the parton density) around ν. Keeping only leading logarithms, we obtain We show in appendix A that with such a y behaviour, the Bessel function J 0 (y∆) in (8.7) can be replaced by an upper cutoff Θ(b 0 − y∆), up to relative corrections in ∆ 2 /ν 2 and up to subleading logarithms α n+1 s (ν) log n−2 (ν/∆). We then have for ∆ ≪ ν. With the scale choice (6.1) our result (8.10) corresponds to the expressions of the perturbative splitting DPD in equation (18) of [12] and in equation (8) of [15]. Note that at leading logarithmic accuracy, the lower integration limit on k in [12] is equivalent to ∆ in the region ∆ ∼ Λ that dominates the 2v1 cross section.
To leading logarithmic accuracy, the results for σ 2v1 obtained by both BDFS [12,13] and RS [14,15] are therefore consistent with our approach (provided of course we make a JHEP06(2017)083 suitable ansatz for F int ). In equation (6) of [15] an upper cutoff min(Q 1 , Q 2 ) is put on the ∆ integration, but this is of no concern since as discussed above, the 2v1 term is dominated by ∆ ∼ Λ.
The situation is different for σ 1v1 , which RS compute using the form (8.10) we just derived for ∆ ≪ ν. The integral over ∆ in the cross section then requires an upper cutoff, which RS choose as min(Q 1 , Q 2 ), as we do for our cutoff ν in y space. The dominant integration region for the 1v1 cross section is ∆ ∼ ν in their approach and in ours, given that up to logarithms the splitting DPD is of order unity for both ∆ ∼ Λ and ∆ ∼ ν. For ∆ ∼ ν our distribution F Φ,spl differs from (8.10) even if we neglect evolution effects, which readily follows from the result (A.3) for the y integral in (8.7). Furthermore, F Φ,spl strongly depends on the choice of Φ in that region. This just highlights that σ 1v1 has no physical meaning by itself. In our formalism, the contribution from ∆ ∼ ν in σ 1v1 is cancelled by the subtraction term σ 1v1,pt , as it should be since the DPS approximation is not valid in that region at all.
The DPS cross section in the approach of BDFS is given entirely by 2v2 and 2v1 terms. To leading logarithmic accuracy, these terms agree with the ones in our approach. The term σ DPS − σ 1v1,pt in our approach has an additional contribution σ 1v1 − σ 1v1,pt , which as already mentioned may be neglected unless DPD evolution enhances the region y ≫ 1/Q (or equivalently ∆ ≪ Q) in σ 1v1 . This possible enhancement is not captured by the BDFS approach.
In the RS approach, the DPS cross section includes 2v2, 2v1 and 1v1 terms, the first two of which agree with ours in the leading logarithmic approximation. If DPD evolution effects are moderate, then the 1v1 term of RS is dominated by large ∆ ∼ Q, strongly depends on the upper cutoff imposed on the ∆ integral, and leads to a double counting problem with σ SPS that is not solved in the RS approach. In the sense we have discussed in section 6.2, one may at best take the cutoff dependence of the 1v1 cross section as an estimate for the size of σ SPS and thus as an indicator for how serious the double counting problem is. Only if DPD evolution enhances the region ∆ ≪ Q so strongly that it dominates the 1v1 term does the RS approach give a valid representation of the overall cross section. The cutoff dependence of the 1v1 cross section is then reduced, as was emphasised in [14].

Quantitative illustrations
We have set up a general formalism for calculating and combining the contributions from SPS, DPS and other mechanisms to the physical cross section. Whilst a detailed phenomenological application is beyond the scope of this work, we wish to give some quantitative illustrations of the formalism. This will allow us to estimate the relative size of contributions in specific channels and kinematics, and to explore the interplay between SPS or DPS terms and the double counting subtractions. We limit ourselves to collinear factorisation, where the construction of a model for DPDs is much simpler than in the TMD case.

JHEP06(2017)083
Convenient quantities for our purpose are double parton luminosities which directly appear in the DPS cross section. For definiteness, we use a sharp cutoff Φ(u) = Θ(u−b 0 ) to remove the UV region. As a model for colour singlet DPDs we take the sum of an intrinsic and a splitting part, as specified in (8.3) and in the paragraph following that equation. The luminosity (9.1) then naturally splits into separate terms for 2v2, 2v1 and 1v1. Colour non-singlet channels will not be considered, except in subsection 9.3, where they can be computed without recourse to a model.

Simplified analytic estimates
We shall see that DPD evolution has important quantitative effects on double parton luminosities. Nevertheless, let us start by deriving simple analytic expressions where evolution is neglected, but where the power behaviour, logarithms and other important factors can be easily identified.
For the intrinsic contribution we make the well-known product ansatz with a Gaussian form for the y dependence. For the splitting contribution we take which is the perturbative expression at O(α s ) times an exponential to suppress the large y region. Sums and labels for the parton type are omitted here. For simplicity we use the same exponential in F int and F spl . Since we ignore evolution for the time being, no scale dependence is taken into account in (9.2) and (9.3). To simplify even further, let us take equal momentum fractions x i =x i = x in both protons. The parton luminosity for the 2v2 contribution then is As mentioned earlier, the UV regulator Φ(yν) only gives rise to a power suppressed effect. The 2v1 luminosity reads

JHEP06(2017)083
This is enhanced over L 2v2 by log(ν 2 /Λ 2 ) but suppressed by α s /(2π). Furthermore, there are different factors regarding the parton densities. We note that the size of the splitting function strongly depends on the parton channel, as is exemplified by the cases P g→qq (1/2) = 1/4, P q→qg (1/2) = 10/3 and P g→gg (1/2) = 27/2. For the 1v1 luminosity we finally obtain We see the power enhancement over the 2v2 and 2v1 terms, as well as the quadratic dependence on the cutoff ν.
It is easy to repeat the above exercise for our alternative regulator function Φ(u) = 1−e −u 2 /4 . The integrals can be performed writing Φ(u) = u 2 1/4 0 dτ e −τ u 2 . For the leading term of the expansion in Λ 2 /ν 2 , one obtains the same result for L 2v2 as in (9.4), whereas in (9.5) and (9.6) one should replace respectively. The leading logarithm log(ν 2 /Λ 2 ) in L 2v1 is independent of the choice for Φ(u), but not the constant term. The corresponding ambiguity is compensated by the DPS subtraction term, as is the change of the numerical prefactor in L 1v1 .

Collinear parton luminosities
We now present a numerical study of the parton luminosities, using a slightly refined model for the input DPDs, and including evolution effects. We start by considering the unpolarised luminosities only, turning to the issue of polarised DPDs and luminosities at the very end of the section.

Setup
We use scale evolution at leading order in α s and take equal scales µ 1 = µ 2 = µ for the two hard processes throughout. In this illustrative study, we fix the number of active quark flavours to n f = 3. This is sufficient for exploring the broad qualitative features of the DPDs and luminosities. Implementation of heavy quark flavours, as required for a realistic phenomenological investigation, would require some dedicated work on flavour thresholds and on memory management in our numerical code, which we postpone to future studies. As in section 9.1, we write each DPD as a sum of intrinsic and splitting pieces: F a 1 a 2 (x 1 , x 2 , y; µ, µ) = F a 1 a 2 ,int (x 1 , x 2 , y; µ, µ) + F a 1 a 2 ,spl (x 1 , x 2 , y; µ, µ) . (9.8) We initialise F int at a low scale µ 0 = 1 GeV using a product ansatz similar to (9.2):

JHEP06(2017)083
For the single PDFs f a (x, µ 0 ), we take the MSTW2008LO distributions [48]. In (9.9) we have multiplied the PDFs by a function of the x i that does not affect the DPD at small x i , but smoothly cuts it off near the kinematic boundary x 1 + x 2 = 1. The function we use is that given in equation (3.12) of [9], with n = 2. We also take a different y dependence for different parton species. For this we use a simplified version of the model in section 4.1 of [39], taking the width h to be x independent (corresponding to h(x 1 , x 2 ) of [39] evaluated at x 1 = x 2 = 10 −3 ) and setting each h with q − indices to be the same as the one with q + . Thus we have For the splitting piece of the DPD we generalise our ansatz in (9.3), choosing an initialisation scale that goes to b 0 /y at small y but freezes to a constant value b 0 /y max when y exceeds a value y max that marks the transition between perturbative and nonperturbative behaviour. This ensures that the single PDF and α s in the splitting expression are never evaluated at too low scales. A suitable choice of scale is Such a prescription is very similar to the b * prescription used in TMD phenomenology [38,49]. Here we take y max = 0.5 GeV −1 , which is one of the values considered in the recent TMD study [50]. Using the same parton dependent Gaussian damping as in (9.13), we have F a 1 a 2 ,spl (x 1 , x 2 , y; µ y , µ y ) The coupling α s (µ y ) is determined by 3 flavour running, starting with the MSTW2008LO value α s (µ 0 ) = 0.68183. The single PDFs are obtained by taking the MSTW2008LO distributions at µ 0 and evolving them according to the DGLAP equations for n f = 3.
To obtain the splitting and intrinsic DPDs at scale µ, as in (9.8), the input forms just discussed must be evolved, starting from µ 0 for F int and from µ y for F spl , according to the homogeneous double DGLAP equations. For this we use a modified version of the code developed in [9]. The modified code works on a grid in the x i , µ and µ y directions (the grid of the original code is in x i and µ only). The grid points in the x i directions are evenly spaced in log(x i /(1 − x i )), whilst those in the µ and µ y directions are evenly spaced in log µ or log µ y . The integrals appearing in the double DGLAP equations are computed from points in the x i grids using Newton-Cotes rules (for details see appendix A of [9]), and evolution from one point in the µ grid to the next is carried out using the Runge-Kutta method.

JHEP06(2017)083
The grid is set up with 88 points in each x i direction, spanning 5 × 10 −5 < x i < 1, with 60 points in the µ direction, spanning µ 0 < µ < 170 GeV, and with 60 points in the µ y direction, spanning b 0 /y max < µ y < 340 GeV. According to the studies made in [9], this suggests an error on the level of a few per cent in the DPD values obtained after evolution, which is tolerable in this first study. The DPDs computed on this grid are used together with an interpolation code to produce numerical values for the investigations below.

Numerical results
We begin with a study where the scale µ is equal to Q 1 = Q 2 = 80 GeV (as in the production of a W boson pair). Taking the collider energy to be √ s = 14 TeV, we set x 1 andx 1 in (9.1) to correspond to central production of the first system and x 2 andx 2 to correspond to the production of the second system with rapidity Y (all rapidities refer to the pp centre of mass). This gives . (9.14) In figure 12 we plot L a 1 a 2 b 1 b 2 (Y ) in the range 0 ≤ Y ≤ 4 for the parton combinations a 1 a 2 b 1 b 2 = uūūu +ūuuū, a 1 a 2 b 1 b 2 = gggg and a 1 a 2 b 1 b 2 = uddu +duud. The first parton combination appears e.g. in ZZ production, the second is important in four-jet production, and the last appears in W + W + . For ease of language, we will refer to these parton combinations as the uū, gg and ud channels, respectively. We split the overall luminosity into contributions from 1v1 (F spl × F spl ), from 2v1 (F spl × F int + F int × F spl ) and from 2v2 (F int × F int ). We vary ν by a factor of 2 up and down around a central value of 80 GeV, in order to see how DPS alone is affected by variation of this cutoff. For each contribution, the line in the plots denotes the luminosity with ν = 80 GeV, whilst the band is generated from the envelope of the functions with ν = 40 GeV, ν = 80 GeV and ν = 160 GeV. In the case of the 2v2 contribution, the ν scale variation is negligible (as expected from basic considerations, see (9.4)), so this appears as a dashed line in each plot.
For the 1v1 contribution of the gg and uū channels, we also plot a band generated by varying the prediction with ν = 80 GeV by a factor 4 up and down. This corresponds to a strictly quadratic cutoff dependence of L, i.e. to the variation of ν by a factor 2 in the naive formula (9.6), where DPD evolution is neglected. Any discrepancy between this band and the actual 1v1 band is therefore due to evolution effects. We do not plot this band for the ud channel: there is no LO splitting process giving ud, so that the scale variation (and central value) from (9.6) is zero in this case.
We immediately notice in figures 12a and b that the 1v1 contribution is generally much larger than the 2v2 and 2v1 contributions, and that it has an enormous ν variation. To obtain a sensible prediction in these cases, one must include the SPS corrections up to the order that includes the double box graph in figure 1a, so that the associated subtraction term can approximately cancel the ν dependence of DPS. We also notice that at central rapidities in the uū channel, the 1v1 ν variation band essentially fills the band corresponding to a quadratic ν dependence, indicating small evolution effects in this case. By contrast, the 1v1 ν variation band for gg is clearly smaller than the naive expectation, which indicates significant evolution effects. We shall explore the reasons for these differences below.

JHEP06(2017)083
In figure 12c, the 1v1 contribution is small compared to the 2v2 piece and has a small ν dependence. This is because, as already mentioned, there is no LO splitting directly giving ud (generation of a ud pair requires at least two steps, such as u → u + g → u + d +d ). In this case, there is less of a need to compute σ SPS and the subtraction term σ 1v1,pt up to the order that contains the lowest-order DPS-type loop (in both amplitude and conjugate). This is fortunate, since this order is two powers of α s higher than for graph 1a (two-step rather than one-step splittings are required in both protons), and the corresponding SPS calculation will not be available for some time.
In the gg channel, the 2v1 and 2v2 contributions are of similar magnitude and shape. This is consistent with what has been found in previous phenomenological studies [13,51] at similar scales and x values. In the uū channel, the magnitude of 2v1 and 2v2 is broadly the same but their shape is somewhat different. In the ud channel, the shape of the two contributions is similar, but 2v1 lies well below 2v2 (and is close to 1v1).
We can examine the degree to which the ν variation in the gg and uū channels is reduced by adding the remaining fixed-order terms in the cross section, even without considering a specific final state or having to compute σ SPS . This is because only σ DPS and the subtraction terms depend on ν. Just as we did for σ DPS , we can factor the two partonlevel cross sections out of σ 1v1,pt and introduce a "subtraction luminosity" L 1v1,pt . This is defined as in (9.1) except that the DPDs are replaced by the fixed-order splitting expression evaluated at scale µ, For a given parton combination a 1 a 2 b 1 b 2 , one can directly subtract L 1v1,pt from the DPS luminosity L DPS and compare the ν dependence of the result to that of L DPS alone. Up to multiplication by the appropriate parton-level cross sections, the ν dependence of L DPS − L 1v1,pt is equal to that of the full cross section. On the other hand, the overall magnitude of L DPS − L 1v1,pt (which can be positive or negative) has no direct significance, since one needs to add the SPS contribution to obtain the physical cross section. We plot L DPS and L DPS − L 1v1,pt for the uū and gg channels in figure 13. The lines are the values for ν = 80 GeV, whilst the bands show the variation when ν is varied by a factor of 2 up and down. Note that for L DPS − L 1v1,pt , the curve with ν = 80 GeV is always at the very top of the scale variation band. We see that the ν scale variation of L DPS − L 1v1,pt is indeed reduced compared to L DPS , with the reduction being much stronger for uū than for gg. The latter is consistent with our previous observation that evolution effects are weaker in the uū channel at central rapidities: if evolution effects are weak, F spl and F spl,fo have a similar y dependence, so that (F spl × F spl − F spl,fo × F spl,fo ) is flat in y space and L DPS − L 1v1,pt varies weakly with ν.
In section 6.2 we argued on generic grounds that the ν variation of the 1v1 contribution to σ DPS should be of the same order as the subtraction term σ 1v1,pt . This is confirmed in figure 13, where the size of L 1v1,pt (evaluated at central ν) can be read off from the distance between the central curve for L DPS and the top of the band for L DPS −L 1v1,pt . This finding allows us to sharpen the argument we already made in the discussion of figure  does a large ν variation of the 1v1 contribution to σ DPS indicate the need to include SPS and the associated subtraction term in the cross section, but the size of the ν variation may even serve as a rough estimate for the size of these terms.
To better understand the behaviour we have seen in the 1v1 luminosities for the uū and gg channels, we now take a closer look at how evolution affects the y dependence of splitting DPDs. In figure 14 we plot F spl (y; µ) against y, setting x 1 = x 2 = 5.7 × 10 −3 so that the uū and gg DPDs are the ones used for the point Y = 0 in figures 12 and 13. We also plot the ug DPD, which mixes with uū and gg under evolution. For comparison we show F spl,fo (y; µ), as well as the initial condition F spl,ini (y) = F spl (y; µ y ) given in (9.13), from which F spl (y; µ) is obtained by evolution. The distributions are plotted in the range b 0 /(160 GeV) < y < b 0 /(40 GeV), i.e. in the range over which the lower integration limit in y is varied when we vary ν between 40 GeV to 160 GeV. In this region, the exponential damping factor in our DPD model is irrelevant, so that F spl,ini and F spl,fo only differ by the scales taken in α s and f a 0 (x 1 + x 2 ).
For the initial conditions F spl,ini we note that the uū and ug distributions are of similar size; the former is initialised by a larger PDF (f g instead of f u ) but has a smaller splitting coefficient P (1/2) as we noted before (9.6). By contrast, the gg distribution is much bigger; here both the initialising PDF and the splitting coefficient are large.
By construction, all three curves in each plot are equal at y = b 0 /(80 GeV), when µ y = µ. In all plots, F spl,ini has a more shallow y behaviour than F spl,fo . This difference is mainly driven by α s , as the PDFs do not vary so much between µ y and µ at momentum fraction x 1 + x 2 = 1.14 × 10 −2 . One observes that the difference between F spl and F spl,ini is more significant for gg and ug than for uū, i.e. that DPD evolution has a much stronger effect on the former two channels. This is to be expected at small x, since the 1/v behaviour JHEP06(2017)083 of the splitting kernels P gg (v) and P gq (v) at small v favours the radiation of a gluon with much smaller momentum than its parent parton, whereas the kernels P qq (v) and P qg (v) giving a quark stay finite for v → 0. An interesting point to note is that, whilst for gg and ug the curves for F spl are more shallow than for F spl,ini , in the case of uū the F spl curve is actually steeper. The latter is surprising since -based on the experience with single PDFs at small x -one may expect that forward evolution for y > b 0 /(80 GeV) would always increase a DPD, and backward evolution for y < b 0 /(80 GeV) would always decrease it. This indeed happens in the gg and ug channels, whilst forward evolution results in a decrease of the uū DPD. The reason for this is that F ug,spl and F gū,spl , which directly feed into F uū,spl during evolution, are comparatively small. In the case of PDFs, f g is much larger than f u and hence can JHEP06(2017)083 drive its small-x evolution although the splitting function P qg does not favour the radiation of low-momentum quarks. The evolution of the gg and ug DPDs is driven by the large distribution F gg,spl and enhanced by the P gg splitting function.
Note that in all three channels, F spl has a smaller slope in y than F spl,fo . This implies that the y integrand for the computation of L DPS − L 1v1,pt is positive for y > b 0 /(80 GeV) and negative for y < b 0 /(80 GeV). Therefore, the L DPS − L 1v1,pt curves for ν = 40 GeV and ν = 160 GeV lie below that for ν = 80 GeV. The curve for the central ν value thus lies at the top of the ν variation band, as we already observed in figure 13.
For uū, the evolution from µ y to µ turns out to have much the same quantitative effect as adjusting the scale in the fixed order expression from µ y to µ, such that F spl and F spl,fo end up extremely close together. This explains why at central rapidities the ν variation of the 1v1 uū contribution in figure 12a coincides almost exactly with the naive prediction from (9.6), which assumes that F spl depends on y like y −2 (as does F spl,fo ). Furthermore, if F spl and F spl,fo are very close to each other, then L DPS − L 1v1,pt is much smaller than L DPS , which we indeed see in figure 13a.
In the gg channel, the modification of the y slope by evolution is significant. In the range of figure 14, the y dependence is changed from y −2 at the starting scale to approximately y −1.4 . One may expect this flattening effect to become even stronger as the x values in the DPDs decrease. The dependence of L DPS on the scale ν should then decrease. This was already anticipated in [14,15], where studies in the double leading logarithm approximation were performed.
Here we investigate this effect using our full LO DGLAP set-up. We consider the gg channel with all x values set equal and study the DPS luminosity as a function of the common x value. We fix the collider energy √ s, so that the scale µ = Q 1 = Q 2 = x √ s also varies with x. We make one plot with √ s = 14 TeV and 5 × 10 −4 < x < 0.2 (corresponding to 7 GeV < µ < 2800 GeV) and another one with higher collider energy √ s = 100 TeV and 10 −4 < x < 0.02 (corresponding to 10 GeV < µ < 2000 GeV). In our numerical code this requires a new DPD grid with larger µ and µ y ranges, which was generated using 60 points in the µ and µ y directions (as in the original grid). The resulting 1v1 luminosity is shown in figure 15, with bands for the ν variation and its naive version as in figure 12. For comparison we also show the 2v1 and 2v2 luminosities. We see that the ν variation does indeed become progressively smaller compared to the central value as x (and µ) decreases. At the lowest x values it is much smaller than the naive expectation. At larger x, where evolution has a smaller effect on the DPDs, the ν variation becomes larger: for √ s = 14 TeV the actual and naive ν bands essentially coincide towards the right of the plot, where values of the order of x ∼ 0.1 are reached. Towards the left of each plot, where x and µ become small, the 2v2 and 2v1 contributions begin to dominate over 1v1. For given µ, this effect is more pronounced for √ s = 100 TeV than for √ s = 14 TeV, since the x values in the former case are smaller. We note that in the uū channel (not shown here) the actual and naive ν bands remain very similar throughout the kinematics of figure 15.
For the small x values where the ν variation is small in figure 15, the gg splitting DPD is significantly shallower in y than the fixed-order y −2 form. For In such a regime, the bulk of the contribution to the 1v1 part of L DPS comes from large distances y ≫ 1/ν, where the DPS approximations used to derive σ DPS are valid. Also, the 1v1 part of σ DPS becomes clearly larger than its ν variation. As we argued earlier, the ν variation should be of the same size as the subtraction term. This is confirmed in figure 16, which is analogous to figure 13b but for the kinematics of the point with µ = 10 GeV in figure 15a. Using our argument JHEP06(2017)083 that σ 1v1,pt is of the same magnitude as the order of σ SPS that contains the lowest-order DPS-type loop, one may in this situation justifiably make predictions that include the DPS piece but omit the order of SPS just specified, as well as the associated subtraction term. This is encouraging, for instance in the context of four-jet production, where the computation of the relevant SPS order (namely NNLO) is well beyond the current state of the art. Note that lower orders in SPS should be computed and included, if possible.
Notice that when the y behaviour of F spl becomes flatter than y −1 , the dominant y region in the 1v1 part of L DPS shifts to values y ∼ 1/Λ, where one cannot compute the splitting DPD in perturbation theory and must rely on a model. Likewise, the 2v1 part of L DPS becomes increasingly sensitive to the region y ∼ 1/Λ as soon as F spl becomes flatter than the fixed-order form y −2 .
Another kinematic regime where the 1v1 contribution to L DPS becomes large compared with its ν variation and compared with L 1v1,pt is when the two hard systems have a large separation in rapidity. This reduction of the ν variation can be seen as Y approaches 4 in figure 12, but only in the uū and not in the gg channel. In both channels, the effect becomes more pronounced once the rapidity separation of the hard systems is increased beyond 4. To illustrate this, we make plots similar to figure 12 but now with one hard system at rapidity Y and the other at rapidity −Y (rather than one at Y and the other at 0), such that a given value of Y corresponds to a rapidity separation of 2Y . The results are shown in figure 17, where we see that the ν variation in the 1v1 contribution is strongly reduced towards the right hand side of the plots, becoming hardly visible in the uū channel. Also notable is the fact that for large Y , the 2v2 contribution in this channel exceeds 1v1, which strongly decreases between Y = 0 and 2.
This reduction in ν dependence can be explained by the fact that at large Y we probe splitting DPDs with one large x and one small x parton. From the point of view of small x logarithms, it is preferable to generate such a configuration by having the 1 → 2 splitting JHEP06(2017)083 at large x, generating directly the large x parton plus a gluon with smaller x, the latter of which splits in a number of stages into smaller x gluons, eventually yielding the small x parton. This increasingly happens with increasing y, since the "evolution distance" between the initial and final scales, µ y and µ, is increased. Thus, evolution again flattens the DPD compared to the naive y −2 expectation and reduces the ν dependence of the DPS luminosities. The effect is particularly drastic in the uū channel, because the lowest-order splitting g → uū is inefficient at generating a pair with very different x values (as P g→qq (v) goes to a constant at small v). Therefore, the repeated splitting mechanism described in the previous paragraph is strongly preferred, even though its last step is penalised by the lack of small-v enhancement in P gq (v). In figure 18 we plot the uū DPD in the x 1 ≪ x 2 configuration relevant for Y = 4 in figure 17a. It has a similar y dependence as the ud DPD in equal kinematics, with the curve for ud crossing zero exactly at y = b 0 /µ (by construction), and the curve for uū crossing zero rather close to this point. In the former case, the lowestorder 1 → 2 splitting process is forbidden, whereas in the latter it is numerically almost irrelevant. The situation for x 1 ≫ x 2 is the same.
To end this section, we study polarised distributions and luminosities. We limit ourselves to the splitting part F spl of the DPD and to the 1v1 contribution to L DPS . In fact, we have little guidance for modelling the intrinsic part F int , where a product ansatz as in (9.9) makes no sense. We note that in [39], different ansätze were made for polarised DPDs, and the effects of evolution on these distributions were investigated. It was found that, in many cases, evolution to higher scales leads to a suppression of the polarised DPD with respect to its unpolarised counterpart. JHEP06(2017)083 We initialise the polarised F spl at µ y using the expression in (9.13) with the unpolarised 1 → 2 splitting functions replaced by their polarised counterparts (given in appendix B of [39]). These distributions are then evolved to the scale µ using the appropriate polarised double DGLAP equations (the required polarised splitting functions are collected in appendix A of [52]). The settings we used for the evolution code are the same as specified in section 9.2.1.
We consider the same scenario as in figure 12, i.e. the production of two systems with Q 1 = Q 2 = 80 GeV at √ s = 14 TeV, one with rapidity 0 and the other with rapidity Y , and now include polarised 1v1 contributions. In figure 19a, we reproduce figure 12a but include 1v1 luminosities for the polarised parton combinations a 1 a 2 b 1 b 2 = ∆u∆ū∆ū∆u + ∆ū∆u∆u∆ū and δuδūδūδu + δūδuδuδū. For brevity we refer to these combinations as ∆u∆ū and δuδū in the following, recalling that ∆q and δq indicate longitudinal and transverse quark polarisation, respectively. These polarised luminosities appear in the DPS cross section for double Z production, see section 3.1 of [53] (combinations like uū∆ū∆u also appear but are not shown here). To avoid cluttering the plots, we omit the unpolarised 2v1 luminosity, reminding the reader that it is of the same magnitude as 2v2. In figure 19b we repeat the exercise for the pure gluon channel, including the 1v1 luminosities for the polarised combinations a 1 a 2 b 1 b 2 = ∆g∆g∆g∆g and δg δg δg δg, referred to as ∆g∆g and δg δg, where ∆g and δg respectively denote longitudinal and linear gluon polarisation.
We see that the 1v1 luminosities for ∆u∆ū and δuδū essentially coincide with the one for uū at central rapidities, with ν variation bands of very similar size. These observations can be explained by the fact that at Y = 0 we initialise all three DPDs using the same expression up to a minus sign, with P g→qq = |P g→∆q∆q | = |P g→δq δq | at v = 1/2, and that evolution has a rather weak effect on all three DPDs for central rapidities. This close agreement extends out to higher values of Y for the uū and ∆u∆ū curves, owing to the fact that P g→qq = |P g→∆q∆q | for general v. By contrast, the δuδū curve clearly falls below JHEP06(2017)083  Figure 19. As in figure 12, but with two polarised 1v1 luminosity bands included in each case (and the unpolarised 2v1 band omitted). the other two for larger Y , because P g→qq > |P g→δq δq | for v = 1/2.
In figure 19b, the δg δg luminosity lies well below the one for ∆g∆g, which in turn lies somewhat below the one for gg. This is mostly driven by the differences in initialisation expressions for the DPDs: at v = 1/2 (relevant at Y = 0) the relevant splitting function satisfy for instance P g→gg : P g→∆g∆g : P g→δg δg = 9 : 7 : 1. Ignoring evolution effects, one would then expect the gg luminosity at Y = 0 to be roughly 1.5 times bigger than the one for ∆g∆g, which in turn should be roughly 50 times bigger than the one for δg δg. This expectation is quite close to the actual luminosity ratios, which are ∼ 4 and ∼ 60, respectively. The remaining difference is due to evolution effects, which increase gg rather considerably, increase ∆g∆g to a lesser extent, and hardly change δgδg.
Overall, we see that the polarised 1v1 luminosities can be of the same order of magnitude as the unpolarised ones, so that one must in general take into account all possible polarisation combinations in the DPS contribution, together with the SPS term and subtraction (where the subtraction will also contain both unpolarised and polarised contributions).

Production of two scalars
In this section, we study an explicit example to test our general argument that σ 1v1,pt is of similar size as the corresponding perturbative order of the SPS cross section. The process we investigate is an artificial one, chosen for ease of computation. We consider the hypothetical production of two identical massive scalar bosons φ, which couple to quarks with a Yukawa term cφqq, where c is a constant. For simplicity we take only one light quark flavour; including further light flavours just multiplies all following results (SPS and subtraction) by n f . We will not compare the subtraction term to the full SPS contribution, but rather to the gg initiated part of σ SPS , which is the piece that contains JHEP06(2017)083 the lowest order DPS-type gg → φφ box diagram. This piece is gauge invariant and can hence be meaningfully considered by itself.
In keeping with the notation of the paper, let us denote the mass of each produced boson φ by Q. Let Y be the rapidity of the diboson system in the pp centre-of-mass frame, s its squared invariant mass, and β = 1 − 4Q 2 /ŝ the velocity of one boson in the diboson centre-of-mass frame. The gg → φφ part of the SPS cross section can be written as where q max = βQ 1 − β 2 and the gg → φφ matrix element squared, |A gg→φφ | 2 , is averaged over the spin and colour of the incoming gluons. The momentum fractions x andx in the gluon distributions are obtained from the constraintsŝ = xxs and Y = 1 2 log(x/x). The matrix element squared for gg → φφ can be obtained from the gg → HH matrix element squared given in [54] by making the replacement (9.17) In particular, |A gg→φφ | 2 in (9.16) is related to the terms "gauge1" and "gauge2" in [54] according to The right hand side of (9.18) is given in [54] for general quark mass m q . We evaluate this expression numerically for very small m q , using analytic m q → 0 approximations where this is necessary to avoid numerical instabilities. Now we turn to the subtraction term. The two elementary qq → φ cross sections in this term contain the kinematic constraint 2πδ(x ixi s − Q 2 ), which fixes the momentum fractions x i andx i entering the two hard subprocesses at given x,x and β up to a two-fold ambiguity. The DPS subtraction term is given by where we sum over all possible colour, spin, and quark number interference/correlation possibilities. The index R on F denotes the colour channel (R = 1 for colour singlet, R = 8 for colour octet), and in the sum over a 1 a 2 b 1 b 2 we sum over both unpolarised (e.g. qqqq) and polarised (e.g. ∆q∆q∆q∆q) combinations. By {F → I} we denote the JHEP06(2017)083 same expression, but with the quark number diagonal distributions replaced by the quark number interference ones (see section 2 of [5]). The splitting kernels for the relevant spin combinations are [5] 1 P g→qq (v) = − 1 P g→∆q∆q (v) = 1 2 v 2 + (1 − v) 2 , 1 P g→δqδq (v) = − δ jj ′ v(1 − v) (9.20) and 8 P = − 1 P/ N 2 c − 1, where j, j ′ are indices for transverse quark polarisation. The term with interference DPDs I gives the same contribution as that with F , since the corresponding diagrams are simply related by reversing the direction of fermion flow in one of the two quark loops, and this does not change the expression for the diagram. The squared subprocess amplitudes, including an average over colour in the initial state, read where j, j ′ are the indices for transverse quark polarisation and where we have used the spin projection operators given in equation (2.90) of [5]. Inserting (9.20) and (9.21) into (9.19), we obtain dσ 1v1,pt dY dβ = xx g(x) g(x) 1 128π 1 Note that both (9.16) and (9.22) contain a product of gluon PDFs evaluated at the same x values, g(x) g(x). For the comparison, we can divide the common PDF factor out of the two expressions, in order to avoid having to use an explicit parameterisation. We also divide out various factors appearing in both expressions and compare the quantity Σ(β) = dσ dY dβ It is straightforward to show that Σ is a function of the variable β only. We plot Σ for SPS and the subtraction term in figure 20, where the curve for Σ 1v1,pt corresponds to the central choice ν = Q of the cutoff scale. We see that Σ 1v1,pt (β) is indeed generally of the same order of magnitude as Σ SPS (β). The agreement is perhaps surprisingly good, given that Σ 1v1,pt (β) involves integrating a low q approximation to the matrix element squared outside the region where the approximation is valid. The agreement gets worse towards the endpoints β → 0 and β → 1, which correspond to the high energy limit and the threshold limit, respectively. It is to be expected that the agreement is especially bad at these points, since in the subtraction term we effectively assume that the integration over the squared transverse momentum q 2 of the scalar particles goes from zero to values of order Q 2 . For both β → 0 and β → 1 this assumption becomes a poor one: near threshold the phase space in q 2 shrinks to zero, whilst in the high energy limit, q 2 can go up to sizeŝ ≫ Q 2 .

Summary
Consistently incorporating the perturbative splitting of one parton into two is a highly nontrivial problem for the theoretical description of double parton scattering. DPS graphs JHEP06(2017)083 in which such splittings occur in both protons (1v1 graphs) overlap with loop corrections to single parton scattering. Another type of graph, typically referred to as 2v1, in which one parton pair arises from a perturbative splitting, and the other pair is an "intrinsic" one already existing at the nonperturbative level, overlaps with twist-four contributions to the cross section. Finally there is an overlap between DPS contributions where a splitting occurs in both protons only in the amplitude or its conjugate, and SPS/DPS interference graphs.
We have presented a scheme to compute DPS and to consistently merge its contribution to the cross section with SPS and the other terms just mentioned. The scheme works in a similar manner for collinear and for TMD factorisation. Ultraviolet divergences that arise from perturbative splitting in a naive treatment of DPS are regulated by a cutoff function Φ(yν) in transverse position space. This avoids modification of the position space DPDs, which are defined via operator matrix elements in close analogy to single-parton distributions. In collinear factorisation, these DPDs hence evolve according to a homogeneous DGLAP equation, whilst their TMD counterparts satisfy a generalisation of the renormalisation group equations for single-parton TMDs. No modification of hard scattering cross sections computed for standard collinear or TMD factorisation is necessary in our scheme. Collins-Soper type equations describe the rapidity evolution of transversemomentum dependent DPDs and of collinear DPDs in colour non-singlet channels [21].
The problem of double counting between DPS and other contributions -notably between DPS and SPS -is solved by subtraction terms as specified in (4.16) and (4.22), which are obtained in a simple way from σ DPS by replacing the DPDs with their appropriate short-distance limits. This paves the way for using the scheme at higher orders in α s , with calculations being considerably simpler for the subtraction terms compared with the full hard scattering process at the corresponding order. With a suitable choice of starting conditions and scales, specified in section 6.1, the DPS part of the cross section correctly resums DGLAP logarithms that are not included in the fixed order twist-four contributions.

JHEP06(2017)083
Our scheme is naturally formulated with position-space DPDs F (x i , y), but it is possible to relate the Fourier transform of F (x i , y) Φ(yν) to DPDs F (x i , ∆) that are renormalised in transverse momentum space and satisfy an inhomogeneous DGLAP equation rather than a homogeneous one. This relation has the form of a perturbative matching equation, see (7.15), and is somewhat similar to the matching between PDFs defined in different schemes such as the MS and the DIS scheme. The momentum space representation also allows us to show that for the 2v1 contribution to DPS our scheme is equivalent to the ones in [12,13] and in [14,15] to leading logarithmic accuracy.
For collinear DPDs, one can make a model ansatz consisting of two terms which in the limit y ≪ 1/Λ respectively give the perturbative splitting and the intrinsic part of the distribution. With such an ansatz, the DPS cross section naturally splits into 1v1, 2v1 and 2v2 terms, where 2v2 refers to contributions in which the parton pairs from both protons are intrinsic. A crucial question is how large DPS is compared with SPS at the perturbative order where graphs contribute to both mechanisms. This is especially acute in collinear factorisation, where DPS is power suppressed with respect to SPS. Note that only in very few channels (notably pair production of electroweak gauge bosons) SPS calculations are available at the required order. We argue in section 6.2 that in our scheme the variation with ν of the 1v1 term in DPS provides an order-of magnitude estimate for the SPS contribution σ SPS (at the appropriate perturbative order), as it involves the same PDFs, overall coupling constants and kinematic region (small y, corresponding to large transverse momenta and virtualities of internal lines). An alternative estimate is provided by the double counting subtraction term σ 1v1,pt , which by construction is dominated by small y. For the hypothetical process of scalar boson pair production from two gluons, we have shown that the latter estimate works well within about a factor of two, provided that one stays away from the extreme kinematic regions where the velocity β of a boson in the boson pair rest frame is close to 0 or 1.
We constructed explicit (collinear) DPD input forms using the model ansatz just described, restricting ourselves to three quark flavours for simplicity and ease of implementation. These inputs were then numerically evolved to other scales using a code that implements the homogeneous double DGLAP equation. We used the resultant DPDs to compute so-called DPS luminosities (DPS cross sections omitting the process-dependent hard parts) and plotted these under various conditions. We observed that generically, the 1v1 contributions to the luminosity (both unpolarised and polarised) are comparable to or larger than the 2v1 and 2v2 contributions, with a large dependence on the cutoff ν. This demonstrates that, when including DPS in a cross section calculation, one must in general include σ SPS up to the order that contains DPS-like double box graphs, together with the associated subtraction term (with unpolarised and polarised partons). Otherwise, one would have an uncertainty on the overall cross section that is as large as, or larger than the DPS term itself. We also confirmed that the ν variation of the 1v1 DPS contribution is indeed comparable to the central value of the associated double counting subtraction term, so that either of them may be used as an estimate for the SPS contribution.
We identified several processes and scenarios where the ν variation of the 1v1 DPS luminosity is considerably smaller than the central value. As we argued above, one may