Double parton distributions out of bounds in colour space

We investigate the positivity of double parton distributions with a non-trivial dependence on the parton colour. It turns out that positivity is not preserved by leading-order evolution from lower to higher scales, in contrast to the case in which parton colour is summed over. We also study the positivity properties of the distributions at small distance between the two partons, where they can be computed in terms of perturbative splitting kernels and ordinary parton densities.


Introduction
Parton distribution functions and related quantities are crucial ingredients for computing hadronic cross sections at high energies, and they are the main quantities that describe the structure of hadrons at the level of quarks and gluons. It is hence important to know and understand their general properties. One of these properties is positivity. For ordinary parton distributions (PDFs) this is just the statement f a (x) ≥ 0 if the parton a and the hadron are unpolarised. In the polarised case, this generalises to a set of inequalities, namely the well-known Soffer bounds on polarised distributions [1]. Corresponding bounds have been formulated for transverse-momentum dependent distributions (TMDs) in [2], for impact parameter distributions in [3], and for double parton distributions (DPDs) in [4]. The latter appear in the description of double parton scattering and contain a wealth of information about correlations between two partons in a hadron; for a recent review we refer to the monograph [5]. DPDs have a non-trivial dependence not only on the polarisation of the partons, but also on their colour, and corresponding positivity bounds have been derived in [6]. Positivity bounds can be of considerable practical value. They may be used as constraints in fits of PDFs, and in the context of spin physics, ansätze that saturate certain bounds are often used to estimate the maximal allowed size of spin asymmetries. A corresponding strategy for DPDs in spin or colour space appears all the more attractive because our current knowledge of these distributions is very incomplete. Whether positivity bounds on parton distributions actually hold turns out to be a non-trivial question. At leading order (LO) accuracy, the positivity of PDFs can quite directly be deduced from the positivity of cross sections, whilst the situation is more involved at next-to-leading order (NLO) in the strong coupling [7,8]. When derived along such lines, positivity holds for renormalisation scales µ that are high enough for the approximations of leading-twist dominance and of the perturbative expansion to be valid. To formulate a corresponding approach for DPDs would be complicated due to the large number of involved degrees of freedom (two pairs of partons in each colliding hadron and two hard-scattering subprocesses), and we shall not pursue such an avenue here.
Intuitively, the positivity bounds on parton distributions are a consequence of their partonmodel interpretation as number densities or linear combinations of number densities. At a more formal level, one may use light-cone quantisation and write appropriate linear combinations of distributions as squared operator matrix elements that are summed over unobserved degrees of freedom. Equivalently, one can represent the distributions in terms of light-cone wave functions. This has for instance been done for ordinary PDFs in [9] and for impact parameter distributions in [10]. A corresponding representation holds for DPDs (see [11] for the momentum space version, which can readily be adapted to the case of definite transverse parton position). A limitation of this approach is that it does not account for the renormalisation of ultraviolet divergences in the matrix elements (at least not in customary schemes such as MS) nor for subtleties related with Wilson lines or the definition of light-cone gauge at infinite light-cone distances. It has long been realised that ultraviolet subtractions can in principle invalidate the positivity of distributions. A detailed discussion and examples can be found in the recent paper [12]. An important result is that LO DGLAP evolution to higher scales conserves the positivity of PDFs, both in the unpolarised sector [13,14] and in the polarised one [15,16]. This means that if PDFs satisfy the positivity bounds at a certain scale µ, these bounds remain valid at higher scales when the PDFs are evolved at leading order. Discussions for NLO evolution can be found in [17,18]. Conversely, experience shows that PDFs eventually turn negative when evolved down to very low scales. Examples for this can for instance be found in [19]. In [4] it was shown that positivity of spin dependent but colour summed DPDs is conserved by LO DGLAP evolution to higher scales. The first goal of the present paper is to investigate whether the same holds for the bounds derived in [6] for DPDs with non-trivial colour dependence. In addition to DGLAP evolution, we will also consider Collins-Soper evolution in a rapidity variable, which appears when the parton colours are not summed over. Following [6] we will limit ourselves to unpolarised partons throughout this work. For small transverse distances y between the two partons, DPDs can be computed in terms of ordinary parton densities and kernels for the perturbative splitting of one parton into the two observed partons and (at higher orders) additional unobserved ones. Using the results of the recent two-loop calculation in [20], we can investigate to which extent positivity of DPDs in colour space is realised in the small y limit. This is the second goal of our work. This paper is organised as follows. In section 2, we specify the two parametrisations for the colour structure of DPDs used in this work, specify the property of positivity, and discuss the perturbative splitting mechanism for DPDs at LO accuracy. In section 3, we analyse whether Collins-Soper evolution from smaller to larger scales conserves positivity of DPDs, and in section 4 we perform a corresponding analysis for LO DGLAP evolution. The perturbative splitting mechanism at NLO accuracy is investigated in section 5. Our results are summarised in section 6, and some technical formulae are collected in an appendix.

Colour structure of DPDs
In this section, we discuss the general colour structure of quark and gluon DPDs and state the hypothesis of positivity in colour space. Throughout this work, we consider distributions for unpolarised partons. As illustrated in figure 1, the colour structure of a DPD can be described in terms of four colour indices, one for each parton field in its definition. We generically write F r 1 r 1 r 2 r 2 a 1 a 2 , where a 1 and a 2 denote the two parton flavours. The colour indices are in the fundamental or adjoint representation as appropriate, with r 1 and r 2 referring to the partons in the amplitude of the scattering process and r 1 and r 2 to the partons in the conjugate amplitude. As described in [21], the definition of a DPD involves the hadronic matrix element of two twist-two operators, Figure 1: Assignment of colour labels for a quark-antiquark distribution (a) and a quarkgluon distribution (b). The dashed vertical line indicates the final state cut of the scattering process in which the distributions appear.
as well as a soft factor given as the matrix element of Wilson line operators in the vacuum. Both matrix elements contain ultraviolet divergences that need to be renormalised. One can take different renormalisation scales µ 1 and µ 2 for the two partons, and the dependence of the DPD on these scales is given by DGLAP equations, which will be discussed in section 4. The soft factor removes rapidity divergences in the hadronic matrix element, in a similar way as in the definition of transverse-momentum dependent distributions [22]. This leads to a dependence of the DPD on a rapidity scale ζ p , which is described by a Collins-Soper equation as discussed in section 3.
The s and t channel colour bases. The four colour indices of a DPD must be coupled to an overall colour singlet. This can be achieved by first coupling the colour of two parton pairs to an irreducible representation and then coupling these two representations to an overall singlet. Depending on the choice of parton pairs, we consider two bases for the colour coupling. In the s channel basis, we pair the partons in the amplitude and in the conjugate amplitude. The projection on irreducible representations can then be written as for (a 1 a 2 ) = (qq), (qq) with colour indices r 1 , r 2 , r 1 , r 2 in the fundamental or adjoint representation as appropriate.
The multiplicity m(R) of the representation R is m(1) = 1, (1) follows the choice made in [6] and corresponds to eq. (6) in [23]. We denote the conjugate of a representation R by R, where it is understood that some representations like the singlet and the octet are their own conjugate. The matrix P R R in (1) projects the colour of the parton pair a 1 a 2 on the representation R in the amplitude and on the representation R in the conjugate amplitude. Its explicit form is given in the appendix. Quark-antiquark distributions F RR qq are defined as in (1), but with transposed indices r 2 r 1 in P R R , which ensures that covariant indices are always contracted with contravariant ones. Likewise, the definition of F RR qq has transposed indices r 2 r 1 . In the t channel basis, we pair the partons with momentum fractions x 1 and x 2 . Following [21,24], we write  Table 1: Combinations of colour representations in the s channel distributions F RR a 1 a 2 and the t channel distributions R 1 R 2 F a 1 a 2 . Two adjoint indices can couple to a symmetric (S) or an antisymmetric (A) octet. For two-gluon distributions, the colour combinations are identical in the s and t channels. Interchanging a 1 ↔ a 2 implies interchanging R 1 ↔ R 2 while keeping RR unchanged. (2) projects the colour indices r i r i of parton a i on the representation R i for i = 1, 2. For distributions with antiquarks, one needs to transpose the corresponding colour indices in P R 1 R 2 . The correct ordering is hence r 1 r 1 if a 1 =q and r 2 r 2 if a 2 =q. Combining the definition (1) with the completeness relation (63) and the explicit form of the singlet projector P 11 , one readily finds that for all parton combinations, where the sum runs over all relevant colour representations R. Throughout this work, we fix N = 3 for the number of colours. The colour factors used in later results thus have the values C F = 4/3, C A = 3, and T F = 1/2. The accessible colour representations for different parton combinations are given in table 1. We note that a quark and a gluon can couple to 6 (see e.g. table 24 in [25]) rather than to 6 as stated in equation (8c) of [23]. As shown in [21], the t channel distributions R 1 R 2 F are real valued, except for the decuplet sector in the pure gluon case, in which case one has ( 1010 F gg ) * = 1010 F gg . In the s channel basis, this translates into all distributions F RR being real, except for the mixed octet combinations, where one finds (F AS gg ) * = F SA gg .
Density interpretation and positivity. The parton model interpretation of DPDs can be obtained in the same way as for single parton distributions by expressing the field operators in terms of creation and annihilation operators in light-cone quantisation, neglecting all complications from Wilson lines and from renormalisation. Details can for instance be found in [22]. The t channel distributions are normalised such that 11 F a 1 a 2 (x 1 , x 2 , y) is the probability density for finding partons a 1 and a 2 with momentum fractions x 1 and x 2 at a transverse distance y from each other, with the density measure being dx 1 dx 2 d 2 y. The colours and polarisations of both partons are summed over in 11 F a 1 a 2 . Correspondingly, the s channel distribution F RR a 1 a 2 is the probability density for finding the parton pair in one of the m(R) states of the colour representation R. This provides an intuitive interpretation of the relation (3).
The positivity property for DPDs in full colour space is then the statement that which of course implies the weaker condition 11 F a 1 a 2 ≥ 0. Note that we define "positivity" as including the value zero. Note that in the pure gluon channel, the distributions in the s channel basis include the cases F AS gg and F SA gg , which correspond not to densities but to interference terms in colour space (and which may be complex valued as noted above). Given the large number of accessible colour channels in that case, we will not consider two-gluon DPDs in the remainder of this work, concentrating on the pure quark-antiquark sector and on mixed quark-gluon or antiquarkgluon distributions.
Basis transformations. Whilst the s channel basis is natural for considering positivity, the evolution of DPDs in the renormalisation and rapidity scales is much simpler in the t channel basis. We hence need the explicit transformations between the two representations. In the pure quark sector, the transformations between s and t channel bases read The transformation matrix for two antiquarks is Mqq = M qq . For gq and gq distributions, one has Perturbative splitting at leading order. If the transverse distance y between the two partons is small, DPDs can be computed in terms of a perturbative splitting process and ordinary PDFs. Example graphs for the perturbative splitting are shown in figure 2. This mechanism is interesting in our context because it generates a non-trivial colour dependence. Let us take a closer look at the splitting process at one-loop order, postponing the discussion of two-loop accuracy to section 5. The splitting formula at leading order reads [21,24] where f a 0 (x, µ) is the PDF for parton a 0 and we defined a s (µ) = α s (µ) 2π (10) and Note that (9) is an approximation for small y and receives corrections suppressed by a power of a s or y 2 Λ 2 , where Λ is a hadronic scale. The splitting kernels 11 V (1) a 1 a 2 ,a 0 (z) for the colour singlet channel are equal to the usual LO DGLAP splitting functions without the distributional parts (plus prescription and delta function) at z = 1. One therefore has 11 V (1) a 1 a 2 ,a 0 (z) > 0. This implies that 11 F a 1 a 2 > 0 at the scale µ where the LO splitting formula (9) is evaluated, provided of course that the PDFs are positive. To avoid large a 2 s corrections, one should take µ ∼ 1/y. The LO splitting kernels for other colour channels are proportional to 11 V (1) a 1 a 2 ,a 0 (z), and it is easy to transform them to the s channel colour basis. The result is where R 0 is the colour representation of the initial parton a 0 of the splitting process (with R 0 = A for g → gg). The parton pair a 1 a 2 must be in the representation R 0 because the LO splitting graphs are disconnected between the amplitude and conjugate amplitude (see figure 2(a)). The result (12) then follows from the relation (3). With 11 F a 1 a 2 > 0 one hence finds positivity of DPDs in colour space when taking the LO splitting approximation. We will see in section 5 whether this still holds at NLO.

Collins-Soper evolution
In this section, we investigate how Collins-Soper evolution affects the positivity of DPDs.
Collins-Soper evolution of DPDs does not mix different colour representations in the t channel basis, where one has Here we displayed all arguments of the functions. The Collins-Soper kernel R 1 J depends only on the multiplicity of R 1 (which is equal to the multiplicity of R 2 ) but not on the parton types. Note that colour singlet distributions in the t channel are ζ p independent, i.e. 1 J = 0. For all parton combinations except gg, the only non-trivial kernel needed is hence the one for the colour octet. Remarkably, this kernel satisfies the exact relation [26] 8 where K g (y, µ) is the Collins-Soper kernel for the evolution of single-gluon TMDs. Let us discuss the sign of the Collins-Soper kernel, which will be important in the following. The renormalisation group equation for the Collins-Soper kernel is solved by with a positive anomalous dimension that is proportional to the cusp anomalous dimension for adjoint Wilson lines: At given y one can hence always achieve a negative 8 J by taking the scales µ 1 and µ 2 sufficiently high.
To make a more specific statement, we first consider small distances y, where one can compute the kernel in perturbation theory and obtains where b 0 = 2e −γ ≈ 1.12 and γ is the Euler-Mascheroni constant. Bearing in mind that there are higher-order terms in (17), we see that the transition from positive to negative 8 J happens at µ around b 0 /y, as long as y remains in the perturbative regime. Not much is known about 8 J or K g for y in the nonperturbative domain. The situation is different for the Collins-Soper kernel K q for quark TMDs. Several phenomenological extractions find that K q (y, µ) is negative for large y, see e.g. [27] (figure 6), [28], and [29] (figure 23). Furthermore a number of lattice determinations, covering a distance range between about 0.1 fm and 0.8 fm, find that K q (y, µ) < 0 at µ = 2 GeV, see figure 7 in [30], figure 5 in [31], and figure 8 in [32]. We find it plausible to assume a qualitatively similar behaviour of K g (y, µ) and K q (y, µ) as functions of y (at perturbatively small y, one actually has K g /K q ≈ C A /C F up to corrections of order α 4 s , see footnote 10 in [21]). Under this assumption, we conclude that 8 J(y, µ 1 , µ 2 ) < 0 for scales µ 1 and µ 2 sufficiently larger than max(b 0 /y, 2 GeV).
Collins-Soper evolution in the s channel. In the s channel, different colour representations mix under Collins-Soper evolution. Starting with the qq channel and using the basis transform given in the previous section, we get Restoring all arguments, we find that the evolution equation is solved by with the matrix exponential where we abbreviate The Collins-Soper equation for the other parton combinations we consider has the same form as (18) with appropriate changes in the colour labels and matrices. In the quark-antiquark case, one hasĴ and for gq distributions we find where α is always given by (22). We furthermore getĴqq =Ĵ qq andĴ gq =Ĵ gq and corresponding equalities for the evolution matrices U a 1 a 2 . The evolution equations for distributions F q 1 q 2 , F q 1q2 , etc. with unequal flavours involve the same matrices as their counterparts for equal flavours.
We see that for all parton combinations a 1 a 2 except gg (which we do not consider) all elements of the evolution matrix U a 1 a 2 (α) are positive for α > 0. This is the case for forward evolution (ζ p > ζ 0 ), provided that 8 J < 0, which is the case when µ 1 and µ 2 are sufficiently large. Under this condition, Collins-Soper evolution to higher scales thus preserves positivity. With the notation the Collins-Soper equation in the s channel basis reads for all parton combinations considered here. Using that (Ĵ a 1 a 2 ) 2 =Ĵ a 1 a 2 , we can write its solution in the form Writing the relation (3) at rapidity scale ζ p and using that 11 We can now discuss the behaviour of Collins-Soper evolution for large negative α, which is relevant when evolving backward with 8 J < 0, and when evolving forward at scales µ 1 and µ 2 that are so low that 8 J > 0. In the regime where the factor e −α in (27) is much larger than 1, the condition (28) implies that the evolved distributions F RR a 1 a 2 (ζ p ) must be positive in some colour channels and negative in others. An exception to this statement is the case where the initial conditions satisfy (Ĵ a 1 a 2 F a 1 a 2 (ζ 0 )) RR = 0 for all R, so that all distributions are independent of ζ p . In the t channel basis, this is tantamount to all distributions other than 11 F a 1 a 2 being zero. Apart from this special case, evolution to large negative α always leads to a violation of positivity.

DGLAP evolution
In this section, we investigate how leading-order DGLAP evolution affects the positivity of DPDs. We limit ourselves to the evolution of two-quark and quark-antiquark distributions, which are the simplest cases as far as mixing and the number of colour channels are concerned. We first discuss evolution in µ 1 at fixed µ 2 and ζ p and then evolution in all three scales simultaneously.

Evolution in the scale of one parton
Let us consider evolution in the renormalisation scale of one parton, which we take to be the first one without loss of generality. For the parton combinations of interest, the LO evolution equations in the t channel basis read where the second parton a is a quark or an antiquark, and where R = 1 if R 1 = 1 and In the second equation we have made use of the charge conjugation relations betweenqg and qg splitting kernels; this gives a sign factor ε 2 (A) = −1 for the antisymmetric octet. The evolution equations involve the Mellin convolution whose lower integration boundary reflects the support region of DPDs in the momentum fractions: F (z 1 , z 2 , . . .) is zero for z 1 + z 2 > 1. The evolution kernels can be written as where b = q, g, γ q (µ) = 3C F a s (µ) + O(a 2 s ), 1 γ J (µ) = 0, and 8 γ J (µ) is given in (16). The z dependent part of (31) involves the familiar splitting functions and the colour factors Whilst our analysis is limited to the LO approximation of the splitting kernels, our results do not depend on whether one uses the LO or the NLO approximation for the anomalous dimension 8 γ J (µ), which is associated with Sudakov double logarithms. Taking such anomalous dimensions at NLO corresponds to next-to-leading logarithmic (NLL) approximation (see e.g. table 1 in [33] for single hard scattering and section 6.6 in [21] for DPS). If one takes 8 γ J (µ) at two-loop accuracy, one may also want to use the two-loop rather than the one-loop running of α s (µ). Our arguments in the present work do not depend on that choice. We will shortly need to know the sign of the Mellin convolutions (30).
On the other hand, the plus-prescription for P qq involves a negative term proportional to F (x 1 , x 2 , . . .): It can therefore have any sign, even if F (z 1 , z 2 , . . .) ≥ 0 for all z 1 , z 2 . To illustrate this, let us consider the DPDs computed with the LO splitting formula (9). We take the PDFs of the CT14lo PDF set [34], using the LHAPDF interface [35] via ManeParse [36]. We evaluate the DPD at µ = b 0 /y = 10 GeV and verify that the PDFs are positive at that scale. For the strong coupling, we use the value α s (10 GeV) = 0.178 provided by the PDF set. As shown in figure 3, the convolution of P qq with 11 F qq so obtained is indeed negative in a large region of the momentum fractions. The evolution equations in the s channel basis are derived along the same lines as the Collins-Soper equation in (18). If the first parton is a quark, we get where the function arguments of P and F are as those of P and F in (29). The colour mixing matrices read whilstĴ qq andĴ qq are given in (19) and (23), respectively. With we find that the evolution equations for Fqq and Fq q are respectively obtained from (35) and (36) by interchanging q ↔q in the DPDs and swapping their representation labels, i.e. F 33 qq → F 33 qq , F 33 gq → F 33 gq , etc. The splitting kernels, anomalous dimensions, and colour mixing matrices remain the same.  The convolution x 1 x 2 P (z, µ) ⊗ x 1 11 F qq (z, x 2 , y, µ, µ) with P (z, µ) = a s (µ) P qq (z) and 11 F qq computed from the LO splitting formula (9) at µ = b 0 /y = 10 GeV. A weighting factor x 1 x 2 is included to keep details visible at higher momentum fractions.
The evolution equations have the same form for distributions with unequal flavours, i.e. one may replace F qq → F q 1 q 2 and F gq → F gq 2 in (35), or F qq → F q 1q2 and F gq → F gq 2 in (36), whilst keeping the splitting kernels and colour mixing matrices unchanged. Before analysing the effect of the evolution equations (35) and (36) on positivity, let us recall the situation for colour singlet distributions in the t channel basis. The LO evolution equation where a is a quark or an antiquark. If 11 F qa and 11 F ga are non-negative for all momentum fractions, the terms involving P qg or γ q in (40) are non-negative as well and hence conserve positivity. As we have shown in figure 3, the first term in (40) can be negative. However, the part that is responsible for a negative sign is proportional to 11 F qa (x 1 , x 2 , . . .) itself, as is easily seen in (34). This negative contribution hence decreases in magnitude when 11 F qa (x 1 , x 2 , . . .) approaches zero from above, and closer inspection shows that it cannot lead to a violation of positivity. This is shown in more detail in appendix B of [4]. Overall, positivity is hence conserved for LO evolution of 11 F qa to higher scales, and the same can be shown for all other parton combinations. We now analyse the sign of the different terms on the r.h.s. of (35) and (36) under the assumption that the s channel distributions on the r.h.s. are non-negative for all momentum fractions. For brevity, we write D 1 F = ∂F/∂ log µ 2 1 .

The scale variation
A negative contribution from the term with a R cannot lead to a violation of positivity, as just discussed. However, the contribution from the term with b R can remain large and negative even if F RR qa approaches zero. This term can therefore lead to a zero crossing of the distribution as one evolves to higher scales.
2. The matrices P qg,a have no negative elements, so that the terms with P qg are all nonnegative.
3. The terms with γ q are all non-negative.
4. The contributions with 8 γ J to D 1 F 33 qq and D 1 F 66 qq have opposite sign, as well as those to D 1 F 11 qq and D 1 F 88 qq . This follows from the form ofĴ qq andĴ qq in (19) and (23). Whether these contributions can violate positivity depends on the sign of log µ 2 1 /(x 2 1 ζ p ).
It follows that LO DGLAP evolution of DPDs is not guaranteed to conserve their positivity in colour space, in contrast to the colour singlet distributions 11 F a 1 a 2 . This is one of our main results. As a numerical illustration, let us take the initial conditions provided by the perturbative splitting mechanism at LO. According to (12), we have F 11 qq = F 66 gq = F 1515 gq = 0 at the scale µ where the splitting formula is evaluated. Inserting this in the evolution equation (36) and taking the LO approximation (16) of 8 γ J , we obtain at the point µ 1 = µ 2 = µ. We evaluate the r.h.s. numerically with the same settings as in figure 3 and the choice which is natural for the perturbative splitting mechanism [20]. With this choice, the term with ζ p is zero at x 1 = x 2 . At that point, a negative value of (41) must hence be due to the term with P qq . We see in figure 4 that there are regions in x 1 and x 2 for which D 1 F 11 qq < 0 at the scale µ. With F 11 qq = 0 at that scale, positivity is thus explicitly violated by evolution to a higher scale for the first parton.

Simultaneous evolution in all scales
When computing double parton scattering cross sections, the choice of scales µ 1 , µ 2 and ζ p is driven by the kinematics of the process. In particular, taking µ 1 = µ 2 is natural for processes with two hard scales of very different size. On the other hand, the simplest setting for the physical interpretation of a DPD is with all relevant scales set equal. We therefore investigate the evolution of DPDs with a common renormalisation scale µ = µ 1 = µ 2 for both partons and the rapidity scale given by ζ p = µ 2 /(x 1 x 2 ) as in (42). Let us briefly comment on the factor x 1 x 2 in (42). As explained in [20,21], the definition of ζ p involves the rapidity regulator and the plus-momentum of the target proton. By contrast, x 1 x 2 ζ p refers to the regulator and the plus-momenta of the two extracted partons, and in this sense is more closely related to the scales µ 1 and µ 2 that refer to the renormalisation of the operators associated with the two partons. This motivates our choice (42), along with the fact that x 1 x 2 ζ p is the combination of variables appearing at higher orders in the perturbative splitting formula for DPDs, see (44) and (45). Combining the DGLAP equations for the first and second parton with the Collins-Soper equation (18) yields With ζ p chosen as in (42), the term with 8 γ J in (35) has cancelled against the corresponding term from the evolution equation in µ 2 . The sign of the r.h.s. of (43) can be analysed along the same lines as in the previous subsection: 1. The terms with P qq can lead to a violation of positivity because of the negative part in the plus-distribution and of the positive off-diagonal entries in the matrix P qq,q .  Figure 4: The right-hand side of (41) with qq = uū, evaluated at µ = b 0 /y = 10 GeV and ζ p = µ 2 /(x 1 x 2 ) with LO splitting DPDs. The individual terms going with P qq , P qg , and γ J are shown as well. The PDFs used in the splitting formula are specified in the text below (34). A weighting factor x 1 x 2 is included to keep details visible at higher momentum fractions. colour channel gives a positive contribution. This is consistent with our findings for Collins-Soper evolution in section 3.
The leading order expression (17) of 8 J(y, µ, µ) contains an explicit logarithm log(µ 2 y 2 ). In the leading double logarithmic approximation, one therefore has to keep only the last term in the evolution equation (43), which then reduces to the LO Collins-Soper equation with µ 2 = x 1 x 2 ζ p . The evolution equation for F qq has the same form as (43) and involves the colour mixing matrices P qq,q , P qg,q , andĴ qq . Its discussion proceeds in full analogy. In summary, we find that evolution to higher scales can violate positivity of DPDs in colour space, both when one evolves in the renormalisation scale of one parton and when one evolves in all scales simultaneously.

DPDs from parton splitting at two-loop accuracy
In section 2 we saw that the perturbative splitting mechanism gives DPDs that satisfy positivity if the splitting is computed at LO, i.e. at one-loop accuracy. This is not surprising, since the LO splitting formula in the s channel basis can be written as a squared matrix element (in the mixed representation of definite plus-momentum and transverse position for the observed partons). Starting from two loops, the splitting formula has explicit logarithms of the renormalisation scale µ and of the rapidity parameter ζ p , which respectively result from subtractions for ultraviolet and rapidity divergences. As a consequence of these subtractions, positivity is no longer guaranteed. It is then natural to ask whether the resulting DPDs at small y violate positivity in colour space, and if so, by how much. We address this question in the present section, using the results of the recent two-loop calculation in [20]. The generalisation of the LO splitting formula (9) to higher orders has the form and the special convolution Here we write and recall that according to (11). The analogue of (44) for distributions in the s channel basis is readily obtained using the transformations in section 2. We limit our attention to the quark-antiquark sector and consider the distributions Among these, only F 88 uū is nonzero at order a s , whilst all others start at order a 2 s .
Parton combinations appearing first at two loops. At order a 2 s the distributions F ud and F ud receive a contribution only from the kernels V qq ,q or V qq ,q , which respectively correspond to the graphs in figure 5(a) and 5(d) and further graphs with identical topology. As a consequence, distributions for different colour representations are proportional to each other, as specified in equation (4.40) of [20]. This leads to the relations in (49). Furthermore, the splitting kernels for F ud and F ud are proportional to each other, because 11 V qq ,q = 11 V qq ,q at two-loop accuracy. Using the basis transform (5) and the proportionality factors between t channel singlet and octet kernels in equation (4.35) of [20], one obtains with V 33 qq ,q = 2 9
Here the convolution V (z 2 , z 1 ) ⊗ f (z) is defined as in (46) with u and 1 − u interchanged on the r.h.s. In the last term of (51) we used the relation Vq q ,q = V qq ,q from charge conjugation invariance. Notice that the kernel 11 V qq ,q includes an overall factor a 2 s .
The distributions F 33 uu and F 66 uu receive contributions from the kernels V qq ,q and V v qq,q with different weights and are therefore not proportional to each other. Using equation (3.3) in [20], we find that the relevant kernels read and are to be convolved with f u (z). At order a 2 s , the kernel 11 V qq ,q and hence all kernels in (52) and (53) require ultraviolet renormalisation but no subtraction of rapidity divergences. As a consequence, they depend linearly on the renormalisation group logarithm L y . A natural choice of scale in the fixed-order formula (44) is µ = µ y , so that L y = 0. In a numerical study of the two-gluon distributions RR F gg in [20], we found that the size of a 2 s corrections relative to the a s term is moderate for µ = µ y but grows substantially for µ = µ y /2 or µ = 2µ y . In the present work, we therefore consider a smaller amount of variation and take µ = 1.2µ y or µ = µ y /1.2 as alternative scales, which corresponds to L y ≈ ±0. 36. We checked that the size of a 2 s corrections in the two-gluon sector remains moderate for these values. The DPDs depend on µ as specified by the relevant DGLAP equations. When comparing the fixed-order formula (44) for different µ, one thus sees evolution effects truncated at the lowest order in a s . In figures 6 and 7, we plot the kernels 11 V qq ,q , V 33 qq,q , and V 66 qq,q . We see that they are negative over wide ranges of z and u, although they are obtained from a sum of graphs that correspond to the squared amplitude for q → qq q or q → qqq. This illustrates that the subtraction of ultraviolet divergences can indeed lead to a negative result. We note that the kernels go to large positive values for z → 1, which can be traced back to terms diverging like log 1/(1 − z) in that limit. Finally, we observe that the kernels increase with µ for all z and u. We now investigate to which extent the negative regions in the splitting kernels lead to negative DPDs. To this end, we evaluate the splitting formula for the distributions in (49) and (50) with the central PDFs from the CT14nlo set, having checked that the PDFs are positive in the kinematics of interest. As we did in section 4, we fix y such that µ y = 10 GeV. In figure 8 we show F 33 ud and F 11 ud for µ = µ y /1.2 and µ = µ y . The curves are scaled such that the difference between them originates from the different PDFs in (51) and not from the different normalisation of the kernels in (52). We observe a strong effect of scale evolution: for x 1 ∼ x 2 the distributions are negative at the smaller scale but positive at µ = µ y (and also at µ = 1.2µ y , which is not shown in the figure). We also note that the negative values in figure 8(c) are tiny compared with the size of the same distributions at other values of x 1 for the same x 2 . The distributions F 33 uu and F 66 uu are shown in figure 9, with a scaling factor such that the difference between the curves is due to the contribution of 11 V v qq,q to the kernels in (53). The situation for µ = µ y /1.2 (not shown in the figure) is qualitatively similar to the one at µ = µ y , where we find negative values at x 1 ∼ x 2 for F 33 uu but not for F 66 uu . At the higher scale µ = 1.2µ y , all values are positive. As in figure 8(c), the negative values in figure 9(c) are tiny compared with the size of the distribution at other momentum fractions. In this sense, the violations of positivity we have shown so far may be regarded as minor.
Quark-antiquark distributions. The last two distributions in (50) are for a quark and an antiquark of equal flavour. F 11 uū receives contributions from all three kernels in the bottom row of figure 5 and from the real two-loop graphs for the splitting g → qq, an example of which is shown in figure 2(b). F 88 uū receives contributions from the same graphs, from the LO graph in figure 2(a), and from virtual two-loop graphs such as the one in figure 2(c). The      uū at x 1 = x 2 at two-loop level, evaluated with x 1 x 2 ζ p = µ 2 for two values of µ. The curves labelled "g" are for the splitting g → qq, and the curve labelled "q" is for the sum of all splitting contributions initiated by a quark or an antiquark. The labels "regular", "[1 − z] + ", and "δ(1 − x)" refer to the different parts of the kernels V 11 [2,0] + L y V 11 [2,1] discussed below equation (55).
virtual graphs depend on the colour of the observed partons in the same way as the LO graph and hence do not contribute to F 11 uū . The distributions F 11 uū and F 88 uū require both ultraviolet renormalisation and the subtraction of rapidity divergences. The latter appears in the splitting process g → qq, whose kernels have a more complicated structure than the ones considered so far. Up to order a 2 s , they can be written as qq,g (u) = u 2 + (1 − u) 2 2 is the LO splitting kernel appearing in (9) and contains the double logarithms associated with rapidity divergences. We note that (55) is obtained with the standard definition of the MS scheme, and that the term π 2 /6 is absent if one instead uses the definition proposed by Collins in section 3.2.6 of [22]. The coefficients V [2,0] and V [2,1] in (54) are smooth functions of u but distributions in z. They consist of a regular part, a part proportional to the plus-distribution 1/[1 − z] + , and a term proportional to δ(1 − z). The regular part is a smooth function of z but may have a log(1 − z) singularity for z → 1, similar to what we saw for the pure quark kernels in figures 6 and 7.
In figure 10 we show the different contributions to F 11 uū for x 1 = x 2 . The regular part of the g → qq kernel turns out to be negative for all values of z and u and results in a large negative contribution to the DPD. At µ = 1.2µ y , another negative contribution comes from the part of the kernel that goes with the plus-distribution 1/[1 − z] + . Adding up all contributions, one obtains the distribution shown in the two upper rows of figure 11. The negative contributions dominate for x 1 = x 2 up to a few 0.01, whereas for larger momentum fractions the positive contributions from quark or antiquark splitting gradually takes over. In the two lower rows of figure 11, we see that the regions of negative F 11 uū are centred around x 1 ∼ x 2 , with values that are small compared with the size of the distribution at other values of x 1 at the same x 2 . This is the same phenomenon that we observed earlier for F ud , F ud , and F 33 uu . However, in the present case the negative values around x 1 ∼ x 2 decrease with µ, so that the violation of positivity in F 11 uū becomes more pronounced as µ becomes larger. In figure 11 we also see that a mild variation of the rapidity parameter ζ p has a rather small effect in the kinematics considered here. The choice of ζ p is only relevant if µ = µ y , because L ζ is multiplied by L y in (55). Let us finally turn to the distribution F 88 uū . Evaluated at one-loop accuracy, this distribution is positive, so that negative values at two-loop level can only appear when the contribution of order a 2 s is larger in size than the one of order a s . In such a case, one may of course worry whether the unknown contributions of yet higher orders will change the sign of the distribution again. This situation is qualitatively different from the one for the distributions discussed so far, where the terms of order a 2 s give the first non-vanishing contribution. In figures 12 and 13 we show F 88 uū for different values of the momentum fractions. We always take x 1 x 2 ζ p = µ 2 , bearing in mind that according to (54) the effect of varying ζ p is eight times smaller for F 88 uū than it is for F 11 uū . As can be seen in figure 12, the difference between the distributions at µ = µ y to µ = 1.2µ y is rather small, in contrast to what we found for the other distributions discussed so far. This is not surprising, because for F 88 uū the relative effect of changing the scale from µ a to µ b is of order a s log(µ a /µ b ), whereas it is of order log(µ a /µ b ) for distributions that receive their first nonzero contribution at two loops. In all panels of figures 12 and 13, we show the result obtained with either the one-loop kernel or with the sum of one-and two-loop kernels. In both cases we take the same NLO PDFs, so that the difference between the LO and the LO+NLO curves directly shows the impact of the two-loop kernels. We find that the size of the a 2 s corrections is often moderate but becomes large in several kinematic situations.
1. As discussed in section 4.3 of [20], the two-loop corrections for g → qq and q → q q are enhanced at small x 1 + x 2 . This is seen in the left panels of figures 12(a), 12(b), 13(b), and 13(c). As follows from equations (4.46) and (4.48) in [20], the enhanced corrections provide a positive contribution to F 88 uū . An all-order resummation of the enhanced corrections using techniques from small-x factorisation may be possible, but details of this have not been worked out.
2. As explained in section 4.4 of [20], the splitting graph for u → uū in figure 5(d) leads to a behaviour of the DPD like 1/(1 − u) ≈ x 1 /x 2 for x 2 x 1 , which is absent at order a s . This explains why in figure 13(a) the LO+NLO result for the scaled DPD x 1 x 2 F 88 uū goes to a finite value when x 2 x 1 , whereas the LO result goes to zero. It also explains the huge relative NLO corrections seen in the right panel of figure 13(b). In the latter case, the splitting process u → uū is further enhanced by the fact that the u quark distribution becomes the dominant PDF with increasing x.
The appearance of an additional power 1/(1 − u) is unique for the step from LO to NLO in this channel and will not repeat itself at yet higher orders.

For x 2
x 1 we find a large negative two-loop contribution to F 88 uū , which is seen in the right panel of figure 13(c) and more clearly in figure 13(d), where the LO+NLO result becomes just slightly negative. This can be traced back to a negative term − 5a 2 qq,g (u) log 2 u (56)  in the kernel V 88 qq,g , which is enhanced by two powers of log u compared with the LO expression.
It would require further analysis to understand whether this type of enhancement repeats itself at yet higher orders and, if so, whether it can be resummed to all orders. We therefore cannot say whether the negative values seen in figure 13(d) will disappear when higher order terms are included.
Notice that figures 13(c) and figure 13(d) refer to the lowest scale µ = µ y /1.2 considered in this study. The corresponding plots for µ = µ y or 1.2µ y show large negative NLO corrections as well, but the LO+NLO curves do not become negative any more.
We remark in passing that the enhanced two-loop contributions described in points 1 and 3 scale like the LO kernel for g → qq and hence cancel out in F 11 uū . By contrast, the 1/(1 − u) behaviour discussed in point 2 appears both in F 88 uū and F 11 uū .
Colour summed distributions Let us finally consider the colour summed distributions 11 F ud , 11 F ud , 11 F uu , and 11 F uū . The relations (3) and (49) imply that 11 F ud and 11 F ud are negative whenever their s channel counterparts are. By contrast, we find that 11 F uu and 11 F uū remain positive for the kinematic settings shown in figures 9 to 13 .

Summary
In the context of the parton model, DPDs in the s channel colour basis are probability densities for finding two partons in a definite colour state (with specified longitudinal momenta and specified transverse distance from each other). This leads to the expectation that s channel DPDs for unpolarised partons should be non-negative. If it is satisfied, this positivity property provides valuable constraints on the colour dependence of DPDs, along with a strategy to model them by saturating the positivity bounds at a certain scale.
In the present work, we investigate whether evolution of DPDs to higher scales preserves the positivity property under the assumption that it holds at the starting scale. We limit ourselves to unpolarised partons, and we exclude two-gluon DPDs from our consideration because their colour structure is much more involved than the one of pure quark or quarkgluon distributions. DPDs in the s channel colour basis are subject to Collins-Soper evolution in the rapidity parameter ζ p , and this evolution includes mixing between different colour channels. Provided that the renormalisation scales µ 1 and µ 2 associated with the two partons are large enough compared with 1/y, the Collins-Soper kernel 8 J(y, µ 1 , µ 2 ) for DPDs is negative. Under this condition, evolution to higher ζ p preserves positivity. Conversely, backward evolution to sufficiently small ζ p eventually leads to negative distributions (except for the special case in which all t channel distributions other than 11 F a 1 a 2 are zero and hence all s channel distributions are independent of ζ p ). We next consider the DGLAP equations for evolution in one of the scales µ 1 or µ 2 , with the evolution kernels taken at LO. We find that evolution to higher scales is not guaranteed to preserve positivity: there are initial conditions that satisfy positivity but lead to negative s channel distributions at higher scales. This is due to the convolution of a plus distribution in the evolution kernel with a DPD different from the one being evolved. There is no contribution of this type in the LO evolution equations for polarised colour summed DPDs 11 F a 1 a 2 , which conserve positivity in the same way as the LO evolution of polarised PDFs [4]. In a numerical illustration, we choose initial conditions where certain s channel DPDs are zero and see that they turn negative at slightly higher scales. We also study joint DGLAP and Collins-Soper evolution in the common scale µ = µ 1 = µ 2 = x 1 x 2 ζ p and find that positivity is not preserved, for the same reasons as above. At small inter-parton distance y, the initial conditions for DPD evolution can be computed using the perturbative splitting mechanism, setting µ = µ 1 = µ 2 ≈ 1/y and using a fixedorder truncation of the DPD splitting kernels. It is easy to see that at order a s one obtains DPDs that satisfy positivity in colour space. At order a 2 s this no longer holds: in a numerical study we obtain negative values for colour channels in which the distributions are zero at order a s , and also for the distribution F 88 uū , which is nonzero at order a s . Negative values are also found for the colour summed distributions 11 F ud and 11 F ud . In several cases, the considered distributions have no ζ p dependence at order a 2 s , and negative values of them can be uniquely traced back to the subtraction of ultraviolet divergences implied in the definition of twist-two operators. The explicit form of this subtraction is given in section 2.6 of [20]. The negative values we find for the distributions are small compared with the size of the same distributions at other values of the momentum fractions x 1 and x 2 . In this sense, the violations of positivity we have seen may be regarded as "relatively small". In view of this, we should also caution that negative values obtained with the splitting formula at order a 2 s may turn into positive ones when yet higher orders are included. Note that the violations of positivity just discussed refer to DPDs defined with MS renormalisation of twist-two operators. By contrast, the violation of positivity by forward DGLAP evolution described earlier occurs at LO and is hence not specific to the MS scheme. We conclude that the positivity of DPDs in full colour space cannot be taken for granted and can be violated in physically realistic settings. Using positivity as a guide for modelling DPDs at large y may still be an option when there is a lack of better information. It should however be done with due caution, and one should check whether the chosen initial conditions give positive distributions when evolved to higher scales.

A Colour space projectors
In this appendix, we list the colour space projectors that appear in the definitions (1) and (2) of DPDs in the s and t channel bases. In the pure quark and the mixed quark-gluon sector, we have and P ij i j For completeness, we also give the projectors for the pure gluon sector, although they are not used in the present work.
In all cases, we have set the number of colours to N = 3. Further projectors are obtained by exchanging the representation labels and the corresponding indices: P r 1 r 2 r 1 r 2 RR = P r 1 r 2 r 1 r 2 R R .
One readily verifies the completeness relations R P r 1 r 2 r 2 r 1 RR = δ r 1 r 1 δ r 2 r 2 for P RR in (57) , R P r 1 r 2 r 1 r 2 RR = δ r 1 r 1 δ r 2 r 2 for P RR in (59), (60), (61) , where in each case the sum runs over all available representations. The multiplicity of a representation R can be computed from the trace m(R) =    P r 1 r 2 r 2 r 1 RR for P RR in (57), P r 1 r 2 r 1 r 2