Soft factors for double parton scattering at NNLO

We show at NNLO that the soft factors for double parton scattering (DPS) for both integrated and unintegrated kinematics, can be presented entirely in the terms of the soft factor for single Drell-Yan process, i.e. the transverse momentum dependent (TMD) soft factor. Using the linearity of the logarithm of TMD soft factor in rapidity divergences, we decompose the DPS soft factor matrices into a product of matrices with rapidity divergences in given sectors, and thus, define individual double parton distributions at NNLO. The rapidity anomalous dimension matrices for double parton distributions are presented in the terms of TMD rapidity anomalous dimension. The analysis is done using the generating function approach to web diagrams. Significant part of the result is obtained from the symmetry properties of web diagrams without referring to explicit expressions or a particular rapidity regularization scheme. Additionally, we present NNLO expression for the web diagram generating function for Wilson lines with two light-like directions.


Introduction
The effects of double parton scattering (DPS), i.e. the scattering with two partons of a hadron participating in the hard subprocess, are usually expected to be small in comparison to a single parton scattering contribution. However, at very high energies the effect of multiple parton interactions increases and presents an important part of the total cross section, see e.g. [1][2][3]. It is known that DPS processes can form strong background for Higgs searches [4,5], as well as, be dominant channel for particular reactions, e.g. in the double Drell-Yan process p + p → W + + W + [6]. Therefore, the practical interest to DPS processes is constantly increasing.
From the theoretical side, the DPS processes are studied rather weakly. One of the reasons is the cumbersome kinematic structure of DPS. The double parton distributions (DPDs), the analogs of parton distribution functions for DPS, are functions of many variables: two momentum fractions and three transverse coordinates (or one transverse coordinate in the integrated case), say nothing of dependencies on two factorization scales. Additionally, DPDs have reach polarization and color structure, and even the leading order factorization formula for unpolarized and integrated double Drell-Yan involves more than dozen presumably independent DPDs [7,8]. Nonetheless, during recent years there was significant progress in the theoretical understanding of DPS processes, due to the formulation of appropriate factorization theorems [7][8][9].
Apart from increased number of various functions, the DPS factorization theorems resemble the factorization formula for transverse momentum dependent (TMD) processes, see e.g. [10]. It is not accidental since the dominant field modes are the same for TMD processes and DPS processes. This analogy grants the possibility to re-use the TMD experience during consideration of DPDs. For example, at NLO all the evolution properties for DPDs can be presented via corresponding evolution properties of TMD distributions [7,8,11].
In this work, we concentrate on the study of DPS soft factors, which are essential part of DPS factorization theorems. Soft factors represent the underlying interaction of soft gluons and contain the mixture of rapidity divergences related to both hadrons. This substructure should be decomposed into the parts with rapidity divergences belonging to a given hadron. Only after such decomposition a finite, i.e. meaningful, parton distributions can be defined. Naturally, the decomposition introduces the rapidity parameter. The anomalous dimension for the rapidity parameter scaling also can be deduced from the soft factor. Therefore, the study of the soft factor is an important part of the study of DPS factorization theorems.
At NLO the soft factors are nearly trivial objects. This order is given by single gluon exchange diagrams only. Therefore, DPS soft factors at NLO scatter into NLO TMD soft factors [7,8]. At NNLO many non-trivial aspects of perturbative expansion arise. The most important one is that simultaneous interaction of several Wilson lines becomes possible, and thus, one can expect highly interesting dynamics. However, the difficulty of consideration also grows. For example, the properties of TMD soft factor, although known for a long time, have been explicitly demonstrated at NNLO only recently [12].
The soft factors for DPS are rather involved objects composed of four half-infinite light-like cusps of Wilson lines positioned at four different points in the transverse plane and connected in all possible ways. Consideration of such an object within a usual diagrammatic is a serious calculation problem, mostly due to confusing combinatoric of the color flow. The structure of perturbative series is exceptionally simplified within the generating function approach for web diagrams, formulated in [13,14]. Within this approach, one should calculate the generating function, which is unique for a given geometry (in the case of TMD-like soft factors, the only important point is two light-like directions). Various matrix elements such as TMD soft factor, DPS soft factors, are obtained by a projection operation on the generating function. In this way, the usually difficult diagrammatic combinatoric is reduced to a couple of lines of simple algebraic manipulations.
One of the most attractive features of the generating function approach is an efficient organization of the expression. In particularly, it allows to avoid the calculation of whole sectors of diagrams, showing their equivalence with lower perturbative orders. As we demonstrate in this work, the consideration of the generating function at NNLO immediately shows the possibility to present any DPS soft factor in terms of TMD soft factors at this order. This fact is not trivial, since the NNLO expression contains products of TMD soft factor, but does not contain new functions. On the level of diagrams it implies that diagrams in particular combinations cancel each other, while in other combinations scatter into one-loop integrals. To find such combinations can be a tricky task, but they reveal automatically within the generating function approach.
In this paper we demonstrate that at NNLO DPS soft factors are given by a simple combination of TMD soft factors. Having at hands DPS soft factors at NNLO we study the structure of rapidity divergences and present the rapidity evolution equations at NNLO. It appears that some important results can be obtained without referring to expressions for diagrams. For example, the factorization of rapidity divergences for DPD soft factor at NNLO appears to be direct consequence of the rapidity factorization for TMD soft factor. We also perform the explicit calculation within the δ-regularization scheme [12] and confirm the results of the general analysis. The expression for the NNLO generating function for web diagrams presented here for the first time can be also used in other applications.
In the section 2 we review the derivation of factorization formula for double Drell-Yan process following articles [7-9, 15, 16]. This section is mostly needed to introduce the compact notation and necessary details about DPS soft factors. In the section 3.1 we give a short introduction to generating function approach for web diagrams. In sections 3.2 and 3.3 we discuss the details of evaluation of the generating function for DPS soft factors at NLO and NNLO, respectively. The particular form of projecting operators for DPS soft factors is given in section 3.4. In sections 3.5,3.6,3.7 we perform the projection operations and obtain the expression for DPS soft factors in terms of TMD soft factors. The origin of the simple structure of DPS soft factors is discussed in 3.6. In section 4 we discuss the influence of NNLO expression on the DPS factorization theorem. In particular, we show the factorization of rapidity divergences and define individual DPDs in sec.4.1. The rapidity evolution equations at NNLO are given in sec.4.2. The consideration is done in the most general case of unintegrated DPS. The important case of integrated DPS is obtained from these results and presented in the sec.5.
Technical details of the evaluation are collected in the set of appendices. In the appendix A we compare our notation with the notation used in [7] and [8]. The explicit expressions for diagrams, as well as, their analysis are given in appendix B. The expressions for basic loop integrals that participate in the generating function are collected in appendix C.

Factorization of double parton scattering
In this section we review some aspects of the double parton scattering factorization. The main aim of this section is to introduce notation and make connections with previous works. The consideration presented here is very superficial, and is an extraction from [7-9, 15, 16], to which we refer for the proofs. To get access to the leading order factorized cross-section we use the SCET II technique. The consideration is in many aspects similar to the consideration of Drell-Yan process at moderate transverse momentum [8,10,[17][18][19] (the so-called TMD factorization). Within this context, our main attention is devoted to the geometry of the double-parton scattering and to the color flow. Thus, we skip many important questions of DPS redirecting the reader to the literature [7-9, 15, 16].

Leading order factorization for double-Drell-Yan
The cross-section of double Drell-Yan process is given by the following matrix element [7,8] where dσ {µ} is a leptonic tensor and J µ ∼qγ µ q is the quark-to-vector boson current. Throughout the text an index enclosed in curly brackets denotes the set of indices of the same kind, e.g. here dσ {µ} = dσ µ1µ2µ3µ4 . As usual, we define two light-like vectors n andn along the largest components of P 1 and P 2 correspondingly, with n·n = 1. The vector decomposition reads v µ = n µ v − +n µ v + +v µ , where v are transverse components (P 1 = P 2 = 0) and v 2 > 0. The phase-space element dX denotes the complete phase-space of produced bosons, i.e.
in the reference frame), s = (P 1 + P 2 ) 2 is center-mass-energy and Y i = ln(q + i /q − i )/2 is the rapidity of produced boson. The large components of momenta are where Q is a generic large scale. One of the coordinates, say z 4 , can be set to zero due to translation invariance, but we keep it explicit for homogeneity of notation and for later convenience. Also in the following we often use the shorthand notation for set of arguments (b 1 , b 2 , b 3 , b 4 ) as b 1,2,3,4 (the order indices is important). In the following, although we introduce the notation convenient for our study, we try to be close to the notation and normalizations of [15].
Integrated double Drell-Yan process attracts even more practical interest. The integrated crosssection has the phase-space element dX = dx 1 dx 1 dx 2 dx 2 . It can be obtained from the unintegrated cross-section (2.1) by the integration over the transverse momenta q 1,2 . Consequently the expressions for the integrated anomalous dimensions, soft factors and over elements can be obtained from the unintegrated ones. In the following sections we consider only the general case of unintegrated kinematics. The expressions for the integrated case are collected in the section 5.
Following the SCET II factorization procedure we consider the quark field in the background gluon field, separating soft and collinear modes [17,[20][21][22], where a, a ′ are color indices, j is spinor index, W aa ′ n [z 1 , z 2 ] is a (soft) Wilson line from the point z 1 to z 2 , and field ξ n is the "large" component of quark field along vector n. The explicit definition of Wilson line is 3) The relation inverse to (2.2) is obtained by applying corresponding "large-component" projector ξ a n,j (z) = W aa ′ n [z, −∞]P n jj ′ q a ′ j ′ (z), (2.4) and similar for anti-quark andn components. Here P n = γ + γ − is the projector in n-direction. The Wilson line W has the same formal definition as W , but instead of soft gluons it consists of collinear ones. Substituting the field decomposition (2.2) into the matrix element (2.1) one obtains a large set of terms. The central point of the SCET approach is that at the leading order of factorization and in the absence of Glauber interaction (which has been proved in [9]), the field ξ does not interact with soft-gluons, and soft-gluon can be split up into separate matrix element. Then the cross-section is presented in the form vi,vi=n,n (2.5) where we extract the Lorentz structures from currents and absorb them into the tensor σ {ij} . The symbol Λ denotes a light-like cusp of half-infinite Wilson lines located at position z, (2.6) Figure 1. Graphical representation of contributions into the DPS factorization formula (2.7). The first terms in lines of (2.7) are presented. The out-going (incoming) to blobs fermion line represents quark (antiquark) field ξ (ξ). The double-lines represents the soft Wilson lines.
Note, that Λ ab vv (z) = δ ab . The dots in (2.5) denote the terms suppressed by powers of s [21,22]. The hard matching coefficients of the vector currents to SCET fields are hidden inside the function σ {ij} . The separation of the hard part introduces the renomalization scales µ's for each hard sub-process. In the following hard renormalization scales are taken equal to µ for brevity.
The further consideration is based on the following assumptions that are correct at the leading order of factorization in the region Λ 2 QCD ≪ q 2 1 ∼ q 2 2 ≪ s [7,8,10,17,21]: • Soft radiation does not resolve collinear scales, therefore, soft Wilson lines can be expanded at light-cone origin, Λ(z) → Λ(b), where b is the transverse component of z; • The "large" components of quark fields couple only to the hadron with corresponding momentum, i.e. ξ n couples to hadron P 1 , while ξn couples to hadron P 2 ; • The "large" component of the quark field does not resolve the scales in perpendicular direction, therefore, it can be expanded in that direction, i.e. ξ n (z) → ξ n (0, z − , b).
Using these assumptions one can compute the leading contribution to double-Drell-Yan process. The factorized expression contains various terms with usual parton distributions, double-partondistribution and combination that mix with each other within the operator-product expansion (see [23] for the leading order analysis). The separation of these terms from each other is an involved procedure (for theoretical development see [7,16,23,24]). In this article we are interested in the study of soft factors responsible only for multi-parton scattering. Therefore, we skip the discussion on the mixture between various matrix elements and consider only the DPS contribution. The DPS part of the cross-section corresponds to terms with simultaneous radiation of two distinct partons from the hadron. Applying leading order factorization restrictions to expression (2.5) and extracting the DPS contributions we obtain where • =↓=nn, and • =↑= nn (note, that order of arguments and indices of the function S is opposite to their order in the matrix element, i.e. graphical). The "arrow" notation is the visual representation of color-flow between light-cone infinities, i.e. if one writes all indices related to n∞ as down indices and all indices relatedn∞ as up indices, the arrows indicate the order of connection, see fig.1. The functions F are double parton distributions (DPDs) and given by expressions (2.10) The similarity sign implies possible normalization factor, which depends on the spinor structure. The distributionsF are obtained by changing components ± → ∓, n →n, and hadron states The visual representation of the terms in cross-section (2.7) is given in fig.1, it also illustrates the "arrow" notation for soft factors.
The following steps of classification consists in the Fiertz decomposition of spinor and color structures. As a result of this procedure one gets a large set of various DPDs with different polarization properties [7,8,25]. However, the details of Lorentz structure are inessential for the study of soft factor, while the color structure should be considered in details.
In fact, the factorization theorems (2.7,5.1) are not complete, in the sense that they consist of individually singular objects (DPDs and soft factors). They suffer from rapidity divergences, and are not entirely defined. Moreover, the soft factors mix the rapidity divergences related to different sectors of integration. The standard procedure implies that a soft factor can be presented as a product of factors with rapidity divergences from different momentum sectors. Then combining these factors with appropriate parton distributions one defines an "individual" parton distribution, which are finite and can be used in the phenomenology. Generally, it is unclear (although always implied) whenever it is possible or not to perform the rapidity-factorization procedure and define non-singular DPDs. In the section.4 we demonstrate that such procedure can be done at least at NNLO.

Color decomposition
The color structure of factorized expressions (2.7,5.1) is rather cumbersome. The notation introduced in (2.7-2.11) specially visualizes the color flow. The indices a and b denote the color adjusted to the antiquark and quark respectively. The subindex of color index designates the position of field in transverse plane (see fig.1). In this way, the SU (N c ) gauge transformation transforms DPD such that it is left with respect to indices a and right with respect to indices b. For example, where all matrices U are located at light-cone infinities. Consequently, the soft factor transforms in conjugated way by eight matrices U . In a non-singular gauge the transformation at light-cone infinites can be reduced to unity. In this way, DPDs and soft factors are gauge invariant objects independently. However, there is the global SU (N c ) rotation of quarks that still transform DPDs. Since global SU (N c ) is a symmetry of QCD, the DPD matrix elements select only the singlet contributions. There are two singlets in 3 ⊗3 ⊗ 3 ⊗ 3 that can be extracted as following here we use the normalization for singlet parts suggested in [7], which is different from the normalization used in [8]. The singlet parts of conjugated distributions are defined in similar manner. Substituting these expressions into (2.7,5.1) we obtain where the DPDs F and I are 2-component vectors F = {F 1 , F 8 }, and soft factors are 2 × 2matrices, which explicit form we present later. The notation (2.17) implies the presentation ofF as a "column", while F as a "row". This defines the order of matrix multiplication in the following sections.
Before we proceed to the definition of soft-factor matrices, let us discuss the symmetries of components. The obvious symmetry of a soft factor is the exchange of Λ's order under the sign of T-ordering. It implies S a1b1,a2b2,a3b3,a4b4 As the consequence of the Lorentz invariance, the directions n andn can be exchanged within the soft factor independently of the rest expression. Therefore, a soft factor is equal the soft factor with all arrows turned upside-down, where• =↓ (↑) for • =↑ (↓). Due to these symmetries the soft factors are related to each other. There are only two independent matrices S ↑↓↑↓ and S ↑↑↓↓ (as discussed later these two soft factors are also related by (2.34)). The rest soft factors are expressed as The soft factor matrix has four components, the Wilson lines of which are connected either by δ's, or by SU (N c ) generators. In turn, the product of generators t A t A can be expressed as products of δ's using Fiertz identities. For practical reasons, it is convenient to consider the products of Wilson lines contracted by δ's only. There are five independent structures that can appear A visual representation of these soft-factors is given in fig.2. The normalization is chosen such that at the leading perturbative order all soft factors are unity, The topology of the components is different, so the components S [1,2,3] form two Wilson loops, while S [4,5] are single Wilson loops. Further simplification of the structure can be made using features of the soft factor geometry special for the Drell-Yan process. The Wilson lines in the matrix element (2.8) are all positioned on the past light cone. Therefore, the distance between any two fields within (2.8) is space-like (or light-like if fields belong to the same Wilson line). It allows to rewrite the T-ordered product of Wilson lines as a usual product of Wilson lines, using the micro-causality relation. However, it is more convenient to organize Wilson lines as a single T-ordered product. We have Such presentation is distinctive feature of Drell-Yan kinematic, and is not possible for, say, double semi-inclusive deep-inelastic scattering (SIDIS). The representation (2.30) suggests higher symmetry of soft factor. Namely, the arguments can be freely exchanged preserving the topology of color-connection. Therefore, we need to conciser only soft-factors of different topology: a single Wilson-loop (we choose S [4] ), and a double Wilson-loop (we choose S [1] ). The rest are related to the chosen in the following way We also find that the soft factor S ↑↓↓↑ can be expressed via S ↑↓↑↓ as Therefore, we can consider only the case of S ↑↓↑↓ , while the results for other channels can be obtained by permuting vectors b i . Using the symmetries of soft factor (2.18-2.19) and the notation for independent components (2.24,2.27), we present the 2 × 2-matrix S ↑↓↑↓ in the form Here, for compactness we use the shorthand notation for the argument The rest of the soft factor matrices can be obtained via (2.20-2.23,2.34). At the leading order of perturbation theory soft factor matrices reduces to identity matrices

Generating function for web diagrams
The straightforward evaluation of functions S [i] requires a calculation of many diagrams, most of which are equivalent under permutation of parameters and change of color factors. Such a consideration would be very inefficient and contains many potential places for a mistake. A more effective approach is to evaluate the generating function for web diagrams, which is common for all soft factors, and project out the appropriate soft factor. The theoretical description of the approach can be found in [13,14]. In this section we describe only the basics of generating function approach needed for this particular calculation. The generating function approach is based on the well-known fact that the perturbative series for vacuum average of some operator sources is an exponent of the connected diagrams. This property immediately leads to exponentiation theorem for Wilson lines for Abelian gauge theories [13]. For a non-Abelian gauge theories one has an additional difficulty coming from the necessity to disentangle the color structure. The disentangling can be done in the general form [14]. In this way, one sees that significant part of diagrams that appear in the usual perturbation expansion (as well as, in the classical Wilson loop exponentiation diagrammatic [26,27]) are composed from the smaller-loop diagrams.
The power of the generating functions approach is that evaluated ones the generating function can be easily used to obtain the perturbative expression for any color topology. Therefore, the generating function that we present later can be used to obtain all DPS soft factors (2.24-2.28), as well as, TMD soft factor and soft factors for multi parton scattering with six, and more operators Λ. Moreover, the approach allows one to consider the exponentiated expression directly in the matrix form.
The starting point of the construction is to carry out the color structure of the Wilson line. The effective way to do so is to introduce the "scalar-reduction" of a Wilson line connecting points x and y [14] where t A is a gauge-group generators, θ A are c-number variables, and V is functional of gauge fields. The particular form of V needed for our calculation is given in (3.4), while the general form can be found in [13,14]. In this expression the matrix structure is carried by the first exponent in the product. The first exponent does not contain any fields and thus, does not participate in the function integration. The Wilson line flowing in the opposite direction can be presented in the form Within the considered task, we have a simple geometry of Wilson lines. All of them are straight (along n orn), and continue from b i to infinity (or in opposite direction). Let us enumerate these segments by number j (for DPS soft factor j = 1, .., 8). The j'th Wilson line can be presented in the form where v j =(n orn), the variable r j denoted the direction r j = 1(−1) if the color flow to (from) light-cone infinity. The operator V , which discribes a half infinite Wilson-line is given by [14] V In the following we omit the arguments of operator V for brevity. The functional V is an infinite series of path-ordered gauge field commutators. Here we truncate the series at g 2 order which is enough for NNLO calculation. The functionals V have many nice properties within the perturbation theory, especially, for light-like paths (see [14]). We use these properties for extra check of loopcalculation. Using the expression (3.3) we obtain the following expression for a generic soft factor where we explicitly write all color indices. The right-hand-side of expression (3.5) is a product of the color projector and the colorless matrix element.
In the case of a Wilson loop it is natural to rename the indices j such that they are ordered along the Wilson loop. The color projection operator in this case has especially simple form of P-ordered exponent of derivative operators, a discrete analog of Wilson line, where P-ordering is made with respect to index j. In the case of double Wilson loop the projector is product of two traces where subsets L1 and L2 of indices j are separately ordered along the first and the second loops. The matrix element on the right-hand-side of (3.5) can be presented as where W is the generating function for web diagrams. It is given by the sum of amplitudes where O is the connected part of the vacuum expectation value of the operator O. We have dropped the term linear in θ since it does not contribute in light-like kinematics. The dots denote the matrix elements with higher number of V . The first dropped contribution is ∼ V 4 , which is of order g 6 and hence NNNLO.
To obtain the expression for a soft factor we need to evaluate the correlator of two and three operators V at g 4 order. The different soft factors are obtained by application of different the color-projection operator.

Generating function at order g 2
In this section we discuss the evaluation of generating function at order g 2 in details. The calculation is almost trivial but we use it for the clarification of notation. Since V ∼ g, only the correlator of two operators V contributes at g 2 order. It is given by a single shown in fig.3.
The soft factor diagrams contains rapidity divergences. To regularize them we use the δregularization scheme as it is defined in [12]. Roughly speaking, the δ-regularization consists in the multiplication of every gluon field A(σn) within a Wilson line by the exponent factor e −δσ . So the gluon field is forced to zero at light-cone infinity. The positive infinitesimal parameters δ are different for n-andn− infinities, δ + and δ − respectively.
The Feynman rule for a single-gluon interaction with operator V A i is given by where k i = (v i · k), δ i = δ + (δ − ) for v i = n(n) and the momentum is incoming to vertex. The expression for diagram shown in fig.3 is where The following observations are useful also for many two-loop diagrams.
• The loop integrals are invariant under separate change of sign for transverse and light-like components. Therefore, the (overall) sign in the exponent(s) is irrelevant.
• The diagram is proportional to v ij , which is zero if both v i and v j along the same direction. Therefore, under the sign of integral we can set i = n and j =n without loss of generality. In this way, the loop-integral is independent on the vectors v i and v j , while this dependence is given solely by the prefactor v ij . Using these observations we combine terms of the generating function at g 2 into very simple form where a s = g 2 /(4π) 2 is QCD perturbative parameter. The explicit expression for K (0) is presented in (C.1). The transverse distances b ij within the perturbative expansion of soft factors are formally unrestricted. However, for the large values of b ij the logarithm contributions in (3.12) became large, and violate the convergence of the perturbative series. Therefore, practically the values of b ij should be restricted as a s (µ) ln(b 2 ij µ 2 ) < 1.

Generating function at order g 4
The diagrams contributing to the generating function at NNLO are presented in fig.4. In many aspects the calculation repeats the one loop calculation and details are presented in appendix B. To evaluate the diagrams one needs the Feynman rules for the radiation of two gluon from the effective vertex. It is where all momenta are incoming to vertex. Let us write the general form for the generating function at NNLO. It can be found from the symmetries of matrix elements. First, due to the global SU (N c ) symmetry a matrix ele- is invariant under permutation of operators V . And finally, due to the Lorentz invariance and power counting the functions W can have the scalar products v ij only as a prefactor. Combining together these observations the generating function at NNLO can be parametrized in the terms of two functions where We also confirm this form by direct calculation presented in appendix B. Note, that at NNNLO the structure of the generating function is reacher. It contains a term proportional to d ABC and various terms of the form f ABα f αCD and f ABα d αCD .
The evaluation of diagrams contributing to W [b] is nearly in one-to-one correspondence with the evaluation of similar diagrams in the case of TMD soft-factor made in [12]. The explicit expression for functions W within the δ-regularization are where base loop-integrals are given in appendix C. One can see that the complete NNLO expression for generating function contains only five basis integrals K (a) , I ′′ A , I C1 , I C2 , J and R, that are given in (C.1,C.3,C.4,C.5,C.6,C.9). These loop integrals can be compared and agree the loop integrals evaluated in [12].

Action of projection operator
As it was shown any soft factor with topology of a single Wilson loop can be obtained by the action of operator (3.6) on the generating function (3.14). The result of the action reads where sums over j are strictly ordered, i.e. j 1 < j 2 < ... < j n . Here and later, we use notation C F = (N 2 c − 1)/2N c and C A = N c where it is convenient. To derive the expression (3.19) we have used that W ij = W ji and that r 2 j = 1. The curly brackets on the indices of the third term denote the anti-symmetrization over permutation of indices (with 1/3! prefactor). One can see that the color factors which appear in the exponent corresponds to the color-connected parts of diagrams, in accordance of the exponentiation theorem for a Wilson loop [26,27].
The topology of the double Wilson-loop (composed of loops L1 and L2) is described by the projection operator (3.7). Applying (3.7) on the generating function (3.14) we obtain where S L1 and S L2 are soft factors evaluated on a single loop and defined in (3.19), and indices i and j belong to loops L1 and L2, respectively, and are strictly ordered along loops. The last two lines of (3.20) represent the interaction of Wilson loops. The leading order of between-loops interaction is given by double-gluon exchange, and, thus, is NNLO in coupling constant.

TMD soft factor
Before we proceed further it is instructive to calculate the TMD soft factor. Within the δregularization the TMD soft factor has been evaluated at NNLO in [12], and successfully used for description of TMD parton distribtions and TMD fragmentation functions at NNLO in [29,30].
Recently, TMD soft factor has been evaluated at NNLO in the rapidity regularization [31]. The TMD soft factor (in Drell-Yan kinematics) is given by the following matrix element In this case the index j runs from 1 to 4 and the parameters of the TMD soft factor are Evaluating the expression (3.19) with these parameters we obtain To present (3.24) in such a form, we have used that These relations follow from the expressions of loop integrals (B.11). Generally speaking, these symmetries are presented only at NNLO, but they are not crucial for further development and used only for visual simplification of the result. Substituting the explicit expressions for loop integrals into (3.24) we have checked that result (3.24) coincides with one presented in [12].

Soft factor S [4]
Let us consider the soft factor S [4] . We choose the enumeration of Wilson segments such that it starts at the point indicated in fig.2 as b 1 , and follows the arrows of color flow. Therefore, we have r j = (−1) j+1 , and Evaluating the formula (3.19) with these parameters we obtain an expression in terms of functions W (3.15,3.16). Considering this expression one can recognize the entries in the form of (3.24). It appears that the soft factor S [4] can be conveniently written via the TMD soft factor only. In terms of the function σ introduced in (3.24) the soft factor S [4] is remarkably simple To derive it we have used relations (3.25). An additional check is granted by the fact that setting two subsequent Λ's at the same point we obtain the TMD soft-factors. Namely, One can see that there is no direct three-lines correlations at this order. The multi-line correlation appears only as a product of pairwise interactions (the second line in (3.27)). In fact, this is a general feature of multi-particle soft factors and can be seen on the diagram level. Let us describe this effect at NNLO, although some part of statements can be generalized to arbitrary order.
At NNLO one has only two topologies of diagrams that connect three Wilson line. They are shown in fig.5. Every such diagram has a partner that is "reflected" upside-down (compare top and bottom diagrams in fig.5). The loop-integrals for pairs of such diagrams are the same, due to Lorentz invariance. However, the difference between diagrams can appear because of the different color connection and directions of Wilson lines. It is clear that the "reflected" diagrams necessary have opposite general sign due to the direction of Wilson lines, i.e. the combination r i r j r k changes sing under the "reflection". If three Wilson lines belong to different Λ's (see diagrams A and B in fig.5) then the "reflected" diagrams have the same color factor. Therefore, the diagrams that connects three different Λ's cancel each other. For the diagrams that connect two Λ's (see diagrams C in fig.5) the color factor also (together with the general sign) changes the sign under "reflection", and thus, these diagrams doubled in the final result. Therefore, the diagrams that connect three Λ's drop from the soft factors. As consequence DPS soft factor can be expressed in via TMD soft factors only.
The discussed cancellations are not transparent within a usual diagrammatic. The diagrams of type B in fig.5 are not presented in the usual diagrammatic. Instead there is a set of diagrams with two gluons connected to the Wilson line in different orders. The corresponding "reflected" diagrams have different loop-integrals, and cannot be directly compared. To perform comparison one should split these diagrams to symmetric (that part reduces to one-loop integrals) and anti-symmetric parts (that part is irreducible). The effective vertex V 2 represents the anti-symmetric contribution. This contribution cancels in the sum of mirrored diagrams. The symmetric combination is reducible and reveals after action of projection operator. Therefore, the generating function approach is effective tool for consideration of many Wilson lines configurations.

Soft factor S [1]
Let us consider the soft factor S [1] . We start the enumeration for the first loop from the point indicated in fig.2 as b 1 , and for the second loop from the point indicated as b 3 . Therefore, we have r j = (−1) j+1 , and r i = (−1) i+1 , and Evaluating expression (3.20) with these parameters, with accordance of discussion given in the previous section, we obtain expression in the terms of the TMD soft factor only. Using the function σ introduced in (3.24) we obtain An additional check is granted by the fact that setting two subsequent Λ's at the same point we obtain the TMD soft-factor. Namely, Substituting the expressions (3.27) and (3.33) for components S [1,4] in to the matrix (2.35), we obtain the explicit expression for S ↑↓↑↓ . The rest soft factors are obtained by the permutation of arguments as discussed in the section 2.2. We have checked that at NLO these expressions coincides with ones calculated in [7].

Recombination of rapidity divergences
As we have discussed in sec.2, the factorization formula (2.17) is incomplete, in the sense, that the collinear matrix elements F and soft-factor contain rapidity divergences. To complete the factorization formula and to construct a well-defined DPDs one has to disentangle rapidity divergences of n− andn-soft Wilson lines of the soft factor. These divergences recombine with the divergences arising in corresponding collinear matrix elements.
The DPS cross-section has a matrix structure (compare (4.2) and (2.17)). Therefore, the rapidity factorization procedure has to be done in the matrix form. In other words DPS soft factor should be factorized onto the product of matrices with the appropriate rapidity divergences. It implies the following expression Here we use a shorthand notation It is important to note, that there is a freedom in the decomposition (4.5) and in the definition of matrix (4.6). We have used that freedom to make diagonal terms pure functions of σ ± .

Evolution with rapidity parameter
The decoupling of rapidity divergences gives rise to the dependence on rapidity parameter ζ. The evolution equations with respect to the rapidity parameter are generally known as CSS equations (Collins-Soper-Sterman) [33]. Having at hands the soft factor one can extract the anomalous dimensions for the rapidity evolution (rAD). Let us remind this procedure in the case of TMD factorization.
The complete definition of a TMD contains several factors. The composition of these factors could be different within different formulations, compare e.g. [10,18,19,32]. However, the final result for "observables" i.e. anomalous dimensions and coefficient functions, is independent on a scheme. It has been recently confirmed at NNLO by direct calculations in different schemes [30,31,34]. Here we use the formulation based on δ-regularization [12,30]. Within δ-regularization a TMD reads where we drop all arguments except the argument related to rapidity parameters. The factor Z is the ultraviolet renormalization constant, which is dependent on ζ via cusp-logarithms, but δindependent. The factor S 1/2 is the part of the soft factor that comes from the cross section after the procedure of rapidity divergences separation (4.2,4.3,4.4). The soft factor in the denominator is the zero-bin subtractions [22,35], which are equal to the soft factor in the δ-regularization. Finally, F unsub is the TMD matrix element, that depends on δ only. In the products of these factors the regularization parameter δ cancels, and the TMD is rapidity-divergences-free. Then according to the definition of rAD one has The last term consists of terms that are divergent in ǫ. It is needed only for the cancellation of ǫ-divergences in the first term. Therefore, one can consider only the first term with all ǫ-divergences dropped, and thus, can obtain rAD solely from the soft factor. Substituting the expression for the TMD soft factor (4.3) we obtain where the mark f.p. denotes the selection of finite part in ǫ. Substituting the explicit expressions (3.17,3.18) into (3.24,4.3,4.10) we obtain the well-known result where Γ 1 = 1 and Γ 2 = C A 67 9 − π 2 3 − 10 9 N f are coefficients of the cusp anomalous dimension, β 0 = 11 3 C A − 2 3 N f is LO QCD β-function, and L(bµ) = ln(b 2 µ 2 e 2γE /4). This coefficient can be found in many papers, see e.g. the collection of formulae in [30]. Recently, rAD has been evaluated at NNNLO in [36].
In the case of DPS the TMD scheme can be used with minimal changes. We should only take care of the matrix structure of DPD. The analog of expression (4.7) in DPD case is (4.12) where S −1 is the inverse matrix of soft factor responsible for zero-bin subtractions, s is defined in (4.5), and Z is the ultraviolet renormalization matrix. The order of matrices is essential and follows from the scheme of divergence recombination [30]: the rapidity divergences are canceled prior to the ultraviolet renormalization and the zero bin subtraction are the part of DPD matrix element. 1 The order we set matrices in the expression (4.12) corresponds to our definition ofF as a "column". For the DPD F ,which is a "row", the whole composition should be transposed.
Defining the rAD matrix in the similar way we obtain . (4.14) To obtain the last equality we have commuted the matrix Z to the right and dropped singular in ǫ terms. It appears, that the expression for the rAD matrix can be found without referring to explicit expressions for loop integrals. Substituting the matrix σ − in the form (4.6) into the equation (4.14), we obtain the expression which consists entirely of function A in the form (4.10). Therefore, the rADs for DPDs can be expresses via the TMD rAD (4.11) where we drop the argument µ from the function D for brevity. The rADs for F qq , Fqq, as well as, for I qq and Iq q can be obtained by permutation of vectors b according to (2.20-2.23,2.34) There are some elementary checks of these expressions. First of all, these expressions do not contain δ-dependence, it cancel in the product of matrices (4.14). Second, the matrices D are symmetric matrices (although the matrix s is not), it implies that DPDs F and DPDsF evolve by the same equations, which in turn is the requirement of Lorentz invariance. The final result for the rADs at NNLO is the same (in pattern) as at the NLO [7]. In this way, the only difference from the result for NLO rapidity evolution given in [7], is that TMD rAD D should be taken at NNLO. This conclusion is very natural since any three Wilson line correlations disappear. However, this pattern does not hold at NNNLO where correlations of four Wilson lines appear. These correlators does not cancel in the sum of diagrams and give rise to a new "quadropole" contribution. 2

Expressions for integrated kinematics
In this section we collect the results for the integrated double-parton scattering. All presented results are obtained by from the expressions of the previous sections.
The integrated double-Drell-Yan cross-section is obtained from the integrated one (2.7) by the integration over momenta q 1,2 . In the regime Λ 2 QCD ≪ s the factorized expression is Contrary to the usual intuition for single Drell-Yan processes, the integrated soft factor matrices are not unity matrices. Correspondingly the integrated DPS contains the rapidity divergences, which are associated with vector y (5.1). The transverse vector y is not observable, and intrinsic for the factorization formula (5.1).
The soft factors in the integrated case can be obtained from the unintegrated ones by the relation (5.3). Using that σ ± (0) = 0 we find S ↑↓↑↓ (y) = S ↓↑↓↑ (y) = S ↑↑↓↓ (y) = S ↓↓↑↑ (y) = exp 0 0 0 CA CF σ(y) , The relations (5.4-5.5) must hold to high orders of perturbation theory and, in fact, represents the so-called Casimir scaling for a Wilson loop (i.e. expressions for Wilson loop of different representations differ only by the overall Casimir eigenvalue). Indeed, gluing together Λ in 88 element we obtain a single Λ in the adjoint representation, or TMD soft factor in the adjoint representation. However, it is possible that at four-or higher loop order this relation would be violated by the qubic Casimirs terms.
The anomalous dimension matrices for the integrated case are obtained by the same procedure as for integrated. Considering (4.14) with soft factors (5.4-5.5) we obtain D qq (y) = Dqq(y) = Dq q (y) = D qq (y) = 0 0 0 CA CF D(y) , One can see that the same expression can be obtained from (4.15-4.16) if we assume D(0) = 0.

Conclusion
We have considered the soft factor for the leading order factorization formula of the double parton scattering (DPS) in the perturbative regime (such that ln(b 2 Q 2 ) is not large, where b is any transverse separation within the soft factor). We have shown that at NNLO the soft factors are expressed entirely via the soft factor for transverse momentum dependent (TMD) factorization. This simplification happens due to the exact cancellation of the non-trivial part of three-Wilson lines interaction, while the trivial part of the three-Wilson line interaction can be presented via the pairwise interaction. Therefore, the DPS soft factor, that is generically a function of four transverse vectors, at NNLO reduces to products of simple functions of a single variable.
Using the fact that the logarithm of TMD soft factor is a linear function of rapidity divergence the DPS soft factor matrix can be split into a product of matrices that contain only the rapidity divergences in appropriate sectors (4.5). Adjusting the matrices to the unsubtracted double parton distributions (DPDs)(in other words, singular DPDs) we define DPDs that are free of rapidity divergences. It validates the DPD factorization theorem at NNLO. Using explicit expression for the soft factor matrices we extract the expressions for the rapidity anomalous dimension (rAD) matrices. At NNLO, the rAD matrices are expressed via rAD for TMD distributions (4.15,5.8). It appears that the expression for NNLO rAD follows the pattern of NLO expression with the only substitution of the two-loop anomalous dimension. Such simple pattern of DPS soft factor does not hold at higher orders. Starting from the NNNLO new functions appear. These functions represent the simultaneous interaction of four Wilson lines (quadrupole terms).
To express the DPS soft factor via TMD soft factors we used the generating function for web diagrams [13,14]. This is the first application of the method to a previously unknown object. This calculation shows high efficiency of the generating function approach. In fact, to obtain many results one does not need any explicit expression for diagrams but only the symmetry properties of Wilson lines which became transparent in the generating function approach. Using explicit expression for the generating function in the δ-regularization (presented in appendix B) we confirm the results of [12]. As a side result, we present the generating function for web diagrams at NNLO for Wilson lines with two light-like directions, which can be used to obtain any matrix element of Wilson lines with similar geometry, e.g. a soft factor for multi-parton scattering.
The presented calculation has been done in the Drell-Yan kinematics. In this case, all Wilson lines of the soft factor can be collected in the single T-ordered product. This trick reduces the number of diagrams and simplifies the calculation. In the case of Wilson lines pointing in different time directions (e.g. double semi-inclusive deep-inelastic scattering (SIDIS) kinematics), one should additionally consider the diagrams with real-gluon-exchange. However, we do not expect any special contribution from these diagrams at NNLO. The point is that the TMD soft factor is the same in both kinematics. This property is known as the universality of the TMD soft factor, see [12,38]. Therefore, we conclude that the DPD soft factor is also universal at least at NNLO.
We stress that relations between DPD soft factors and TMD soft factors are based on the geometry of the construction and color algebra only. Therefore, the results of the article are independent of the regularization scheme (although we imply some basic "good" properties of a regularization, such as non-violation of the exponentiation theorem).

Acknowledgments
The author gratefully acknowledges V.Braun and I.Scimemi for numerous stimulating discussions, and M.Diehl for useful comments. Also, the author acknowledges Erwin Schrödinger International Institute for Mathematics and Physics (ESI) for the support during the programme "Challenges and Concepts for Field Theory and Applications in the Era of LHC Run-2", which stimulates this research.

A Relation between normalizations
There are several sets of notations used for DPDs. We write the explicit relation of DPDs considered in this work to those introduced in [8] and in [7]. Since we did not consider the Lorentz structure we leave it hidden and discuss only the color decomposition and the order of transverse arguments.

B.3 Diagram 11
The evaluation of diagram 11 is straightforward, and in all aspects repeats the one done in [12]. Therefore, we present only the final expression where the last line is the counter term.

B.4 Diagram 211
The expression for diagram is The second term in the brackets is related to the first one by the replacement k ↔ l. Therefore, we need to consider only the following integral here k 2 = k 2 + i0. In this notation the diagram reads where expression for integral J is given in (C.6). Recalling that the generating function contains the sum of the diagrams with operator V 2 on different lines we obtain the contribution to W ijk in the form This representation is totally antisymmetric under the permutation of indices ijk. Taking into account the antisymmetry of prefactor θ A i θ B j θ C k if ABC we write it in the simple form, which is used in (3.16) At b ij = b ik = 0 the integral J turns to zero in accordance to symmetries of the generating function.

B.5 Diagram 111
The expression for the diagram 111 is Contrary to the previous diagrams within diagram 111 there are multiple possibilities to associate vectors v i , v j , and v k with n andn. For example, the term with v ij requires only v i = v j , while does not fix v k . The convenient way to consider this diagram is to write explicitly both possibilities for vector v k (i.e. v k = v i and v k = v j ) as v ij (v ik + v jk ). Then the numerator can be simplified Using this trick and change of variables (alike k → −k − l) in some terms we obtain the following expressions To obtain expression (B.10) we have also used the symmetries of the integral The expression (B.10) is totally anti-symmetric under the permutation of indices ijk. Using it we rewrite the diagram 111 in the form The expression for integral R is given in (C.9).

C Expression for loop integrals
It appears that one needs less number of base loop-integrals in comparison to the calculation presented in [12], but the integrals are of more general structure. In fact, the most part of integrals evaluated in [12] are particular cases of a few general integrals. The loop-integrals presented in the following paragraphs are all taken by the following strategy: • The integral over convenient light-cone components, say k + and l + are taken by residues closing the integration contour in lower-or upper-half of the complex plane. It also restricts the integrations over opposite light-cone components.
• The obtained propagators are expanded in Mellin-Barnes (MB) representation such that transverse components are separated from the light-cone ones.
• The integration over the transverse components is reduced to the scalar massless two-loop integrals in 2 − 2ǫ dimensions. Except the integral R these scalar integrals are product of Γ-functions only.
• The integrals over the rest light-cone components are table integrals of hypergeometric type. At this point of evaluation the expression for integral has form of single-or double-MB integral over product of Γ-function and possibly additional hypergeometric function. Note, that along all evaluation the regularization parameter ǫ > 0.
• The MB integration is taken by closing the contour. This integration is significantly simplified by the fact that we need only the leading small-δ asymptotic. It allows to consider only the poles in the vicinity of zero for many series of poles. The droped terms correspond to the power suppressed terms in the small-δ expansion.

C.2 Three point integrals
Three point integrals appear in the diagrams that connect three Wilson lines and contribute into W ijk . Generally, these integrals cannot be presented in the closed form.
The integral that appears in diagrams of group 211 is where L δ1 = ln(δ δ δB B B 1 e 2γE ), and where the integration contour passes between series of poles. The integrals J with exchanged variables are related to each other which can be checked explicitly using expression (C.6). If one of the variables zero, the expression can be simplified J(0, b) = 2 (4π) d δ δ δ −2ǫ Γ 4 (ǫ)Γ 2 (1 − ǫ) − B B B ǫ δ δ δ −ǫ F ǫ − 2B B B 2ǫ Γ 2 (−ǫ)Γ(−2ǫ)Γ(ǫ)Γ(1 + ǫ) . (C. 8) In this limit the integral has been evaluated in [12] with the same result, which can be seen by relation J(0, b) = −I ′ A . The integral that appears in diagram 111 is The R sing is a subintegral that contains the soft singularities and double-rapidity singularities and can be evaluated explicitly The regular part contains only single rapidity divergence, and has the following integral represen-tation R reg (b 1 , b 2 ) = 2 (4π) d is the derivative of hypergeometric function over the second index (it can be expressed as 3 F 2 -function). The integral over y is regular at ǫ → 0, and can be integrated order-by-order of ǫ-expansion. One can check that R(b 1 , b 2 ) = R(b 2 , b 1 ).
Within our calculation we do not need the complete expression for R but only its limiting cases, These cases can be related to the integrals evaluated in [12] as R(b, b) = I ′ C3 + I ′′ C1 , R(b, 0) = −I ′′